Search

Understanding the nature of dark matter in the Universe is an important goal of modern cosmology. A key method for probing this distribution is via weak gravitational lensing mass-mapping - a challenging ill-posed inverse problem where one infers the convergence field from observed shear measurements. Upcoming stage IV surveys, such as those made by the Vera C. Rubin Observatory and Euclid satellite, will provide a greater quantity and precision of data for lensing analyses, necessitating high-fidelity mass-mapping methods that are computationally efficient and that also provide uncertainties for integration into downstream cosmological analyses. In this work we introduce MMGAN, a novel mass-mapping method based on a regularised conditional generative adversarial network (GAN) framework, which generates approximate posterior samples of the convergence field given shear data. We adopt Wasserstein GANs to improve training stability and apply regularisation techniques to overcome mode collapse, issues that otherwise are particularly acute for conditional GANs. We train and validate our model on a mock COSMOS-style dataset before applying it to true COSMOS survey data. Our approach significantly outperforms the Kaiser-Squires technique and achieves similar reconstruction fidelity as alternative state-of-the-art deep learning approaches. Notably, while alternative approaches for generating samples from a learned posterior are slow (e.g. requiring $\sim$10 GPU minutes per posterior sample), MMGAN can produce a high-quality convergence sample in less than a second.

We present an accelerated pipeline, based on high-performance computing techniques and normalizing flows, for joint Bayesian parameter estimation and model selection and demonstrate its efficiency in gravitational wave astrophysics. We integrate the Jim inference toolkit, a normalizing flow-enhanced Markov chain Monte Carlo (MCMC) sampler, with the learned harmonic mean estimator. Our Bayesian evidence estimates run on $1$ GPU are consistent with traditional nested sampling techniques run on $16$ CPU cores, while reducing the computation time by factors of $5\times$ and $15\times$ for $4$-dimensional and $11$-dimensional gravitational wave inference problems, respectively. Our code is available in well-tested and thoroughly documented open-source packages, ensuring accessibility and reproducibility for the wider research community., Comment: accepted to NeurIPS 2024 workshop on Machine Learning and the Physical Sciences

We demonstrate the successful use of scattering representations without further compression for simulation-based inference (SBI) with images (i.e. field-level), illustrated with a cosmological case study. Scattering representations provide a highly effective representational space for subsequent learning tasks, although the higher dimensional compressed space introduces challenges. We overcome these through spatial averaging, coupled with more expressive density estimators. Compared to alternative methods, such an approach does not require additional simulations for either training or computing derivatives, is interpretable, and resilient to covariate shift. As expected, we show that a scattering only approach extracts more information than traditional second order summary statistics., Comment: 9 pages, 2 figures, accepted by NeurIPS workshop on Machine Learning and the Physical Sciences

Scattering transforms are a new type of summary statistics recently developed for the study of highly non-Gaussian processes, which have been shown to be very promising for astrophysical studies. In particular, they allow one to build generative models of complex non-linear fields from a limited amount of data. In the context of upcoming cosmological surveys, the extension of these tools to spherical data is necessary. We develop scattering transforms on the sphere and focus on the construction of maximum-entropy generative models of astrophysical fields. The quality of the generative models, both statistically and visually, is very satisfying, which therefore open up a wide range of new applications for future cosmological studies., Comment: Contribution to the 2024 Cosmology session of the 58th Rencontres de Moriond. For details, please refer to the full article arXiv:2407.07007

Scattering transforms are a new type of summary statistics recently developed for the study of highly non-Gaussian processes, which have been shown to be very promising for astrophysical studies. In particular, they allow one to build generative models of complex non-linear fields from a limited amount of data, and have also been used as the basis of new statistical component separation algorithms. In the context of upcoming cosmological surveys, such as LiteBIRD for the cosmic microwave background polarization or Rubin-LSST and Euclid for study of the large scale structures of the Universe, the extension of these tools to spherical data is necessary. We develop scattering transforms on the sphere and focus on the construction of maximum-entropy generative models of several astrophysical fields. We construct, from a single target field, generative models of homogeneous astrophysical and cosmological fields, whose samples are quantitatively compared to the target fields using common statistics (power spectrum, pixel probability density function and Minkowski functionals). Our sampled fields agree well with the target fields, both statistically and visually. These generative models therefore open up a wide range of new applications for future astrophysical and cosmological studies; particularly those for which very little simulated data is available. We make our code available to the community so that this work can be easily reproduced and developed further., Comment: submitted to A&A, 13 pages, 10 figures, code available at https://github.com/astro-informatics/s2scat

Understanding the large-scale structure of the Universe and unravelling the mysteries of dark matter are fundamental challenges in contemporary cosmology. Reconstruction of the cosmological matter distribution from lensing observables, referred to as 'mass-mapping' is an important aspect of this quest. Mass-mapping is an ill-posed problem, meaning there is inherent uncertainty in any convergence map reconstruction. The demand for fast and efficient reconstruction techniques is rising as we prepare for upcoming surveys. We present a novel approach which utilises deep learning, in particular a conditional Generative Adversarial Network (cGAN), to approximate samples from a Bayesian posterior distribution, meaning they can be interpreted in a statistically robust manner. By combining data-driven priors with recent regularisation techniques, we introduce an approach that facilitates the swift generation of high-fidelity, mass maps. Furthermore, to validate the effectiveness of our approach, we train the model on mock COSMOS-style data, generated using Colombia Lensing's kappaTNG mock weak lensing suite. These preliminary results showcase compelling convergence map reconstructions and ongoing refinement efforts are underway to enhance the robustness of our method further., Comment: 2 pages, 1 figure, submitted for conference proceedings for '58th Recontres de Moriond', 2024

We advocate for a new paradigm of cosmological likelihood-based inference, leveraging recent developments in machine learning and its underlying technology, to accelerate Bayesian inference in high-dimensional settings. Specifically, we combine (i) emulation, where a machine learning model is trained to mimic cosmological observables, e.g. CosmoPower-JAX; (ii) differentiable and probabilistic programming, e.g. JAX and NumPyro, respectively; (iii) scalable Markov chain Monte Carlo (MCMC) sampling techniques that exploit gradients, e.g. Hamiltonian Monte Carlo; and (iv) decoupled and scalable Bayesian model selection techniques that compute the Bayesian evidence purely from posterior samples, e.g. the learned harmonic mean implemented in harmonic. This paradigm allows us to carry out a complete Bayesian analysis, including both parameter estimation and model selection, in a fraction of the time of traditional approaches. First, we demonstrate the application of this paradigm on a simulated cosmic shear analysis for a Stage IV survey in 37- and 39-dimensional parameter spaces, comparing $\Lambda$CDM and a dynamical dark energy model ($w_0w_a$CDM). We recover posterior contours and evidence estimates that are in excellent agreement with those computed by the traditional nested sampling approach while reducing the computational cost from 8 months on 48 CPU cores to 2 days on 12 GPUs. Second, we consider a joint analysis between three simulated next-generation surveys, each performing a 3x2pt analysis, resulting in 157- and 159-dimensional parameter spaces. Standard nested sampling techniques are simply unlikely to be feasible in this high-dimensional setting, requiring a projected 12 years of compute time on 48 CPU cores; on the other hand, the proposed approach only requires 8 days of compute time on 24 GPUs. All packages used in our analyses are publicly available., Comment: 14 pages, 6 figures. Accepted for publication in the Open Journal of Astrophysics. Codes available at https://github.com/alessiospuriomancini/cosmopower, https://github.com/dpiras/cosmopower-jax, https://github.com/astro-informatics/harmonic/

With the next generation of interferometric telescopes, such as the Square Kilometre Array (SKA), the need for highly computationally efficient reconstruction techniques is particularly acute. The challenge in designing learned, data-driven reconstruction techniques for radio interferometry is that they need to be agnostic to the varying visibility coverages of the telescope, since these are different for each observation. Because of this, learned post-processing or learned unrolled iterative reconstruction methods must typically be retrained for each specific observation, amounting to a large computational overhead. In this work we develop learned post-processing and unrolled iterative methods for varying visibility coverages, proposing training strategies to make these methods agnostic to variations in visibility coverage with minimal to no fine-tuning. Learned post-processing techniques are heavily dependent on the prior information encoded in training data and generalise poorly to other visibility coverages. In contrast, unrolled iterative methods, which include the telescope measurement operator inside the network, achieve state-of-the-art reconstruction quality and computation time, generalising well to other coverages and require little to no fine-tuning. Furthermore, they generalise well to realistic radio observations and are able to reconstruct the high dynamic range of these images.

We present the learned harmonic mean estimator with normalizing flows - a robust, scalable and flexible estimator of the Bayesian evidence for model comparison. Since the estimator is agnostic to sampling strategy and simply requires posterior samples, it can be applied to compute the evidence using any Markov chain Monte Carlo (MCMC) sampling technique, including saved down MCMC chains, or any variational inference approach. The learned harmonic mean estimator was recently introduced, where machine learning techniques were developed to learn a suitable internal importance sampling target distribution to solve the issue of exploding variance of the original harmonic mean estimator. In this article we present the use of normalizing flows as the internal machine learning technique within the learned harmonic mean estimator. Normalizing flows can be elegantly coupled with the learned harmonic mean to provide an approach that is more robust, flexible and scalable than the machine learning models considered previously. We perform a series of numerical experiments, applying our method to benchmark problems and to a cosmological example in up to 21 dimensions. We find the learned harmonic mean estimator is in agreement with ground truth values and nested sampling estimates. The open-source harmonic Python package implementing the learned harmonic mean, now with normalizing flows included, is publicly available., Comment: 14 pages, 8 figures, harmonic code available at https://github.com/astro-informatics/harmonic

Directional wavelet dictionaries are hierarchical representations which efficiently capture and segment information across scale, location and orientation. Such representations demonstrate a particular affinity to physical signals, which often exhibit highly anisotropic, localised multiscale structure. Many physically important signals are observed over spherical domains, such as the celestial sky in cosmology. Leveraging recent advances in computational harmonic analysis, we design new highly distributable and automatically differentiable directional wavelet transforms on the $2$-dimensional sphere $\mathbb{S}^2$ and $3$-dimensional ball $\mathbb{B}^3 = \mathbb{R}^+ \times \mathbb{S}^2$ (the space formed by augmenting the sphere with the radial half-line). We observe up to a $300$-fold and $21800$-fold acceleration for signals on the sphere and ball, respectively, compared to existing software, whilst maintaining 64-bit machine precision. Not only do these algorithms dramatically accelerate existing spherical wavelet transforms, the gradient information afforded by automatic differentiation unlocks many data-driven analysis techniques previously not possible for these spaces. We publicly release both S2WAV and S2BALL, open-sourced JAX libraries for our transforms that are automatically differentiable and readily deployable both on and over clusters of hardware accelerators (e.g. GPUs & TPUs)., Comment: code available on the sphere at https://github.com/astro-informatics/s2wav and on the ball at https://github.com/astro-informatics/s2ball

Next-generation radio interferometers like the Square Kilometer Array have the potential to unlock scientific discoveries thanks to their unprecedented angular resolution and sensitivity. One key to unlocking their potential resides in handling the deluge and complexity of incoming data. This challenge requires building radio interferometric imaging methods that can cope with the massive data sizes and provide high-quality image reconstructions with uncertainty quantification (UQ). This work proposes a method coined QuantifAI to address UQ in radio-interferometric imaging with data-driven (learned) priors for high-dimensional settings. Our model, rooted in the Bayesian framework, uses a physically motivated model for the likelihood. The model exploits a data-driven convex prior, which can encode complex information learned implicitly from simulations and guarantee the log-concavity of the posterior. We leverage probability concentration phenomena of high-dimensional log-concave posteriors that let us obtain information about the posterior, avoiding MCMC sampling techniques. We rely on convex optimisation methods to compute the MAP estimation, which is known to be faster and better scale with dimension than MCMC sampling strategies. Our method allows us to compute local credible intervals, i.e., Bayesian error bars, and perform hypothesis testing of structure on the reconstructed image. In addition, we propose a novel blazing-fast method to compute pixel-wise uncertainties at different scales. We demonstrate our method by reconstructing radio-interferometric images in a simulated setting and carrying out fast and scalable UQ, which we validate with MCMC sampling. Our method shows an improved image quality and more meaningful uncertainties than the benchmark method based on a sparsity-promoting prior. QuantifAI's source code: https://github.com/astro-informatics/QuantifAI., Comment: 30 pages, 14 figures, 10 tables, code available at https://github.com/astro-informatics/QuantifAI

Many areas of science and engineering encounter data defined on spherical manifolds. Modelling and analysis of spherical data often necessitates spherical harmonic transforms, at high degrees, and increasingly requires efficient computation of gradients for machine learning or other differentiable programming tasks. We develop novel algorithmic structures for accelerated and differentiable computation of generalised Fourier transforms on the sphere $\mathbb{S}^2$ and rotation group $\text{SO}(3)$, i.e. spherical harmonic and Wigner transforms, respectively. We present a recursive algorithm for the calculation of Wigner $d$-functions that is both stable to high harmonic degrees and extremely parallelisable. By tightly coupling this with separable spherical transforms, we obtain algorithms that exhibit an extremely parallelisable structure that is well-suited for the high throughput computing of modern hardware accelerators (e.g. GPUs). We also develop a hybrid automatic and manual differentiation approach so that gradients can be computed efficiently. Our algorithms are implemented within the JAX differentiable programming framework in the S2FFT software code. Numerous samplings of the sphere are supported, including equiangular and HEALPix sampling. Computational errors are at the order of machine precision for spherical samplings that admit a sampling theorem. When benchmarked against alternative C codes we observe up to a 400-fold acceleration. Furthermore, when distributing over multiple GPUs we achieve very close to optimal linear scaling with increasing number of GPUs due to the highly parallelised and balanced nature of our algorithms. Provided access to sufficiently many GPUs our transforms thus exhibit an unprecedented effective linear time complexity., Comment: 30 pages, 7 figures, accepted by Journal of Computational Physics, code available at https://github.com/astro-informatics/s2fft

Cosmic strings are linear topological defects that may have been produced during symmetry-breaking phase transitions in the very early Universe. In an expanding Universe the existence of causally separate regions prevents such symmetries from being broken uniformly, with a network of cosmic string inevitably forming as a result. To faithfully generate observables of such processes requires computationally expensive numerical simulations, which prohibits many types of analyses. We propose a technique to instead rapidly emulate observables, thus circumventing simulation. Emulation is a form of generative modelling, often built upon a machine learning backbone. End-to-end emulation often fails due to high dimensionality and insufficient training data. Consequently, it is common to instead emulate a latent representation from which observables may readily be synthesised. Wavelet phase harmonics are an excellent latent representations for cosmological fields, both as a summary statistic and for emulation, since they do not require training and are highly sensitive to non-Gaussian information. Leveraging wavelet phase harmonics as a latent representation, we develop techniques to emulate string induced CMB anisotropies over a 7.2 degree field of view, with sub-arcminute resolution, in under a minute on a single GPU. Beyond generating high fidelity emulations, we provide a technique to ensure these observables are distributed correctly, providing a more representative ensemble of samples. The statistics of our emulations are commensurate with those calculated on comprehensive Nambu-Goto simulations. Our findings indicate these fast emulation approaches may be suitable for wide use in, e.g., simulation based inference pipelines. We make our code available to the community so that researchers may rapidly emulate cosmic string induced CMB anisotropies for their own analysis., Comment: code available at https://github.com/astro-informatics/stringgen

Proximal nested sampling was introduced recently to open up Bayesian model selection for high-dimensional problems such as computational imaging. The framework is suitable for models with a log-convex likelihood, which are ubiquitous in the imaging sciences. The purpose of this article is two-fold. First, we review proximal nested sampling in a pedagogical manner in an attempt to elucidate the framework for physical scientists. Second, we show how proximal nested sampling can be extended in an empirical Bayes setting to support data-driven priors, such as deep neural networks learned from training data., Comment: 9 pages, 4 figures

Computing the marginal likelihood (also called the Bayesian model evidence) is an important task in Bayesian model selection, providing a principled quantitative way to compare models. The learned harmonic mean estimator solves the exploding variance problem of the original harmonic mean estimation of the marginal likelihood. The learned harmonic mean estimator learns an importance sampling target distribution that approximates the optimal distribution. While the approximation need not be highly accurate, it is critical that the probability mass of the learned distribution is contained within the posterior in order to avoid the exploding variance problem. In previous work a bespoke optimization problem is introduced when training models in order to ensure this property is satisfied. In the current article we introduce the use of normalizing flows to represent the importance sampling target distribution. A flow-based model is trained on samples from the posterior by maximum likelihood estimation. Then, the probability density of the flow is concentrated by lowering the variance of the base distribution, i.e. by lowering its "temperature", ensuring its probability mass is contained within the posterior. This approach avoids the need for a bespoke optimisation problem and careful fine tuning of parameters, resulting in a more robust method. Moreover, the use of normalizing flows has the potential to scale to high dimensional settings. We present preliminary experiments demonstrating the effectiveness of the use of flows for the learned harmonic mean estimator. The harmonic code implementing the learned harmonic mean, which is publicly available, has been updated to now support normalizing flows., Comment: 9 pages, 6 figures. arXiv admin note: text overlap with arXiv:2111.12720

A new golden age in astronomy is upon us, dominated by data. Large astronomical surveys are broadcasting unprecedented rates of information, demanding machine learning as a critical component in modern scientific pipelines to handle the deluge of data. The upcoming Legacy Survey of Space and Time (LSST) of the Vera C. Rubin Observatory will raise the big-data bar for time-domain astronomy, with an expected 10 million alerts per-night, and generating many petabytes of data over the lifetime of the survey. Fast and efficient classification algorithms that can operate in real-time, yet robustly and accurately, are needed for time-critical events where additional resources can be sought for follow-up analyses. In order to handle such data, state-of-the-art deep learning architectures coupled with tools that leverage modern hardware accelerators are essential. We showcase how the use of modern deep compression methods can achieve a $18\times$ reduction in model size, whilst preserving classification performance. We also show that in addition to the deep compression techniques, careful choice of file formats can improve inference latency, and thereby throughput of alerts, on the order of $8\times$ for local processing, and $5\times$ in a live production setting. To test this in a live setting, we deploy this optimised version of the original time-series transformer, t2, into the community alert broking system of FINK on real Zwicky Transient Facility (ZTF) alert data, and compare throughput performance with other science modules that exist in FINK. The results shown herein emphasise the time-series transformer's suitability for real-time classification at LSST scale, and beyond, and introduce deep model compression as a fundamental tool for improving deploy-ability and scalable inference of deep learning models for transient classification., Comment: 16 pages, 11 figures

Inspired by recent interest in geometric deep learning, this work generalises the recently developed Slepian scale-discretised wavelets on the sphere to Riemannian manifolds. Through the sifting convolution, one may define translations and, thus, convolutions on manifolds - which are otherwise not well-defined in general. Slepian wavelets are constructed on a region of a manifold and are therefore suited to problems where data only exists in a particular region. The Slepian functions, on which Slepian wavelets are built, are the basis functions of the Slepian spatial-spectral concentration problem on the manifold. A tiling of the Slepian harmonic line with smoothly decreasing generating functions defines the scale-discretised wavelets; allowing one to probe spatially localised, scale-dependent features of a signal. By discretising manifolds as graphs, the Slepian functions and wavelets of a triangular mesh are presented. Through a wavelet transform, the wavelet coefficients of a field defined on the mesh are found and used in a straightforward thresholding denoising scheme., Comment: 12 pages, 12 figures

The Segmented Planar Imaging Detector for Electro-Optical Reconnaissance (SPIDER) is an optical interferometric imaging device that aims to offer an alternative to the large space telescope designs of today with reduced size, weight and power consumption. This is achieved through interferometric imaging. State-of-the-art methods for reconstructing images from interferometric measurements adopt proximal optimization techniques, which are computationally expensive and require handcrafted priors. In this work we present two data-driven approaches for reconstructing images from measurements made by the SPIDER instrument. These approaches use deep learning to learn prior information from training data, increasing the reconstruction quality, and significantly reducing the computation time required to recover images by orders of magnitude. Reconstruction time is reduced to ${\sim} 10$ milliseconds, opening up the possibility of real-time imaging with SPIDER for the first time. Furthermore, we show that these methods can also be applied in domains where training data is scarce, such as astronomical imaging, by leveraging transfer learning from domains where plenty of training data are available., Comment: 21 pages, 14 figures

Uncertainty quantification is a crucial step of cosmological mass-mapping that is often ignored. Suggested methods are typically only approximate or make strong assumptions of Gaussianity of the shear field. Probabilistic sampling methods, such as Markov chain Monte Carlo (MCMC), draw samples form a probability distribution, allowing for full and flexible uncertainty quantification, however these methods are notoriously slow and struggle in the high-dimensional parameter spaces of imaging problems. In this work we use, for the first time, a trans-dimensional MCMC sampler for mass-mapping, promoting sparsity in a wavelet basis. This sampler gradually grows the parameter space as required by the data, exploiting the extremely sparse nature of mass maps in wavelet space. The wavelet coefficients are arranged in a tree-like structure, which adds finer scale detail as the parameter space grows. We demonstrate the trans-dimensional sampler on galaxy cluster-scale images where the planar modelling approximation is valid. In high-resolution experiments, this method produces naturally parsimonious solutions, requiring less than 1% of the potential maximum number of wavelet coefficients and still producing a good fit to the observed data. In the presence of noisy data, trans-dimensional MCMC produces a better reconstruction of mass-maps than the standard smoothed Kaiser-Squires method, with the addition that uncertainties are fully quantified. This opens up the possibility for new mass maps and inferences about the nature of dark matter using the new high-resolution data from upcoming weak lensing surveys such as Euclid.

The Vera C. Rubin Observatory's Legacy Survey of Space and Time (LSST) will discover an unprecedented number of supernovae (SNe), making spectroscopic classification for all the events infeasible. LSST will thus rely on photometric classification, whose accuracy depends on the not-yet-finalized LSST observing strategy. In this work, we analyze the impact of cadence choices on classification performance using simulated multi-band light curves. First, we simulate SNe with an LSST baseline cadence, a non-rolling cadence, and a presto-color cadence which observes each sky location three times per night instead of twice. Each simulated dataset includes a spectroscopically-confirmed training set, which we augment to be representative of the test set as part of the classification pipeline. Then, we use the photometric transient classification library snmachine to build classifiers. We find that the active region of the rolling cadence used in the baseline observing strategy yields a 25% improvement in classification performance relative to the background region. This improvement in performance in the actively-rolling region is also associated with an increase of up to a factor of 2.7 in the number of cosmologically-useful Type Ia supernovae relative to the background region. However, adding a third visit per night as implemented in presto-color degrades classification performance due to more irregularly sampled light curves. Overall, our results establish desiderata on the observing cadence related to classification of full SNe light curves, which in turn impacts photometric SNe cosmology with LSST., Comment: 22 pages, 14 figures. Changed to match version accepted by the Astrophysical Journal Supplement Series (accepted 06/02/2023)

No existing spherical convolutional neural network (CNN) framework is both computationally scalable and rotationally equivariant. Continuous approaches capture rotational equivariance but are often prohibitively computationally demanding. Discrete approaches offer more favorable computational performance but at the cost of equivariance. We develop a hybrid discrete-continuous (DISCO) group convolution that is simultaneously equivariant and computationally scalable to high-resolution. While our framework can be applied to any compact group, we specialize to the sphere. Our DISCO spherical convolutions exhibit $\text{SO}(3)$ rotational equivariance, where $\text{SO}(n)$ is the special orthogonal group representing rotations in $n$-dimensions. When restricting rotations of the convolution to the quotient space $\text{SO}(3)/\text{SO}(2)$ for further computational enhancements, we recover a form of asymptotic $\text{SO}(3)$ rotational equivariance. Through a sparse tensor implementation we achieve linear scaling in number of pixels on the sphere for both computational cost and memory usage. For 4k spherical images we realize a saving of $10^9$ in computational cost and $10^4$ in memory usage when compared to the most efficient alternative equivariant spherical convolution. We apply the DISCO spherical CNN framework to a number of benchmark dense-prediction problems on the sphere, such as semantic segmentation and depth estimation, on all of which we achieve the state-of-the-art performance., Comment: 19 pages, 7 figures, accepted by ICLR 2023

We resurrect the infamous harmonic mean estimator for computing the marginal likelihood (Bayesian evidence) and solve its problematic large variance. The marginal likelihood is a key component of Bayesian model selection to evaluate model posterior probabilities; however, its computation is challenging. The original harmonic mean estimator, first proposed by Newton and Raftery in 1994, involves computing the harmonic mean of the likelihood given samples from the posterior. It was immediately realised that the original estimator can fail catastrophically since its variance can become very large (possibly not finite). A number of variants of the harmonic mean estimator have been proposed to address this issue although none have proven fully satisfactory. We present the \emph{learnt harmonic mean estimator}, a variant of the original estimator that solves its large variance problem. This is achieved by interpreting the harmonic mean estimator as importance sampling and introducing a new target distribution. The new target distribution is learned to approximate the optimal but inaccessible target, while minimising the variance of the resulting estimator. Since the estimator requires samples of the posterior only, it is agnostic to the sampling strategy used. We validate the estimator on a variety of numerical experiments, including a number of pathological examples where the original harmonic mean estimator fails catastrophically. We also consider a cosmological application, where our approach leads to $\sim$ 3 to 6 times more samples than current state-of-the-art techniques in 1/3 of the time. In all cases our learnt harmonic mean estimator is shown to be highly accurate. The estimator is computationally scalable and can be applied to problems of dimension $O(10^3)$ and beyond. Code implementing the learnt harmonic mean estimator is made publicly available, Comment: 42 pages, 10 figures, code available at https://github.com/astro-informatics/harmonic

The Vera C. Rubin Observatory will increase the number of observed supernovae (SNe) by an order of magnitude; however, it is impossible to spectroscopically confirm the class for all the SNe discovered. Thus, photometric classification is crucial but its accuracy depends on the not-yet-finalized observing strategy of Rubin Observatory's Legacy Survey of Space and Time (LSST). We quantitatively analyze the impact of the LSST observing strategy on SNe classification using simulated multi-band light curves from the Photometric LSST Astronomical Time-Series Classification Challenge (PLAsTiCC). First, we augment the simulated training set to be representative of the photometric redshift distribution per supernovae class, the cadence of observations, and the flux uncertainty distribution of the test set. Then we build a classifier using the photometric transient classification library snmachine, based on wavelet features obtained from Gaussian process fits, yielding similar performance to the winning PLAsTiCC entry. We study the classification performance for SNe with different properties within a single simulated observing strategy. We find that season length is important, with light curves of 150 days yielding the highest performance. Cadence also has an important impact on SNe classification; events with median inter-night gap <3.5 days yield higher classification performance. Interestingly, we find that large gaps (>10 days) in light curve observations do not impact performance if sufficient observations are available on either side, due to the effectiveness of the Gaussian process interpolation. This analysis is the first exploration of the impact of observing strategy on photometric supernova classification with LSST., Comment: 18 pages, 13 figures. Changed to match version accepted by the Astrophysical Journal Supplement Series (accepted 28/10/2021). Software publicly available at https://github.com/LSSTDESC/snmachine

Bayesian model selection provides a powerful framework for objectively comparing models directly from observed data, without reference to ground truth data. However, Bayesian model selection requires the computation of the marginal likelihood (model evidence), which is computationally challenging, prohibiting its use in many high-dimensional Bayesian inverse problems. With Bayesian imaging applications in mind, in this work we present the proximal nested sampling methodology to objectively compare alternative Bayesian imaging models for applications that use images to inform decisions under uncertainty. The methodology is based on nested sampling, a Monte Carlo approach specialised for model comparison, and exploits proximal Markov chain Monte Carlo techniques to scale efficiently to large problems and to tackle models that are log-concave and not necessarily smooth (e.g., involving l_1 or total-variation priors). The proposed approach can be applied computationally to problems of dimension O(10^6) and beyond, making it suitable for high-dimensional inverse imaging problems. It is validated on large Gaussian models, for which the likelihood is available analytically, and subsequently illustrated on a range of imaging problems where it is used to analyse different choices of dictionary and measurement model.

This work presents the construction of a novel spherical wavelet basis designed for incomplete spherical datasets, i.e. datasets which are missing in a particular region of the sphere. The eigenfunctions of the Slepian spatial-spectral concentration problem (the Slepian functions) are a set of orthogonal basis functions which are more concentrated within a defined region. Slepian functions allow one to compute a convolution on the incomplete sphere by leveraging the recently proposed sifting convolution and extending it to any set of basis functions. Through a tiling of the Slepian harmonic line, one may construct scale-discretised wavelets. An illustration is presented based on an example region on the sphere defined by the topographic map of the Earth. The Slepian wavelets and corresponding wavelet coefficients are constructed from this region and are used in a straightforward denoising example., Comment: 12 pages, 14 figures

Future surveys such as the Legacy Survey of Space and Time (LSST) of the Vera C. Rubin Observatory will observe an order of magnitude more astrophysical transient events than any previous survey before. With this deluge of photometric data, it will be impossible for all such events to be classified by humans alone. Recent efforts have sought to leverage machine learning methods to tackle the challenge of astronomical transient classification, with ever improving success. Transformers are a recently developed deep learning architecture, first proposed for natural language processing, that have shown a great deal of recent success. In this work we develop a new transformer architecture, which uses multi-head self attention at its core, for general multi-variate time-series data. Furthermore, the proposed time-series transformer architecture supports the inclusion of an arbitrary number of additional features, while also offering interpretability. We apply the time-series transformer to the task of photometric classification, minimising the reliance of expert domain knowledge for feature selection, while achieving results comparable to state-of-the-art photometric classification methods. We achieve a logarithmic-loss of 0.507 on imbalanced data in a representative setting using data from the Photometric LSST Astronomical Time-Series Classification Challenge (PLAsTiCC). Moreover, we achieve a micro-averaged receiver operating characteristic area under curve of 0.98 and micro-averaged precision-recall area under curve of 0.87., Comment: Manuscript Accepted to RAS Techniques and Instruments. 15 pages, 12 figures

We develop variational regularization methods which leverage sparsity-promoting priors to solve severely ill posed inverse problems defined on the 3D ball (i.e. the solid sphere). Our method solves the problem natively on the ball and thus does not suffer from discontinuities that plague alternate approaches where each spherical shell is considered independently. Additionally, we leverage advances in probability density theory to produce Bayesian variational methods which benefit from the computational efficiency of advanced convex optimization algorithms, whilst supporting principled uncertainty quantification. We showcase these variational regularization and uncertainty quantification techniques on an illustrative example. The C++ code discussed throughout is provided under a GNU general public license.

Inverse problems defined naturally on the sphere are becoming increasingly of interest. In this article we provide a general framework for evaluation of inverse problems on the sphere, with a strong emphasis on flexibility and scalability. We consider flexibility with respect to the prior selection (regularization), the problem definition - specifically the problem formulation (constrained/unconstrained) and problem setting (analysis/synthesis) - and optimization adopted to solve the problem. We discuss and quantify the trade-offs between problem formulation and setting. Crucially, we consider the Bayesian interpretation of the unconstrained problem which, combined with recent developments in probability density theory, permits rapid, statistically principled uncertainty quantification (UQ) in the spherical setting. Linearity is exploited to significantly increase the computational efficiency of such UQ techniques, which in some cases are shown to permit analytic solutions. We showcase this reconstruction framework and UQ techniques on a variety of spherical inverse problems. The code discussed throughout is provided under a GNU general public license, in both C++ and Python.

Convolutional neural networks (CNNs) constructed natively on the sphere have been developed recently and shown to be highly effective for the analysis of spherical data. While an efficient framework has been formulated, spherical CNNs are nevertheless highly computationally demanding; typically they cannot scale beyond spherical signals of thousands of pixels. We develop scattering networks constructed natively on the sphere that provide a powerful representational space for spherical data. Spherical scattering networks are computationally scalable and exhibit rotational equivariance, while their representational space is invariant to isometries and provides efficient and stable signal representations. By integrating scattering networks as an additional type of layer in the generalized spherical CNN framework, we show how they can be leveraged to scale spherical CNNs to the high-resolution data typical of many practical applications, with spherical signals of many tens of megapixels and beyond., Comment: 18 pages, 6 figures, accepted by ICLR, code at https://www.kagenova.com/products/fourpiAI/

Many problems across computer vision and the natural sciences require the analysis of spherical data, for which representations may be learned efficiently by encoding equivariance to rotational symmetries. We present a generalized spherical CNN framework that encompasses various existing approaches and allows them to be leveraged alongside each other. The only existing non-linear spherical CNN layer that is strictly equivariant has complexity $\mathcal{O}(C^2L^5)$, where $C$ is a measure of representational capacity and $L$ the spherical harmonic bandlimit. Such a high computational cost often prohibits the use of strictly equivariant spherical CNNs. We develop two new strictly equivariant layers with reduced complexity $\mathcal{O}(CL^4)$ and $\mathcal{O}(CL^3 \log L)$, making larger, more expressive models computationally feasible. Moreover, we adopt efficient sampling theory to achieve further computational savings. We show that these developments allow the construction of more expressive hybrid models that achieve state-of-the-art accuracy and parameter efficiency on spherical benchmark problems., Comment: 20 pages, 4 figures, accepted by ICLR, code at https://www.kagenova.com/products/fourpiAI/

Since the 1970s, much of traditional interferometric imaging has been built around variations of the CLEAN algorithm, in both terminology, methodology, and algorithm development. Recent developments in applying new algorithms from convex optimization to interferometry has allowed old concepts to be viewed from a new perspective, ranging from image restoration to the development of computationally distributed algorithms. We present how this has ultimately led the authors to new perspectives in wide-field imaging, allowing for the first full individual non-coplanar corrections applied during imaging over extremely wide-fields of view for the Murchison Widefield Array (MWA) telescope. Furthermore, this same mathematical framework has provided a novel understanding of wide-band polarimetry at low frequencies, where instrumental channel depolarization can be corrected through the new $\delta\lambda^2$-projection algorithm. This is a demonstration that new algorithm development outside of traditional radio astronomy is valuable for the new theoretical and practical perspectives gained. These perspectives are timely with the next generation of radio telescopes coming online., Comment: 4 pages, 1 figure. URSI GASS 2020, Rome, Italy, 29 August - 5 September 2020

A novel spherical convolution is defined through the sifting property of the Dirac delta on the sphere. The so-called sifting convolution is defined by the inner product of one function with a translated version of another, but with the adoption of an alternative translation operator on the sphere. This translation operator follows by analogy with the Euclidean translation when viewed in harmonic space. The sifting convolution satisfies a variety of desirable properties that are lacking in alternate definitions, namely: it supports directional kernels; it has an output which remains on the sphere; and is efficient to compute. An illustration of the sifting convolution on a topographic map of the Earth demonstrates that it supports directional kernels to perform anisotropic filtering, while its output remains on the sphere., Comment: 5 pages, 3 figures

Radio astronomy is transitioning to a big-data era due to the emerging generation of radio interferometric (RI) telescopes, such as the Square Kilometre Array (SKA), which will acquire massive volumes of data. In this article we review methods proposed recently to resolve the ill-posed inverse problem of imaging the raw visibilities acquired by RI telescopes in the big-data scenario. We focus on the recently proposed online reconstruction method [4] and the considerable savings in data storage requirements and computational cost that it yields., Comment: 4 pages; 2 figures; URSI GASS 2020. arXiv admin note: substantial text overlap with arXiv:1712.04462

The deconvolution, or cleaning, of radio interferometric images often involves computing model visibilities from a list of clean components, in order that the contribution from the model can be subtracted from the observed visibilities. This step is normally performed using a forward fast Fourier transform (FFT), followed by a 'degridding' step that interpolates over the uv plane to construct the model visibilities. An alternative approach is to calculate the model visibilities directly by summing over all the members of the clean component list, which is a more accurate method that can also be much slower. However, if the clean components are used to construct a model image on the surface of the celestial sphere then the model visibilities can be generated directly from the wavelet coefficients, and the sparsity of the model means that most of these coefficients are zero, and can be ignored. We have constructed a prototype imager that uses a spherical-wavelet representation of the model image to generate model visibilities during each major cycle, and find empirically that the execution time scales with the wavelet resolution level, J, as O(1.07 J), and with the number of distinct clean components, N_C, as O(N_C). The prototype organises the wavelet coefficients into a tree structure, and does not store or process the zero wavelet coefficients., Comment: 12 pages, 12 figures

We present a new algorithm to perform wide-field radio interferometric image reconstruction, with exact non-coplanar correction, that scales to big-data. This algorithm allows us to image 2 billion visibilities on 50 nodes of a computing cluster for a 25 by 25 degree field of view, in a little over an hour. We build on the recently developed distributed $w$-stacking $w$-projection hybrid algorithm, extending it to include a new distributed degridding algorithm that balances the computational load of the $w$-projection gridding kernels. The implementation of our algorithm is made publicly available in the PURIFY software package. Wide-field image reconstruction for data sets of this size cannot be performed effectively using the allocated computational resources without computational load balancing, demonstrating that our algorithms are critical for next-generation wide-field radio interferometers., Comment: 4 pages, 1 figure

The concept of a recently proposed small-scale interferometric optical imaging device, an instrument known as the Segmented Planar Imaging Detector for Electro-optical Reconnaissance (SPIDER), is of great interest for its possible applications in astronomy and space science. Due to low weight, low power consumption, and high resolution, the SPIDER telescope could replace the large space telescopes that exist today. Unlike traditional optical interferometry the SPIDER accurately retrieves both phase and amplitude information, making the measurement process analogous to a radio interferometer. State of the art sparse radio interferometric image reconstruction techniques have been gaining traction in radio astronomy and reconstruct accurate images of the radio sky. In this work we describe algorithms from radio interferometric imaging and sparse image reconstruction and demonstrate their application to the SPIDER concept telescope through simulated observation and reconstruction of the optical sky. Such algorithms are important for providing high fidelity images from SPIDER observations, helping to power the SPIDER concept for scientific and astronomical analysis., Comment: 4 Pages, 2 Figures, 1 Table

Next generation radio interferometric telescopes are entering an era of big data with extremely large data sets. While these telescopes can observe the sky in higher sensitivity and resolution than before, computational challenges in image reconstruction need to be overcome to realize the potential of forthcoming telescopes. New methods in sparse image reconstruction and convex optimization techniques (cf. compressive sensing) have shown to produce higher fidelity reconstructions of simulations and real observations than traditional methods. This article presents distributed and parallel algorithms and implementations to perform sparse image reconstruction, with significant practical considerations that are important for implementing these algorithms for Big Data. We benchmark the algorithms presented, showing that they are considerably faster than their serial equivalents. We then pre-sample gridding kernels to scale the distributed algorithms to larger data sizes, showing application times for 1 Gb to 2.4 Tb data sets over 25 to 100 nodes for up to 50 billion visibilities, and find that the run-times for the distributed algorithms range from 100 milliseconds to 3 minutes per iteration. This work presents an important step in working towards computationally scalable and efficient algorithms and implementations that are needed to image observations of both extended and compact sources from next generation radio interferometers such as the SKA. The algorithms are implemented in the latest versions of the SOPT (https://github.com/astro-informatics/sopt) and PURIFY (https://github.com/astro-informatics/purify) software packages {(Versions 3.1.0)}, which have been released alongside of this article., Comment: 25 pages, 5 figures

Polarization of radiation is a powerful tool to study cosmic magnetism and analysis of polarization can be used as a diagnostic tool for large-scale structures. In this paper, we present a solid theoretical foundation for using polarized light to investigate large-scale magnetic field structures: the cosmological polarized radiative transfer (CPRT) formulation. The CPRT formulation is fully covariant. It accounts for cosmological and relativistic effects in a self-consistent manner and explicitly treats Faraday rotation, as well as Faraday conversion, emission, and absorption processes. The formulation is derived from the first principles of conservation of phase-space volume and photon number. Without loss of generality, we consider a flat Friedmann-Robertson-Walker (FRW) space-time metric and construct the corresponding polarized radiative transfer equations. We propose an all-sky CPRT calculation algorithm, based on a ray-tracing method, which allows cosmological simulation results to be incorporated and, thereby, model templates of polarization maps to be constructed. Such maps will be crucial in our interpretation of polarized data, such as those to be collected by the Square Kilometer Array (SKA). We describe several tests which are used for verifying the code and demonstrate applications in the study of the polarization signatures in different distributions of electron number density and magnetic fields. We present a pencil-beam CPRT calculation and an all-sky calculation, using a simulated galaxy cluster or a model magnetized universe obtained from GCMHD+ simulations as the respective input structures. The implications on large-scale magnetic field studies are discussed; remarks on the standard methods using rotation measure are highlighted., Comment: 32 pages, 14 figures

We review the measurements of dark energy enabled by observations of the Deep Drilling Fields and the optimization of survey design for cosmological measurements. This white paper is the result of efforts by the LSST DESC Observing Strategy Task Force (OSTF), which represents the entire collaboration, and aims to make recommendations on observing strategy for the DDFs that will benefit all cosmological analyses with LSST. It is accompanied by the DESC-WFD white paper (Lochner et al.). We argue for altering the nominal deep drilling plan to have $>6$ month seasons, interweaving $gri$ and $zy$ observations every 3 days with 2, 4, 8, 25, 4 visits in $grizy$, respectively. These recommendations are guided by metrics optimizing constraints on dark energy and mitigation of systematic uncertainties, including specific requirements on total number of visits after Y1 and Y10 for photometric redshifts (photo-$z$) and weak lensing systematics. We specify the precise locations for the previously-chosen LSST deep fields (ELAIS-S1, XMM-LSS, CDF-S, and COSMOS) and recommend Akari Deep Field South as the planned fifth deep field in order to synergize with Euclid and WFIRST. Our recommended DDF strategy uses $6.2\%$ of the LSST survey time. We briefly discuss synergy with white papers from other collaborations, as well as additional mini-surveys and Target-of-Opportunity programs that lead to better measurements of dark energy., Comment: The LSST DESC response (DDF) to the Call for White Papers on LSST Cadence Optimization. Comments welcome

Cosmology is one of the four science pillars of LSST, which promises to be transformative for our understanding of dark energy and dark matter. The LSST Dark Energy Science Collaboration (DESC) has been tasked with deriving constraints on cosmological parameters from LSST data. Each of the cosmological probes for LSST is heavily impacted by the choice of observing strategy. This white paper is written by the LSST DESC Observing Strategy Task Force (OSTF), which represents the entire collaboration, and aims to make recommendations on observing strategy that will benefit all cosmological analyses with LSST. It is accompanied by the DESC DDF (Deep Drilling Fields) white paper (Scolnic et al.). We use a variety of metrics to understand the effects of the observing strategy on measurements of weak lensing, large-scale structure, clusters, photometric redshifts, supernovae, strong lensing and kilonovae. In order to reduce systematic uncertainties, we conclude that the current baseline observing strategy needs to be significantly modified to result in the best possible cosmological constraints. We provide some key recommendations: moving the WFD (Wide-Fast-Deep) footprint to avoid regions of high extinction, taking visit pairs in different filters, changing the 2x15s snaps to a single exposure to improve efficiency, focusing on strategies that reduce long gaps (>15 days) between observations, and prioritizing spatial uniformity at several intervals during the 10-year survey., Comment: The LSST DESC response (WFD) to the Call for White Papers on LSST Cadence Optimization. Comments welcome

Inverse problems play a key role in modern image/signal processing methods. However, since they are generally ill-conditioned or ill-posed due to lack of observations, their solutions may have significant intrinsic uncertainty. Analysing and quantifying this uncertainty is very challenging, particularly in high-dimensional problems and problems with non-smooth objective functionals (e.g. sparsity-promoting priors). In this article, a series of strategies to visualise this uncertainty are presented, e.g. highest posterior density credible regions, and local credible intervals (cf. error bars) for individual pixels and superpixels. Our methods support non-smooth priors for inverse problems and can be scaled to high-dimensional settings. Moreover, we present strategies to automatically set regularisation parameters so that the proposed uncertainty quantification (UQ) strategies become much easier to use. Also, different kinds of dictionaries (complete and over-complete) are used to represent the image/signal and their performance in the proposed UQ methodology is investigated., Comment: 5 pages, 5 figures

We explore a new Bayesian method of detecting galaxies from radio interferometric data of the faint sky. Working in the Fourier domain, we fit a single, parameterised galaxy model to simulated visibility data of star-forming galaxies. The resulting multimodal posterior distribution is then sampled using a multimodal nested sampling algorithm such as MultiNest. For each galaxy, we construct parameter estimates for the position, flux, scale-length and ellipticities from the posterior samples. We first test our approach on simulated SKA1-MID visibility data of up to 100 galaxies in the field of view, considering a typical weak lensing survey regime (SNR $\ge 10$) where 98% of the input galaxies are detected with no spurious source detections. We then explore the low SNR regime, finding our approach reliable in galaxy detection and providing in particular high accuracy in positional estimates down to SNR $\sim 5$. The presented method does not require transformation of visibilities to the image domain, and requires no prior knowledge of the number of galaxies in the field of view, thus could become a useful tool for constructing accurate radio galaxy catalogs in the future., Comment: 11 pages, 11 figures. Accepted for publication in MNRAS

The Photometric LSST Astronomical Time Series Classification Challenge (PLAsTiCC) is an open data challenge to classify simulated astronomical time-series data in preparation for observations from the Large Synoptic Survey Telescope (LSST), which will achieve first light in 2019 and commence its 10-year main survey in 2022. LSST will revolutionize our understanding of the changing sky, discovering and measuring millions of time-varying objects. In this challenge, we pose the question: how well can we classify objects in the sky that vary in brightness from simulated LSST time-series data, with all its challenges of non-representativity? In this note we explain the need for a data challenge to help classify such astronomical sources and describe the PLAsTiCC data set and Kaggle data challenge, noting that while the references are provided for context, they are not needed to participate in the challenge., Comment: Research note to accompany the https://www.kaggle.com/c/PLAsTiCC-2018 challenge

For the representation of spin-$s$ band-limited functions on the sphere, we propose a sampling scheme with optimal number of samples equal to the number of degrees of freedom of the function in harmonic space. In comparison to the existing sampling designs, which require ${\sim}2L^2$ samples for the representation of spin-$s$ functions band-limited at $L$, the proposed scheme requires $N_o=L^2-s^2$ samples for the accurate computation of the spin-$s$ spherical harmonic transform~($s$-SHT). For the proposed sampling scheme, we also develop a method to compute the $s$-SHT. We place the samples in our design scheme such that the matrices involved in the computation of $s$-SHT are well-conditioned. We also present a multi-pass $s$-SHT to improve the accuracy of the transform. We also show the proposed sampling design exhibits superior geometrical properties compared to existing equiangular and Gauss-Legendre sampling schemes, and enables accurate computation of the $s$-SHT corroborated through numerical experiments., Comment: 5 pages, 2 figures

The standard wide-field imaging technique, the $w$-projection, allows correction for wide-fields of view for non-coplanar radio interferometric arrays. However, calculating exact corrections for each measurement has not been possible due to the amount of computation required at high resolution and with the large number of visibilities from current interferometers. The required accuracy and computational cost of these corrections is one of the largest unsolved challenges facing next generation radio interferometers such as the Square Kilometre Array. We show that the same calculation can be performed with a radially symmetric $w$-projection kernel, where we use one dimensional adaptive quadrature to calculate the resulting Hankel transform, decreasing the computation required for kernel generation by several orders of magnitude, whilst preserving the accuracy. We confirm that the radial $w$-projection kernel is accurate to approximately 1% by imaging the zero-spacing with an added $w$-term. We demonstrate the potential of our radially symmetric $w$-projection kernel via sparse image reconstruction, using the software package PURIFY. We develop a distributed $w$-stacking and $w$-projection hybrid algorithm. We apply this algorithm to individually correct for non-coplanar effects in 17.5 million visibilities over a $25$ by $25$ degree field of view MWA observation for image reconstruction. Such a level of accuracy and scalability is not possible with standard $w$-projection kernel generation methods. This demonstrates that we can scale to a large number of measurements with large image sizes whilst still maintaining both speed and accuracy., Comment: 9 Figures, 19 Pages. Accepted to ApJ

The emerging generation of radio interferometric (RI) telescopes, such as the Square Kilometre Array (SKA), will acquire massive volumes of data and transition radio astronomy to a big-data era. The ill-posed inverse problem of imaging the raw visibilities acquired by RI telescopes will become significantly more computationally challenging, particularly in terms of data storage and computational cost. Current RI imaging methods, such as CLEAN, its variants, and compressive sensing approaches (sparse regularisation), have yielded excellent reconstruction fidelity. However, scaling these methods to big-data remains difficult if not impossible in some cases. All state-of-the-art methods in RI imaging lack the ability to process data streams as they are acquired during the data observation stage. Such approaches are referred to as online processing methods. We present an online sparse regularisation methodology for RI imaging. Image reconstruction is performed simultaneously with data acquisition, where observed visibilities are assimilated into the reconstructed image as they arrive and then discarded. Since visibilities are processed online, good reconstructions are recovered much faster than standard (offline) methods which cannot start until the data acquisition stage completes. Moreover, the online method provides additional computational savings and, most importantly, dramatically reduces data storage requirements. Theoretically, the reconstructed images are of the same fidelity as those recovered by the equivalent offline approach and, in practice, very similar reconstruction fidelity is achieved. We anticipate online imaging techniques, as proposed here, will be critical in scaling RI imaging to the emerging big-data era of radio astronomy., Comment: 14 pages, 5 figures

Uncertainty quantification is a critical missing component in radio interferometric imaging that will only become increasingly important as the big-data era of radio interferometry emerges. Since radio interferometric imaging requires solving a high-dimensional, ill-posed inverse problem, uncertainty quantification is difficult but also critical to the accurate scientific interpretation of radio observations. Statistical sampling approaches to perform Bayesian inference, like Markov Chain Monte Carlo (MCMC) sampling, can in principle recover the full posterior distribution of the image, from which uncertainties can then be quantified. However, traditional high-dimensional sampling methods are generally limited to smooth (e.g. Gaussian) priors and cannot be used with sparsity-promoting priors. Sparse priors, motivated by the theory of compressive sensing, have been shown to be highly effective for radio interferometric imaging. In this article proximal MCMC methods are developed for radio interferometric imaging, leveraging proximal calculus to support non-differential priors, such as sparse priors, in a Bayesian framework. Furthermore, three strategies to quantify uncertainties using the recovered posterior distribution are developed: (i) local (pixel-wise) credible intervals to provide error bars for each individual pixel; (ii) highest posterior density credible regions; and (iii) hypothesis testing of image structure. These forms of uncertainty quantification provide rich information for analysing radio interferometric observations in a statistically robust manner., Comment: 16 pages, 7 figures, see companion article in this arXiv listing

Uncertainty quantification is a critical missing component in radio interferometric imaging that will only become increasingly important as the big-data era of radio interferometry emerges. Statistical sampling approaches to perform Bayesian inference, like Markov Chain Monte Carlo (MCMC) sampling, can in principle recover the full posterior distribution of the image, from which uncertainties can then be quantified. However, for massive data sizes, like those anticipated from the Square Kilometre Array (SKA), it will be difficult if not impossible to apply any MCMC technique due to its inherent computational cost. We formulate Bayesian inference problems with sparsity-promoting priors (motivated by compressive sensing), for which we recover maximum a posteriori (MAP) point estimators of radio interferometric images by convex optimisation. Exploiting recent developments in the theory of probability concentration, we quantify uncertainties by post-processing the recovered MAP estimate. Three strategies to quantify uncertainties are developed: (i) highest posterior density credible regions; (ii) local credible intervals (cf. error bars) for individual pixels and superpixels; and (iii) hypothesis testing of image structure. These forms of uncertainty quantification provide rich information for analysing radio interferometric observations in a statistically robust manner. Our MAP-based methods are approximately $10^5$ times faster computationally than state-of-the-art MCMC methods and, in addition, support highly distributed and parallelised algorithmic structures. For the first time, our MAP-based techniques provide a means of quantifying uncertainties for radio interferometric imaging for realistic data volumes and practical use, and scale to the emerging big-data era of radio astronomy., Comment: 13 pages, 10 figures, see companion article in this arXiv listing

The Large Synoptic Survey Telescope is designed to provide an unprecedented optical imaging dataset that will support investigations of our Solar System, Galaxy and Universe, across half the sky and over ten years of repeated observation. However, exactly how the LSST observations will be taken (the observing strategy or "cadence") is not yet finalized. In this dynamically-evolving community white paper, we explore how the detailed performance of the anticipated science investigations is expected to depend on small changes to the LSST observing strategy. Using realistic simulations of the LSST schedule and observation properties, we design and compute diagnostic metrics and Figures of Merit that provide quantitative evaluations of different observing strategies, analyzing their impact on a wide range of proposed science projects. This is work in progress: we are using this white paper to communicate to each other the relative merits of the observing strategy choices that could be made, in an effort to maximize the scientific value of the survey. The investigation of some science cases leads to suggestions for new strategies that could be simulated and potentially adopted. Notably, we find motivation for exploring departures from a spatially uniform annual tiling of the sky: focusing instead on different parts of the survey area in different years in a "rolling cadence" is likely to have significant benefits for a number of time domain and moving object astronomy projects. The communal assembly of a suite of quantified and homogeneously coded metrics is the vital first step towards an automated, systematic, science-based assessment of any given cadence simulation, that will enable the scheduling of the LSST to be as well-informed as possible., Comment: 312 pages, 90 figures. Browse the current version at https://github.com/LSSTScienceCollaborations/ObservingStrategy, new contributions welcome!

Convergence maps of the integrated matter distribution are a key science result from weak gravitational lensing surveys. To date, recovering convergence maps has been performed using a planar approximation of the celestial sphere. However, with the increasing area of sky covered by dark energy experiments, such as Euclid, the Large Synoptic Survey Telescope (LSST), and the Wide Field Infrared Survey Telescope (WFIRST), this assumption will no longer be valid. We recover convergence fields on the celestial sphere using an extension of the Kaiser-Squires estimator to the spherical setting. Through simulations we study the error introduced by planar approximations. Moreover, we examine how best to recover convergence maps in the planar setting, considering a variety of different projections and defining the local rotations that are required when projecting spin fields such as cosmic shear. For the sky coverages typical of future surveys, errors introduced by projection effects can be of order tens of percent, exceeding 50% in some cases. The stereographic projection, which is conformal and so preserves local angles, is the most effective planar projection. In any case, these errors can be avoided entirely by recovering convergence fields directly on the celestial sphere. We apply the spherical Kaiser-Squires mass-mapping method presented to the public Dark Energy Survey (DES) science verification data to recover convergence maps directly on the celestial sphere., Comment: 18 Pages, 10 Figures, comments welcome, code will made publicly available until then is available on request