Showing 1–50 of 64 results for author: Pereyra, M

Bayesian computation with generative diffusion models by Multilevel Monte Carlo

Authors: Abdul-Lateef Haji-Ali, Marcelo Pereyra, Luke Shaw, Konstantinos Zygalakis

Abstract: Generative diffusion models have recently emerged as a powerful strategy to perform stochastic sampling in Bayesian inverse problems, delivering remarkably accurate solutions for a wide range of challenging applications. However, diffusion models often require a large number of neural function evaluations per sample in order to deliver accurate posterior samples. As a result, using diffusion model… ▽ More Generative diffusion models have recently emerged as a powerful strategy to perform stochastic sampling in Bayesian inverse problems, delivering remarkably accurate solutions for a wide range of challenging applications. However, diffusion models often require a large number of neural function evaluations per sample in order to deliver accurate posterior samples. As a result, using diffusion models as stochastic samplers for Monte Carlo integration in Bayesian computation can be highly computationally expensive. This cost is especially high in large-scale inverse problems such as computational imaging, which rely on large neural networks that are expensive to evaluate. With Bayesian imaging problems in mind, this paper presents a Multilevel Monte Carlo strategy that significantly reduces the cost of Bayesian computation with diffusion models. This is achieved by exploiting cost-accuracy trade-offs inherent to diffusion models to carefully couple models of different levels of accuracy in a manner that significantly reduces the overall cost of the calculation, without reducing the final accuracy. The effectiveness of the proposed Multilevel Monte Carlo approach is demonstrated with three canonical computational imaging problems, where we observe a $4\times$-to-$8\times$ reduction in computational cost compared to conventional Monte Carlo averaging. △ Less

Submitted 23 September, 2024; originally announced September 2024.

Comments: 13 images

Empirical Bayesian image restoration by Langevin sampling with a denoising diffusion implicit prior

Authors: Charlesquin Kemajou Mbakam, Jean-Francois Giovannelli, Marcelo Pereyra

Abstract: Score-based diffusion methods provide a powerful strategy to solve image restoration tasks by flexibly combining a pre-trained foundational prior model with a likelihood function specified during test time. Such methods are predominantly derived from two stochastic processes: reversing Ornstein-Uhlenbeck, which underpins the celebrated denoising diffusion probabilistic models (DDPM) and denoising… ▽ More Score-based diffusion methods provide a powerful strategy to solve image restoration tasks by flexibly combining a pre-trained foundational prior model with a likelihood function specified during test time. Such methods are predominantly derived from two stochastic processes: reversing Ornstein-Uhlenbeck, which underpins the celebrated denoising diffusion probabilistic models (DDPM) and denoising diffusion implicit models (DDIM), and the Langevin diffusion process. The solutions delivered by DDPM and DDIM are often remarkably realistic, but they are not always consistent with measurements because of likelihood intractability issues and the associated required approximations. Alternatively, using a Langevin process circumvents the intractable likelihood issue, but usually leads to restoration results of inferior quality and longer computing times. This paper presents a novel and highly computationally efficient image restoration method that carefully embeds a foundational DDPM denoiser within an empirical Bayesian Langevin algorithm, which jointly calibrates key model hyper-parameters as it estimates the model's posterior mean. Extensive experimental results on three canonical tasks (image deblurring, super-resolution, and inpainting) demonstrate that the proposed approach improves on state-of-the-art strategies both in image estimation accuracy and computing time. △ Less

Submitted 6 September, 2024; originally announced September 2024.

Comments: 24 pages

Do Bayesian imaging methods report trustworthy probabilities?

Authors: David Y. W. Thong, Charlesquin Kemajou Mbakam, Marcelo Pereyra

Abstract: Bayesian statistics is a cornerstone of imaging sciences, underpinning many and varied approaches from Markov random fields to score-based denoising diffusion models. In addition to powerful image estimation methods, the Bayesian paradigm also provides a framework for uncertainty quantification and for using image data as quantitative evidence. These probabilistic capabilities are important for th… ▽ More Bayesian statistics is a cornerstone of imaging sciences, underpinning many and varied approaches from Markov random fields to score-based denoising diffusion models. In addition to powerful image estimation methods, the Bayesian paradigm also provides a framework for uncertainty quantification and for using image data as quantitative evidence. These probabilistic capabilities are important for the rigorous interpretation of experimental results and for robust interfacing of quantitative imaging pipelines with scientific and decision-making processes. However, are the probabilities delivered by existing Bayesian imaging methods meaningful under replication of an experiment, or are they only meaningful as subjective measures of belief? This paper presents a Monte Carlo method to explore this question. We then leverage the proposed Monte Carlo method and run a large experiment requiring 1,000 GPU-hours to probe the accuracy of five canonical Bayesian imaging methods that are representative of some of the main Bayesian imaging strategies from the past decades (a score-based denoising diffusion technique, a plug-and-play Langevin algorithm utilising a Lipschitz-regularised DnCNN denoiser, a Bayesian method with a dictionary-based prior trained subject to a log-concavity constraint, an empirical Bayesian method with a total-variation prior, and a hierarchical Bayesian Gibbs sampler based on a Gaussian Markov random field model). We find that, a few cases, the probabilities reported by modern Bayesian imaging techniques are in broad agreement with long-term averages as observed over a large number of replication of an experiment, but existing Bayesian imaging methods are generally not able to deliver reliable uncertainty quantification results. △ Less

Submitted 13 May, 2024; originally announced May 2024.

MSC Class: 65J22 (Primary); 62F15 (Secondary); 68U10

An evaluation of the BALROG and RoboBA algorithms for determining the position of Fermi/GBM GRBs

Authors: K. Océlotl. C. López, Alan M. Watson, William H. Lee, Rosa L. Becerra, Margarita Pereyra

Abstract: The Fermi/GBM instrument is a vital source of detections of gamma-ray bursts and has an increasingly important role to play in understanding gravitational-wave transients. In both cases, its impact is increased by accurate positions with reliable uncertainties. We evaluate the RoboBA and BALROG algorithms for determining the position of gamma-ray bursts detected by the Fermi/GBM instrument. We con… ▽ More The Fermi/GBM instrument is a vital source of detections of gamma-ray bursts and has an increasingly important role to play in understanding gravitational-wave transients. In both cases, its impact is increased by accurate positions with reliable uncertainties. We evaluate the RoboBA and BALROG algorithms for determining the position of gamma-ray bursts detected by the Fermi/GBM instrument. We construct a sample of 54 bursts with detections both by Swift/BAT and by Fermi/GBM. We then compare the positions predicted by RoboBA and BALROG with the positions measured by BAT, which we can assume to be the true position. We find that RoboBA and BALROG are similarly precise for bright bursts whose uncertainties are dominated by systematic errors, but RoboBA performs better for faint bursts whose uncertainties are dominated by statistical noise. We further find that the uncertainties in the positions predicted by RoboBA are consistent with the distribution of position errors, whereas BALROG seems to be underestimating the uncertainties by a factor of about two. Additionally, we consider the implications of these results for the follow-up of the optical afterglows of Fermi/GBM bursts. In particular, for the DDOTI wide-field imager we conclude that a single pointing is best. Our sample would allow a similar study to be carried out for other telescopes. △ Less

Submitted 30 April, 2024; originally announced April 2024.

Comments: Revision sent to MNRAS after addressing minimal referee comments

Unsupervised Training of Convex Regularizers using Maximum Likelihood Estimation

Authors: Hong Ye Tan, Ziruo Cai, Marcelo Pereyra, Subhadip Mukherjee, Junqi Tang, Carola-Bibiane Schönlieb

Abstract: Imaging is a standard example of an inverse problem, where the task of reconstructing a ground truth from a noisy measurement is ill-posed. Recent state-of-the-art approaches for imaging use deep learning, spearheaded by unrolled and end-to-end models and trained on various image datasets. However, many such methods require the availability of ground truth data, which may be unavailable or expensi… ▽ More Imaging is a standard example of an inverse problem, where the task of reconstructing a ground truth from a noisy measurement is ill-posed. Recent state-of-the-art approaches for imaging use deep learning, spearheaded by unrolled and end-to-end models and trained on various image datasets. However, many such methods require the availability of ground truth data, which may be unavailable or expensive, leading to a fundamental barrier that can not be bypassed by choice of architecture. Unsupervised learning presents an alternative paradigm that bypasses this requirement, as they can be learned directly on noisy data and do not require any ground truths. A principled Bayesian approach to unsupervised learning is to maximize the marginal likelihood with respect to the given noisy measurements, which is intrinsically linked to classical variational regularization. We propose an unsupervised approach using maximum marginal likelihood estimation to train a convex neural network-based image regularization term directly on noisy measurements, improving upon previous work in both model expressiveness and dataset size. Experiments demonstrate that the proposed method produces priors that are near competitive when compared to the analogous supervised training method for various image corruption operators, maintaining significantly better generalization properties when compared to end-to-end methods. Moreover, we provide a detailed theoretical analysis of the convergence properties of our proposed algorithm. △ Less

Submitted 29 July, 2024; v1 submitted 8 April, 2024; originally announced April 2024.

MSC Class: 62C12; 62F15; 65C40; 65J22

A stochastic optimisation unadjusted Langevin method for empirical Bayesian estimation in semi-blind image deblurring problems

Authors: Charlesquin Kemajou Mbakam, Marcelo Pereyra, Jean-François Giovannelli

Abstract: This paper presents a novel stochastic optimisation methodology to perform empirical Bayesian inference in semi-blind image deconvolution problems. Given a blurred image and a parametric class of possible operators, the proposed optimisation approach automatically calibrates the parameters of the blur model by maximum marginal likelihood estimation, followed by (non-blind) image deconvolution by m… ▽ More This paper presents a novel stochastic optimisation methodology to perform empirical Bayesian inference in semi-blind image deconvolution problems. Given a blurred image and a parametric class of possible operators, the proposed optimisation approach automatically calibrates the parameters of the blur model by maximum marginal likelihood estimation, followed by (non-blind) image deconvolution by maximum-a-posteriori estimation conditionally to the estimated model parameters. In addition to the blur model, the proposed approach also automatically calibrates the noise variance as well as any regularisation parameters. The marginal likelihood of the blur, noise variance, and regularisation parameters is generally computationally intractable, as it requires calculating several integrals over the entire solution space. Our approach addresses this difficulty by using a stochastic approximation proximal gradient optimisation scheme, which iteratively solves such integrals by using a Moreau-Yosida regularised unadjusted Langevin Markov chain Monte Carlo algorithm. This optimisation strategy can be easily and efficiently applied to any model that is log-concave, and by using the same gradient and proximal operators that are required to compute the maximum-a-posteriori solution by convex optimisation. We provide convergence guarantees for the proposed optimisation scheme under realistic and easily verifiable conditions and subsequently demonstrate the effectiveness of the approach with a series of deconvolution experiments and comparisons with alternative strategies from the state of the art. △ Less

Submitted 7 March, 2024; originally announced March 2024.

Comments: 43 pages

Statistical modelling and Bayesian inversion for a Compton imaging system: application to radioactive source localisation

Authors: Cecilia Tarpau, Ming Fang, Konstantinos C. Zygalakis, Marcelo Pereyra, Angela Di Fulvio, Yoann Altmann

Abstract: This paper presents a statistical forward model for a Compton imaging system, called Compton imager. This system, under development at the University of Illinois Urbana Champaign, is a variant of Compton cameras with a single type of sensors which can simultaneously act as scatterers and absorbers. This imager is convenient for imaging situations requiring a wide field of view. The proposed statis… ▽ More This paper presents a statistical forward model for a Compton imaging system, called Compton imager. This system, under development at the University of Illinois Urbana Champaign, is a variant of Compton cameras with a single type of sensors which can simultaneously act as scatterers and absorbers. This imager is convenient for imaging situations requiring a wide field of view. The proposed statistical forward model is then used to solve the inverse problem of estimating the location and energy of point-like sources from observed data. This inverse problem is formulated and solved in a Bayesian framework by using a Metropolis within Gibbs algorithm for the estimation of the location, and an expectation-maximization algorithm for the estimation of the energy. This approach leads to more accurate estimation when compared with the deterministic standard back-projection approach, with the additional benefit of uncertainty quantification in the low photon imaging setting. △ Less

Submitted 16 February, 2024; v1 submitted 12 February, 2024; originally announced February 2024.

Scalable Bayesian uncertainty quantification with data-driven priors for radio interferometric imaging

Authors: Tobías I. Liaudat, Matthijs Mars, Matthew A. Price, Marcelo Pereyra, Marta M. Betcke, Jason D. McEwen

Abstract: 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 size… ▽ More 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. △ Less

Submitted 31 July, 2024; v1 submitted 30 November, 2023; originally announced December 2023.

Comments: 30 pages, 14 figures, 10 tables, code available at https://github.com/astro-informatics/QuantifAI

Equivariant Bootstrapping for Uncertainty Quantification in Imaging Inverse Problems

Authors: Julian Tachella, Marcelo Pereyra

Abstract: Scientific imaging problems are often severely ill-posed, and hence have significant intrinsic uncertainty. Accurately quantifying the uncertainty in the solutions to such problems is therefore critical for the rigorous interpretation of experimental results as well as for reliably using the reconstructed images as scientific evidence. Unfortunately, existing imaging methods are unable to quantify… ▽ More Scientific imaging problems are often severely ill-posed, and hence have significant intrinsic uncertainty. Accurately quantifying the uncertainty in the solutions to such problems is therefore critical for the rigorous interpretation of experimental results as well as for reliably using the reconstructed images as scientific evidence. Unfortunately, existing imaging methods are unable to quantify the uncertainty in the reconstructed images in a manner that is robust to experiment replications. This paper presents a new uncertainty quantification methodology based on an equivariant formulation of the parametric bootstrap algorithm that leverages symmetries and invariance properties commonly encountered in imaging problems. Additionally, the proposed methodology is general and can be easily applied with any image reconstruction technique, including unsupervised training strategies that can be trained from observed data alone, thus enabling uncertainty quantification in situations where there is no ground truth data available. We demonstrate the proposed approach with a series of numerical experiments and through comparisons with alternative uncertainty quantification strategies from the state-of-the-art, such as Bayesian strategies involving score-based diffusion models and Langevin samplers. In all our experiments, the proposed method delivers remarkably accurate high-dimensional confidence regions and outperforms the competing approaches in terms of estimation accuracy, uncertainty quantification accuracy, and computing time. △ Less

Submitted 20 October, 2023; v1 submitted 18 October, 2023; originally announced October 2023.

MSC Class: 68T07 ACM Class: G.3

Journal ref: AISTATS 2024 (Oral presentation)

Machine-Learning Enhanced Photometric Analysis of the Extremely Bright GRB 210822A

Authors: Camila Angulo-Valdez, Rosa L. Becerra, Margarita Pereyra, Keneth Garcia-Cifuentes, Felipe Vargas, Alan M. Watson, Fabio De Colle, Nissim Fraija, Nathaniel R. Butler, Maria G. Dainotti, Simone Dichiara, William H. Lee, Eleonora Troja, Joshua S. Bloom, J. Jesús González, Alexander S. Kutyrev, J. Xavier Prochaska, Enrico Ramirez-Ruiz, Michael G. Richer

Abstract: We present analytical and numerical models of the bright long GRB 210822A at $z=1.736$. The intrinsic extreme brightness exhibited in the optical, which is very similar to other bright GRBs (e.g., GRBs 080319B, 130427A, 160625A 190114C, and 221009A), makes GRB 210822A an ideal case for studying the evolution of this particular kind of GRB. We use optical data from the RATIR instrument starting at… ▽ More We present analytical and numerical models of the bright long GRB 210822A at $z=1.736$. The intrinsic extreme brightness exhibited in the optical, which is very similar to other bright GRBs (e.g., GRBs 080319B, 130427A, 160625A 190114C, and 221009A), makes GRB 210822A an ideal case for studying the evolution of this particular kind of GRB. We use optical data from the RATIR instrument starting at $T+315.9$ s, with publicly available optical data from other ground-based observatories, as well as Swift/UVOT, and X-ray data from the Swift/XRT instrument. The temporal profiles and spectral properties during the late stages align consistently with the conventional forward shock model, complemented by a reverse shock element that dominates optical emissions during the initial phases ($T<300$ s). Furthermore, we observe a break at $T=80000$s that we interpreted as evidence of a jet break, which constrains the opening angle to be about $θ_\mathrm{j}=(3-5)$ degrees. Finally, we apply a machine-learning technique to model the multi-wavelength light curve of GRB 210822A using the AFTERGLOWPY library. We estimate the angle of sight $θ_{obs}=(6.4 \pm 0.1) \times 10^{-1}$ degrees, the energy $E_0=(7.9 \pm 1.6)\times 10^{53}$ ergs, the electron index $p=2.54 \pm 0.10$, the thermal energy fraction in electrons $ε_\mathrm{e}=(4.63 \pm 0.91) \times 10^{-5}$ and in the magnetic field $ε_\mathrm{B}= (8.66 \pm 1.01) \times 10^{-6}$, the efficiency $χ= 0.89 \pm 0.01$, and the density of the surrounding medium $n_\mathrm{0} = 0.85 \pm 0.01 cm^{-3}$. △ Less

Submitted 17 November, 2023; v1 submitted 18 September, 2023; originally announced September 2023.

Comments: Resubmitted to MNRAS after moderate revision, 12 pages, 6 figures

Accelerated Bayesian imaging by relaxed proximal-point Langevin sampling

Authors: Teresa Klatzer, Paul Dobson, Yoann Altmann, Marcelo Pereyra, Jesús María Sanz-Serna, Konstantinos C. Zygalakis

Abstract: This paper presents a new accelerated proximal Markov chain Monte Carlo methodology to perform Bayesian inference in imaging inverse problems with an underlying convex geometry. The proposed strategy takes the form of a stochastic relaxed proximal-point iteration that admits two complementary interpretations. For models that are smooth or regularised by Moreau-Yosida smoothing, the algorithm is eq… ▽ More This paper presents a new accelerated proximal Markov chain Monte Carlo methodology to perform Bayesian inference in imaging inverse problems with an underlying convex geometry. The proposed strategy takes the form of a stochastic relaxed proximal-point iteration that admits two complementary interpretations. For models that are smooth or regularised by Moreau-Yosida smoothing, the algorithm is equivalent to an implicit midpoint discretisation of an overdamped Langevin diffusion targeting the posterior distribution of interest. This discretisation is asymptotically unbiased for Gaussian targets and shown to converge in an accelerated manner for any target that is $κ$-strongly log-concave (i.e., requiring in the order of $\sqrtκ$ iterations to converge, similarly to accelerated optimisation schemes), comparing favorably to [M. Pereyra, L. Vargas Mieles, K.C. Zygalakis, SIAM J. Imaging Sciences, 13,2 (2020), pp. 905-935] which is only provably accelerated for Gaussian targets and has bias. For models that are not smooth, the algorithm is equivalent to a Leimkuhler-Matthews discretisation of a Langevin diffusion targeting a Moreau-Yosida approximation of the posterior distribution of interest, and hence achieves a significantly lower bias than conventional unadjusted Langevin strategies based on the Euler-Maruyama discretisation. For targets that are $κ$-strongly log-concave, the provided non-asymptotic convergence analysis also identifies the optimal time step which maximizes the convergence speed. The proposed methodology is demonstrated through a range of experiments related to image deconvolution with Gaussian and Poisson noise, with assumption-driven and data-driven convex priors. Source codes for the numerical experiments of this paper are available from https://github.com/MI2G/accelerated-langevin-imla. △ Less

Submitted 12 January, 2024; v1 submitted 18 August, 2023; originally announced August 2023.

Comments: 34 pages, 13 figures

MSC Class: 65C40; 68U10; 62F15; 65C60; 65J22; 68W25

Understanding the Nature of the Optical Emission in Gamma-Ray Bursts: Analysis from TAROT, COATLI, and RATIR Observations

Authors: R. L. Becerra, A. Klotz, J. L. Atteia, D. Guetta, A. M. Watson, F. De Colle, C. Angulo-Valdez, N. R. Butler, S. Dichiara, N. Fraija, K. Garcia-Cifuentes, A. S. Kutyrev, W. H. Lee, M. Pereyra, E. Troja

Abstract: We collected the optical light curve data of 227 gamma-ray bursts (GRBs) observed with the TAROT, COATLI, and RATIR telescopes. These consist of 133 detections and 94 upper limits. We constructed average light curves in the observer and rest frames in both X-rays (from {\itshape Swift}/XRT) and in the optical. Our analysis focused on investigating the observational and intrinsic properties of GRBs… ▽ More We collected the optical light curve data of 227 gamma-ray bursts (GRBs) observed with the TAROT, COATLI, and RATIR telescopes. These consist of 133 detections and 94 upper limits. We constructed average light curves in the observer and rest frames in both X-rays (from {\itshape Swift}/XRT) and in the optical. Our analysis focused on investigating the observational and intrinsic properties of GRBs. Specifically, we examined observational properties, such as the optical brightness function of the GRBs at $T=1000$ seconds after the trigger, as well as the temporal slope of the afterglow. We also estimated the redshift distribution for the GRBs within our sample. Of the 227 GRBs analysed, we found that 116 had a measured redshift. Based on these data, we calculated a local rate of $ρ_0=0.2$ Gpc$^{-3}$ yr$^{-1}$ for these events with $z<1$. To explore the intrinsic properties of GRBs, we examined the average X-ray and optical light curves in the rest frame. We use the {\scshape afterglowpy} library to generate synthetic curves to constrain the parameters typical of the bright GRB jet, such as energy (${\langle} {E_{0}}{\rangle}\sim 10^{53.6}$~erg), opening angle (${\langle}θ_\mathrm{core}{\rangle}\sim 0.2$~rad), and density (${\langle}n_\mathrm{0}{\rangle}\sim10^{-2.1}$ cm$^{-3}$). Furthermore, we analyse microphysical parameters, including the fraction of thermal energy in accelerated electrons (${\langle}ε_e{\rangle}\sim 10^{-1.37}$) and in the magnetic field (${\langle}ε_B{\rangle}\sim10^{-2.26}$), and the power-law index of the population of non-thermal electrons (${\langle}p{\rangle}\sim 2.2$). △ Less

Submitted 17 August, 2023; v1 submitted 15 August, 2023; originally announced August 2023.

Comments: Accepted for publication in MNRAS, 13 pages and 9 figures

Proximal nested sampling with data-driven priors for physical scientists

Authors: Jason D. McEwen, Tobías I. Liaudat, Matthew A. Price, Xiaohao Cai, Marcelo Pereyra

Abstract: 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 fra… ▽ More 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. △ Less

Submitted 28 July, 2023; v1 submitted 30 June, 2023; originally announced July 2023.

Comments: 9 pages, 4 figures

Weighted Inequalities for $t$-Haar multipliers

Authors: Daewon Chung, Weiyan Huang, Jean Carlo Moraes, María Cristina Pereyra, Brett D. Wick

Abstract: In this paper, we provide necessary and sufficient conditions on a triple of weights $(u, v, w)$ so that the $t$-Haar multipliers $T^t_{w,σ}, t \in \mathbb{R}$, are uniformly (on the choice of signs $σ$) bounded from $L^2(u)$ into $L^2(v)$. These dyadic operators have symbols $s(x; I) = σ_I (w(x)/\langle w \rangle_I )^t$ which are functions of the space variable $x \in \mathbb{R}$ and the frequenc… ▽ More In this paper, we provide necessary and sufficient conditions on a triple of weights $(u, v, w)$ so that the $t$-Haar multipliers $T^t_{w,σ}, t \in \mathbb{R}$, are uniformly (on the choice of signs $σ$) bounded from $L^2(u)$ into $L^2(v)$. These dyadic operators have symbols $s(x; I) = σ_I (w(x)/\langle w \rangle_I )^t$ which are functions of the space variable $x \in \mathbb{R}$ and the frequency variable $I \in \mathcal{D}$, making them dyadic analogues of pseudo-differential operators. Here $\mathcal{D}$ denotes the dyadic intervals, $σ_I = \pm 1$, and $\langle w \rangle_I$ denotes the integral average of $w$ on $I$. When $w\equiv 1$ we have the martingale transform and our conditions recover the known two-weight necessary and sufficient conditions of Nazarov, Treil and Volberg. We also show how these conditions are simplified when $u = v$. In particular, the martingale one-weight and the t-Haar multiplier unsigned and unweighted (corresponding to $σ_I= 1$ and $u = v \equiv 1$) known results are recovered or improved. We also obtain necessary and sufficient testing conditions of Sawyer type for the two-weight boundedness of a single variable Haar multiplier similar to those known for the martingale transform. △ Less

Submitted 28 March, 2023; v1 submitted 25 March, 2023; originally announced March 2023.

Deciphering the unusual stellar progenitor of GRB 210704A

Authors: R. L. Becerra, E. Troja, A. M. Watson, B. O'Connor, P. Veres, S. Dichiara, N. R. Butler, T. Sakamoto, K. O. C. Lopez, F. De Colle, K. Aoki, N. Fraija, M. Im, A. S. Kutyrev, W. H. Lee, G. S. H. Paek, M. Pereyra, S. Ravi, Y. Urata

Abstract: GRB~210704A is a burst of intermediate duration ($T_{90} \sim 1-4$~s) followed by a fading afterglow and an optical excess that peaked about 7 days after the explosion. Its properties, and in particular those of the excess, do not easily fit into the well established classification scheme of GRBs as being long or short, leaving the nature of its progenitor uncertain. We present multi-wavelength ob… ▽ More GRB~210704A is a burst of intermediate duration ($T_{90} \sim 1-4$~s) followed by a fading afterglow and an optical excess that peaked about 7 days after the explosion. Its properties, and in particular those of the excess, do not easily fit into the well established classification scheme of GRBs as being long or short, leaving the nature of its progenitor uncertain. We present multi-wavelength observations of the GRB and its counterpart, observed up to 160 days after the burst. In order to decipher the nature of the progenitor system, we present a detailed analysis of the GRB high-energy properties (duration, spectral lag, and Amati correlation), its environment, and late-time optical excess. We discuss three possible scenarios: a neutron star merger, a collapsing massive star, and an atypical explosion possibly hosted in a cluster of galaxies. We find that traditional kilonova and supernova models do not match well the properties of the optical excess, leaving us with the intriguing suggestion that this event was an exotic high-energy merger. △ Less

Submitted 2 May, 2023; v1 submitted 13 March, 2023; originally announced March 2023.

Comments: Revised version submitted to MNRAS after minor comments, 14 pages, 9 figures

The split Gibbs sampler revisited: improvements to its algorithmic structure and augmented target distribution

Authors: Marcelo Pereyra, Luis A. Vargas-Mieles, Konstantinos C. Zygalakis

Abstract: Developing efficient Bayesian computation algorithms for imaging inverse problems is challenging due to the dimensionality involved and because Bayesian imaging models are often not smooth. Current state-of-the-art methods often address these difficulties by replacing the posterior density with a smooth approximation that is amenable to efficient exploration by using Langevin Markov chain Monte Ca… ▽ More Developing efficient Bayesian computation algorithms for imaging inverse problems is challenging due to the dimensionality involved and because Bayesian imaging models are often not smooth. Current state-of-the-art methods often address these difficulties by replacing the posterior density with a smooth approximation that is amenable to efficient exploration by using Langevin Markov chain Monte Carlo (MCMC) methods. An alternative approach is based on data augmentation and relaxation, where auxiliary variables are introduced in order to construct an approximate augmented posterior distribution that is amenable to efficient exploration by Gibbs sampling. This paper proposes a new accelerated proximal MCMC method called latent space SK-ROCK (ls SK-ROCK), which tightly combines the benefits of the two aforementioned strategies. Additionally, instead of viewing the augmented posterior distribution as an approximation of the original model, we propose to consider it as a generalisation of this model. Following on from this, we empirically show that there is a range of values for the relaxation parameter for which the accuracy of the model improves, and propose a stochastic optimisation algorithm to automatically identify the optimal amount of relaxation for a given problem. In this regime, ls SK-ROCK converges faster than competing approaches from the state of the art, and also achieves better accuracy since the underlying augmented Bayesian model has a higher Bayesian evidence. The proposed methodology is demonstrated with a range of numerical experiments related to image deblurring and inpainting, as well as with comparisons with alternative approaches from the state of the art. An open-source implementation of the proposed MCMC methods is available from https://github.com/luisvargasmieles/ls-MCMC. △ Less

Submitted 3 May, 2023; v1 submitted 28 June, 2022; originally announced June 2022.

Comments: 22 pages, 9 figures, 7 tables. Accepted for publication in SIAM Journal on Imaging Sciences (SIIMS)

MSC Class: 60J22; 65C40; 68U10 (Primary) 62E17; 62F15; 62H10; 65J22; 68W25 (Secondary)

Learned reconstruction methods with convergence guarantees

Authors: Subhadip Mukherjee, Andreas Hauptmann, Ozan Öktem, Marcelo Pereyra, Carola-Bibiane Schönlieb

Abstract: In recent years, deep learning has achieved remarkable empirical success for image reconstruction. This has catalyzed an ongoing quest for precise characterization of correctness and reliability of data-driven methods in critical use-cases, for instance in medical imaging. Notwithstanding the excellent performance and efficacy of deep learning-based methods, concerns have been raised regarding the… ▽ More In recent years, deep learning has achieved remarkable empirical success for image reconstruction. This has catalyzed an ongoing quest for precise characterization of correctness and reliability of data-driven methods in critical use-cases, for instance in medical imaging. Notwithstanding the excellent performance and efficacy of deep learning-based methods, concerns have been raised regarding their stability, or lack thereof, with serious practical implications. Significant advances have been made in recent years to unravel the inner workings of data-driven image recovery methods, challenging their widely perceived black-box nature. In this article, we will specify relevant notions of convergence for data-driven image reconstruction, which will form the basis of a survey of learned methods with mathematically rigorous reconstruction guarantees. An example that is highlighted is the role of ICNN, offering the possibility to combine the power of deep learning with classical convex regularization theory for devising methods that are provably convergent. This survey article is aimed at both methodological researchers seeking to advance the frontiers of our understanding of data-driven image reconstruction methods as well as practitioners, by providing an accessible description of useful convergence concepts and by placing some of the existing empirical practices on a solid mathematical foundation. △ Less

Submitted 14 September, 2022; v1 submitted 11 June, 2022; originally announced June 2022.

Efficient Bayesian computation for low-photon imaging problems

Authors: Savvas Melidonis, Paul Dobson, Yoann Altmann, Marcelo Pereyra, Konstantinos C. Zygalakis

Abstract: This paper studies a new and highly efficient Markov chain Monte Carlo (MCMC) methodology to perform Bayesian inference in low-photon imaging problems, with particular attention to situations involving observation noise processes that deviate significantly from Gaussian noise, such as binomial, geometric and low-intensity Poisson noise. These problems are challenging for many reasons. From an infe… ▽ More This paper studies a new and highly efficient Markov chain Monte Carlo (MCMC) methodology to perform Bayesian inference in low-photon imaging problems, with particular attention to situations involving observation noise processes that deviate significantly from Gaussian noise, such as binomial, geometric and low-intensity Poisson noise. These problems are challenging for many reasons. From an inferential viewpoint, low-photon numbers lead to severe identifiability issues, poor stability and high uncertainty about the solution. Moreover, low-photon models often exhibit poor regularity properties that make efficient Bayesian computation difficult; e.g., hard non-negativity constraints, non-smooth priors, and log-likelihood terms with exploding gradients. More precisely, the lack of suitable regularity properties hinders the use of state-of-the-art Monte Carlo methods based on numerical approximations of the Langevin stochastic differential equation (SDE), as both the SDE and its numerical approximations behave poorly. We address this difficulty by proposing an MCMC methodology based on a reflected and regularised Langevin SDE, which is shown to be well-posed and exponentially ergodic under mild and easily verifiable conditions. This then allows us to derive four reflected proximal Langevin MCMC algorithms to perform Bayesian computation in low-photon imaging problems. The proposed approach is demonstrated with a range of experiments related to image deblurring, denoising, and inpainting under binomial, geometric and Poisson noise. △ Less

Submitted 10 June, 2022; originally announced June 2022.

Comments: 28 pages, 10 figures

MSC Class: 65C40; 68U10 (Primary) 62F15; 65C60; 65J22; 62E17; 62F30; 62H10; 68W25 (Secondary)

On Maximum-a-Posteriori estimation with Plug & Play priors and stochastic gradient descent

Authors: Rémi Laumont, Valentin de Bortoli, Andrés Almansa, Julie Delon, Alain Durmus, Marcelo Pereyra

Abstract: Bayesian methods to solve imaging inverse problems usually combine an explicit data likelihood function with a prior distribution that explicitly models expected properties of the solution. Many kinds of priors have been explored in the literature, from simple ones expressing local properties to more involved ones exploiting image redundancy at a non-local scale. In a departure from explicit model… ▽ More Bayesian methods to solve imaging inverse problems usually combine an explicit data likelihood function with a prior distribution that explicitly models expected properties of the solution. Many kinds of priors have been explored in the literature, from simple ones expressing local properties to more involved ones exploiting image redundancy at a non-local scale. In a departure from explicit modelling, several recent works have proposed and studied the use of implicit priors defined by an image denoising algorithm. This approach, commonly known as Plug & Play (PnP) regularisation, can deliver remarkably accurate results, particularly when combined with state-of-the-art denoisers based on convolutional neural networks. However, the theoretical analysis of PnP Bayesian models and algorithms is difficult and works on the topic often rely on unrealistic assumptions on the properties of the image denoiser. This papers studies maximum-a-posteriori (MAP) estimation for Bayesian models with PnP priors. We first consider questions related to existence, stability and well-posedness, and then present a convergence proof for MAP computation by PnP stochastic gradient descent (PnP-SGD) under realistic assumptions on the denoiser used. We report a range of imaging experiments demonstrating PnP-SGD as well as comparisons with other PnP schemes. △ Less

Submitted 16 January, 2022; originally announced January 2022.

MSC Class: 65K10 (Primary) 65K05; 62F15; 62C10; 68Q25; 68U10; 90C26 (Secondary) 65K10; 65K05; 62F15; 62C10; 68Q25; 68U10; 90C26

GRB 191016A: The onset of the forward shock and evidence of late energy injection

Authors: M. Pereyra, N. Fraija, A. M. Watson, R. L. Becerra, N. R. Butler, F. De Colle, E. Troja, S. Dichiara, E. Fraire-Bonilla, W. H. Lee, A. S. Kutyrev, J. X. Prochaska, J. S. Bloom, J. J. González, E. Ramirez-Ruiz, M. G. Richer

Abstract: We present optical and near-infrared photometric observations of GRB 191016 with the COATLI, DDOTI and RATIR ground-based telescopes over the first three nights. We present the temporal evolution of the optical afterglow and describe 5 different stages that were not completely characterized in previous works, mainly due to scarcity of data points to accurately fit the different components of the o… ▽ More We present optical and near-infrared photometric observations of GRB 191016 with the COATLI, DDOTI and RATIR ground-based telescopes over the first three nights. We present the temporal evolution of the optical afterglow and describe 5 different stages that were not completely characterized in previous works, mainly due to scarcity of data points to accurately fit the different components of the optical emission. After the end of the prompt gamma-ray emission, we observed the afterglow rise slowly in the optical and near-infrared (NIR) wavelengths and peak at around T+1450s in all filters. This was followed by an early decay, a clear plateau from T+5000s to T+11000s, and then a regular late decay. We also present evidence of the jet break at later times, with a temporal index in good agreement with the temporal slope obtained from X-ray observations. Although many of the features observed in the optical light curves of GRBs are usually well explained by a reverse shock (RS) or forward shock(FS), the shallowness of the optical rise and enhanced peak emission in the GRB191016A afterglow is not well-fitted by only a FS or a RS. We propose a theoretical model which considers both of these components and combines an evolving FS with a later embedded RS and a subsequent late energy injection from the central engine activity. We use this model to successfully explain the temporal evolution of the light curves and discuss its implications on the fireball properties. △ Less

Submitted 30 November, 2021; originally announced November 2021.

Comments: 15 pages, 5 figures, 7 Tables, Submitted to MNRAS

Constraints on the electromagnetic counterpart of the Neutron Star Black Hole merger GW200115

Authors: S. Dichiara, R. L. Becerra, E. A. Chase, E. Troja, W. H. Lee, A. M. Watson, N. R. Butler, B. O'Connor, M. Pereyra, K. O. C. López, A. Y. Lien, A. Gottlieb, A. S. Kutyrev

Abstract: We report the results of our follow-up campaign for the neutron star - black hole (NSBH) merger GW200115 detected during the O3 run of the Advanced LIGO and Advanced Virgo detectors. We obtained wide-field observations with the Deca-Degree Optical Transient Imager (DDOTI) covering ~20% of the total probability area down to a limiting magnitude of $w$=20.5 AB at ~23 h after the merger. Our search f… ▽ More We report the results of our follow-up campaign for the neutron star - black hole (NSBH) merger GW200115 detected during the O3 run of the Advanced LIGO and Advanced Virgo detectors. We obtained wide-field observations with the Deca-Degree Optical Transient Imager (DDOTI) covering ~20% of the total probability area down to a limiting magnitude of $w$=20.5 AB at ~23 h after the merger. Our search for counterparts returns a single candidate (AT2020aeo), likely not associate to the merger. In total, only 25 sources of interest were identified by the community and later discarded as unrelated to the GW event. We compare our upper limits with the emission predicted by state-of-the-art kilonova simulations and disfavor high mass ejecta (>0.1$M_{\odot}$), indicating that the spin of the system is not particularly high. By combining our optical limits with gamma-ray constraints from $Swift$ and $Fermi$, we disfavor the presence of a standard short duration burst for viewing angles $\lesssim$15 deg from the jet axis. Our conclusions are however limited by the large localization region of this GW event, and accurate prompt positions remain crucial to improving the efficiency of follow-up efforts. △ Less

Submitted 16 February, 2022; v1 submitted 22 October, 2021; originally announced October 2021.

Comments: 6 pages, 4 figures, accepted for publication ApJ

DDOTI Observations of Gravitational-Wave Sources Discovered in O3

Authors: R. L. Becerra, S. Dichiara, A. M. Watson, E. Troja, N. R. Butler, M. Pereyra, E. Moreno Méndez, F. De Colle, W. H. Lee, A. S. Kutyrev, K. O. C. López

Abstract: We present optical follow-up observations with the DDOTI telescope of gravitational-wave events detected during the Advanced LIGO and Advanced Virgo O3 observing run. DDOTI is capable of responding to an alert in a few minutes, has an instantaneous field of about 69 deg$^{2}$, and obtains $10σ$ upper limits of $w_{\rm lim}=18.5$ to 20.5 AB mag in 1000~s of exposure, depending on the conditions. We… ▽ More We present optical follow-up observations with the DDOTI telescope of gravitational-wave events detected during the Advanced LIGO and Advanced Virgo O3 observing run. DDOTI is capable of responding to an alert in a few minutes, has an instantaneous field of about 69 deg$^{2}$, and obtains $10σ$ upper limits of $w_{\rm lim}=18.5$ to 20.5 AB mag in 1000~s of exposure, depending on the conditions. We observed 54\% (26 out of 48) of the unretracted gravitational-wave alerts and did not find any electromagnetic counterparts. We compare our upper limits to various possible counterparts: the kilonova AT~2017gfo, models of radioactive- and magnetar-powered kilonovae, short gamma-ray burst afterglows, and AGN flares. Although the large positional uncertainties of GW sources do not allow us to place strong constraints during O3, DDOTI observations of well-localized GW events in O4 and beyond could meaningfully constrain models of compact binary mergers. We show that DDOTI is able to detect kilonovae similar to AT~2017gfo up to about 200~Mpc and magnetar-powered kilonovae up to 1~Gpc. We calculate that nearby ($\lesssim$200 Mpc) afterglows have a high chance ($\approx$70\%) to be detected by rapid ($\lesssim$3 hours) DDOTI observations if observed on-axis, whereas off-axis afterglows are unlikely to be seen. Finally, we suggest that long-term monitoring of massive BBH events with DDOTI could confirm or rule out late AGN flares associated with these events. △ Less

Submitted 19 July, 2021; v1 submitted 28 June, 2021; originally announced June 2021.

Comments: 26 pages, 5 figures. Accepted for publication in MNRAS

Proximal nested sampling for high-dimensional Bayesian model selection

Authors: Xiaohao Cai, Jason D. McEwen, Marcelo Pereyra

Abstract: 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 applica… ▽ More 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. △ Less

Submitted 9 September, 2022; v1 submitted 7 June, 2021; originally announced June 2021.

Journal ref: Statistics & Computing, 32, 87 (2002)

Bayesian Imaging With Data-Driven Priors Encoded by Neural Networks: Theory, Methods, and Algorithms

Authors: Matthew Holden, Marcelo Pereyra, Konstantinos C. Zygalakis

Abstract: This paper proposes a new methodology for performing Bayesian inference in imaging inverse problems where the prior knowledge is available in the form of training data. Following the manifold hypothesis and adopting a generative modelling approach, we construct a data-driven prior that is supported on a sub-manifold of the ambient space, which we can learn from the training data by using a variati… ▽ More This paper proposes a new methodology for performing Bayesian inference in imaging inverse problems where the prior knowledge is available in the form of training data. Following the manifold hypothesis and adopting a generative modelling approach, we construct a data-driven prior that is supported on a sub-manifold of the ambient space, which we can learn from the training data by using a variational autoencoder or a generative adversarial network. We establish the existence and well-posedness of the associated posterior distribution and posterior moments under easily verifiable conditions, providing a rigorous underpinning for Bayesian estimators and uncertainty quantification analyses. Bayesian computation is performed by using a parallel tempered version of the preconditioned Crank-Nicolson algorithm on the manifold, which is shown to be ergodic and robust to the non-convex nature of these data-driven models. In addition to point estimators and uncertainty quantification analyses, we derive a model misspecification test to automatically detect situations where the data-driven prior is unreliable, and explain how to identify the dimension of the latent space directly from the training data. The proposed approach is illustrated with a range of experiments with the MNIST dataset, where it outperforms alternative image reconstruction approaches from the state of the art. A model accuracy analysis suggests that the Bayesian probabilities reported by the data-driven models are also remarkably accurate under a frequentist definition of probability. △ Less

Submitted 18 March, 2021; originally announced March 2021.

Bayesian imaging using Plug & Play priors: when Langevin meets Tweedie

Authors: Rémi Laumont, Valentin de Bortoli, Andrés Almansa, Julie Delon, Alain Durmus, Marcelo Pereyra

Abstract: Since the seminal work of Venkatakrishnan et al. in 2013, Plug & Play (PnP) methods have become ubiquitous in Bayesian imaging. These methods derive Minimum Mean Square Error (MMSE) or Maximum A Posteriori (MAP) estimators for inverse problems in imaging by combining an explicit likelihood function with a prior that is implicitly defined by an image denoising algorithm. The PnP algorithms proposed… ▽ More Since the seminal work of Venkatakrishnan et al. in 2013, Plug & Play (PnP) methods have become ubiquitous in Bayesian imaging. These methods derive Minimum Mean Square Error (MMSE) or Maximum A Posteriori (MAP) estimators for inverse problems in imaging by combining an explicit likelihood function with a prior that is implicitly defined by an image denoising algorithm. The PnP algorithms proposed in the literature mainly differ in the iterative schemes they use for optimisation or for sampling. In the case of optimisation schemes, some recent works guarantee the convergence to a fixed point, albeit not necessarily a MAP estimate. In the case of sampling schemes, to the best of our knowledge, there is no known proof of convergence. There also remain important open questions regarding whether the underlying Bayesian models and estimators are well defined, well-posed, and have the basic regularity properties required to support these numerical schemes. To address these limitations, this paper develops theory, methods, and provably convergent algorithms for performing Bayesian inference with PnP priors. We introduce two algorithms: 1) PnP-ULA (Unadjusted Langevin Algorithm) for Monte Carlo sampling and MMSE inference; and 2) PnP-SGD (Stochastic Gradient Descent) for MAP inference. Using recent results on the quantitative convergence of Markov chains, we establish detailed convergence guarantees for these two algorithms under realistic assumptions on the denoising operators used, with special attention to denoisers based on deep neural networks. We also show that these algorithms approximately target a decision-theoretically optimal Bayesian model that is well-posed. The proposed algorithms are demonstrated on several canonical problems such as image deblurring, inpainting, and denoising, where they are used for point estimation as well as for uncertainty visualisation and quantification. △ Less

Submitted 12 January, 2022; v1 submitted 8 March, 2021; originally announced March 2021.

MSC Class: 65K10; 65K05; 65D18; 62F15; 62C10; 68Q25; 68U10; 90C26

Bayesian model selection for unsupervised image deconvolution with structured Gaussian priors

Authors: Benjamin Harroué, Jean-François Giovannelli, Marcelo Pereyra

Abstract: This paper considers the objective comparison of stochastic models to solve inverse problems, more specifically image restoration. Most often, model comparison is addressed in a supervised manner, that can be time-consuming and partly arbitrary. Here we adopt an unsupervised Bayesian approach and objectively compare the models based on their posterior probabilities, directly from the data without… ▽ More This paper considers the objective comparison of stochastic models to solve inverse problems, more specifically image restoration. Most often, model comparison is addressed in a supervised manner, that can be time-consuming and partly arbitrary. Here we adopt an unsupervised Bayesian approach and objectively compare the models based on their posterior probabilities, directly from the data without ground truth available. The probabilities depend on the marginal likelihood or "evidence" of the models and we resort to the Chib approach including a Gibbs sampler. We focus on the family of Gaussian models with circulant covariances and unknown hyperparameters, and compare different types of covariance matrices for the image and noise. △ Less

Submitted 13 October, 2020; originally announced October 2020.

Comments: 5 pages,5 figures, conference paper, journal paper in processing

Maximum likelihood estimation of regularisation parameters in high-dimensional inverse problems: an empirical Bayesian approach. Part II: Theoretical Analysis

Authors: Valentin De Bortoli, Alain Durmus, Ana F. Vidal, Marcelo Pereyra

Abstract: This paper presents a detailed theoretical analysis of the three stochastic approximation proximal gradient algorithms proposed in our companion paper [49] to set regularization parameters by marginal maximum likelihood estimation. We prove the convergence of a more general stochastic approximation scheme that includes the three algorithms of [49] as special cases. This includes asymptotic and non… ▽ More This paper presents a detailed theoretical analysis of the three stochastic approximation proximal gradient algorithms proposed in our companion paper [49] to set regularization parameters by marginal maximum likelihood estimation. We prove the convergence of a more general stochastic approximation scheme that includes the three algorithms of [49] as special cases. This includes asymptotic and non-asymptotic convergence results with natural and easily verifiable conditions, as well as explicit bounds on the convergence rates. Importantly, the theory is also general in that it can be applied to other intractable optimisation problems. A main novelty of the work is that the stochastic gradient estimates of our scheme are constructed from inexact proximal Markov chain Monte Carlo samplers. This allows the use of samplers that scale efficiently to large problems and for which we have precise theoretical guarantees. △ Less

Submitted 13 August, 2020; originally announced August 2020.

Comments: SIIMS 2020 - 30 pages

A search for optical and near-infrared counterparts of the compact binary merger GW190814

Authors: A. L. Thakur, S. Dichiara, E. Troja, E. A. Chase, R. Sanchez-Ramirez, L. Piro, C. L. Fryer, N. R. Butler, A. M. Watson, R. T. Wollaeger, E. Ambrosi, J. Becerra González, R. L. Becerra, G. Bruni, S. B. Cenko, G. Cusumano, Antonino D'Aì, J. Durbak, C. J. Fontes, P. Gatkine, A. L. Hungerford, O. Korobkin, A. S. Kutyrev, W. H. Lee, S. Lotti , et al. (7 additional authors not shown)

Abstract: We report on our observing campaign of the compact binary merger GW190814, detected by the Advanced LIGO and Advanced Virgo detectors on August 14th, 2019. This signal has the best localisation of any observed gravitational wave (GW) source, with a 90% probability area of 18.5 deg$^2$, and an estimated distance of ~ 240 Mpc. We obtained wide-field observations with the Deca-Degree Optical Transien… ▽ More We report on our observing campaign of the compact binary merger GW190814, detected by the Advanced LIGO and Advanced Virgo detectors on August 14th, 2019. This signal has the best localisation of any observed gravitational wave (GW) source, with a 90% probability area of 18.5 deg$^2$, and an estimated distance of ~ 240 Mpc. We obtained wide-field observations with the Deca-Degree Optical Transient Imager (DDOTI) covering 88% of the probability area down to a limiting magnitude of $w$ = 19.9 AB. Nearby galaxies within the high probability region were targeted with the Lowell Discovery Telescope (LDT), whereas promising candidate counterparts were characterized through multi-colour photometry with the Reionization and Transients InfraRed (RATIR) and spectroscopy with the Gran Telescopio de Canarias (GTC). We use our optical and near-infrared limits in conjunction with the upper limits obtained by the community to constrain the possible electromagnetic counterparts associated with the merger. A gamma-ray burst seen along its jet's axis is disfavoured by the multi-wavelength dataset, whereas the presence of a burst seen at larger viewing angles is not well constrained. Although our observations are not sensitive to a kilonova similar to AT2017gfo, we can rule out high-mass (> 0.1 M$_{\odot}$) fast-moving (mean velocity >= 0.3c) wind ejecta for a possible kilonova associated with this merger. △ Less

Submitted 3 November, 2020; v1 submitted 9 July, 2020; originally announced July 2020.

Comments: 17 pages, 11 figures, 5 tables; updated acknowledgement section. Accepted for publication in MNRAS (10 September 2020)

Wasserstein Control of Mirror Langevin Monte Carlo

Authors: Kelvin Shuangjian Zhang, Gabriel Peyré, Jalal Fadili, Marcelo Pereyra

Abstract: Discretized Langevin diffusions are efficient Monte Carlo methods for sampling from high dimensional target densities that are log-Lipschitz-smooth and (strongly) log-concave. In particular, the Euclidean Langevin Monte Carlo sampling algorithm has received much attention lately, leading to a detailed understanding of its non-asymptotic convergence properties and of the role that smoothness and lo… ▽ More Discretized Langevin diffusions are efficient Monte Carlo methods for sampling from high dimensional target densities that are log-Lipschitz-smooth and (strongly) log-concave. In particular, the Euclidean Langevin Monte Carlo sampling algorithm has received much attention lately, leading to a detailed understanding of its non-asymptotic convergence properties and of the role that smoothness and log-concavity play in the convergence rate. Distributions that do not possess these regularity properties can be addressed by considering a Riemannian Langevin diffusion with a metric capturing the local geometry of the log-density. However, the Monte Carlo algorithms derived from discretizations of such Riemannian Langevin diffusions are notoriously difficult to analyze. In this paper, we consider Langevin diffusions on a Hessian-type manifold and study a discretization that is closely related to the mirror-descent scheme. We establish for the first time a non-asymptotic upper-bound on the sampling error of the resulting Hessian Riemannian Langevin Monte Carlo algorithm. This bound is measured according to a Wasserstein distance induced by a Riemannian metric ground cost capturing the Hessian structure and closely related to a self-concordance-like condition. The upper-bound implies, for instance, that the iterates contract toward a Wasserstein ball around the target density whose radius is made explicit. Our theory recovers existing Euclidean results and can cope with a wide variety of Hessian metrics related to highly non-flat geometries. △ Less

Submitted 11 February, 2020; originally announced February 2020.

Comments: 22 pages, 2 tables

Limits on the Electromagnetic Counterpart to S190814bv

Authors: Alan M. Watson, Nathaniel R. Butler, William H. Lee, Rosa L. Becerra, Margarita Pereyra, Fernando Angeles, Alejandro Farah, Liliana Figueroa, Diego González-Buitrago, Fernando Quirós, Jaime Ruíz-Díaz-Soto, Carlos Tejada de Vargas, Silvio J. Tinoco, Tanner Wolfram

Abstract: We derive limits on any electromagnetic counterpart to the compact binary merger S190814bv, whose parameters are consistent with the merger of a black hole and a neutron star. We present observations with the new wide-field optical imager DDOTI and also consider Swift/BAT observations reported by Palmer et al. (2019). We show that Swift/BAT would have detected a counterpart with similar properties… ▽ More We derive limits on any electromagnetic counterpart to the compact binary merger S190814bv, whose parameters are consistent with the merger of a black hole and a neutron star. We present observations with the new wide-field optical imager DDOTI and also consider Swift/BAT observations reported by Palmer et al. (2019). We show that Swift/BAT would have detected a counterpart with similar properties to a typical on-axis short GRB at the 98 per cent confidence level, whereas our DDOTI observations only rule out such a counterpart at the 27 per cent confidence level. Neither have sufficient sensitivity to rule out an off-axis counterpart like GW 170817. We compare the efficiency of Swift/BAT and DDOTI for future observations, and show that DDOTI is likely to be about twice as efficient as Swift/BAT for off-axis events up to about 100 Mpc. △ Less

Submitted 20 January, 2020; v1 submitted 15 January, 2020; originally announced January 2020.

Comments: Accepted by MNRAS. 6 pages. Fixed typos and references

Size and Shape Constraints of (486958) Arrokoth from Stellar Occultations

Authors: Marc W. Buie, Simon B. Porter, Peter Tamblyn, Dirk Terrell, Alex Harrison Parker, David Baratoux, Maram Kaire, Rodrigo Leiva, Anne J. Verbiscer, Amanda M. Zangari, François Colas, Baïdy Demba Diop, Joseph I. Samaniego, Lawrence H. Wasserman, Susan D. Benecchi, Amir Caspi, Stephen Gwyn, J. J. Kavelaars, Adriana C. Ocampo Uría, Jorge Rabassa, M. F. Skrutskie, Alejandro Soto, Paolo Tanga, Eliot F. Young, S. Alan Stern , et al. (108 additional authors not shown)

Abstract: We present the results from four stellar occultations by (486958) Arrokoth, the flyby target of the New Horizons extended mission. Three of the four efforts led to positive detections of the body, and all constrained the presence of rings and other debris, finding none. Twenty-five mobile stations were deployed for 2017 June 3 and augmented by fixed telescopes. There were no positive detections fr… ▽ More We present the results from four stellar occultations by (486958) Arrokoth, the flyby target of the New Horizons extended mission. Three of the four efforts led to positive detections of the body, and all constrained the presence of rings and other debris, finding none. Twenty-five mobile stations were deployed for 2017 June 3 and augmented by fixed telescopes. There were no positive detections from this effort. The event on 2017 July 10 was observed by SOFIA with one very short chord. Twenty-four deployed stations on 2017 July 17 resulted in five chords that clearly showed a complicated shape consistent with a contact binary with rough dimensions of 20 by 30 km for the overall outline. A visible albedo of 10% was derived from these data. Twenty-two systems were deployed for the fourth event on 2018 Aug 4 and resulted in two chords. The combination of the occultation data and the flyby results provides a significant refinement of the rotation period, now estimated to be 15.9380 $\pm$ 0.0005 hours. The occultation data also provided high-precision astrometric constraints on the position of the object that were crucial for supporting the navigation for the New Horizons flyby. This work demonstrates an effective method for obtaining detailed size and shape information and probing for rings and dust on distant Kuiper Belt objects as well as being an important source of positional data that can aid in spacecraft navigation that is particularly useful for small and distant bodies. △ Less

Submitted 31 December, 2019; originally announced January 2020.

Comments: Submitted to Astronomical Journal (revised); 40 pages, 13 figures, 9 tables

Journal ref: The Astronomical Journal, Vol. 159, Issue 4, 130 (27pp); 2020 April

Maximum likelihood estimation of regularisation parameters in high-dimensional inverse problems: an empirical Bayesian approach. Part I: Methodology and Experiments

Authors: Ana F. Vidal, Valentin De Bortoli, Marcelo Pereyra, Alain Durmus

Abstract: Many imaging problems require solving an inverse problem that is ill-conditioned or ill-posed. Imaging methods typically address this difficulty by regularising the estimation problem to make it well-posed. This often requires setting the value of the so-called regularisation parameters that control the amount of regularisation enforced. These parameters are notoriously difficult to set a priori,… ▽ More Many imaging problems require solving an inverse problem that is ill-conditioned or ill-posed. Imaging methods typically address this difficulty by regularising the estimation problem to make it well-posed. This often requires setting the value of the so-called regularisation parameters that control the amount of regularisation enforced. These parameters are notoriously difficult to set a priori, and can have a dramatic impact on the recovered estimates. In this work, we propose a general empirical Bayesian method for setting regularisation parameters in imaging problems that are convex w.r.t. the unknown image. Our method calibrates regularisation parameters directly from the observed data by maximum marginal likelihood estimation, and can simultaneously estimate multiple regularisation parameters. Furthermore, the proposed algorithm uses the same basic operators as proximal optimisation algorithms, namely gradient and proximal operators, and it is therefore straightforward to apply to problems that are currently solved by using proximal optimisation techniques. Our methodology is demonstrated with a range of experiments and comparisons with alternative approaches from the literature. The considered experiments include image denoising, non-blind image deconvolution, and hyperspectral unmixing, using synthesis and analysis priors involving the L1, total-variation, total-variation and L1, and total-generalised-variation pseudo-norms. A detailed theoretical analysis of the proposed method is presented in the companion paper arXiv:2008.05793. △ Less

Submitted 14 August, 2020; v1 submitted 26 November, 2019; originally announced November 2019.

Comments: 37 pages - SIIMS 2020

MSC Class: 62C12; 65C40; 68U10; 62F15; 65J20; 65C60; 65J22

Accelerating proximal Markov chain Monte Carlo by using an explicit stabilised method

Authors: Luis Vargas, Marcelo Pereyra, Konstantinos C. Zygalakis

Abstract: We present a highly efficient proximal Markov chain Monte Carlo methodology to perform Bayesian computation in imaging problems. Similarly to previous proximal Monte Carlo approaches, the proposed method is derived from an approximation of the Langevin diffusion. However, instead of the conventional Euler-Maruyama approximation that underpins existing proximal Monte Carlo methods, here we use a st… ▽ More We present a highly efficient proximal Markov chain Monte Carlo methodology to perform Bayesian computation in imaging problems. Similarly to previous proximal Monte Carlo approaches, the proposed method is derived from an approximation of the Langevin diffusion. However, instead of the conventional Euler-Maruyama approximation that underpins existing proximal Monte Carlo methods, here we use a state-of-the-art orthogonal Runge-Kutta-Chebyshev stochastic approximation that combines several gradient evaluations to significantly accelerate its convergence speed, similarly to accelerated gradient optimisation methods. The proposed methodology is demonstrated via a range of numerical experiments, including non-blind image deconvolution, hyperspectral unmixing, and tomographic reconstruction, with total-variation and $\ell_1$-type priors. Comparisons with Euler-type proximal Monte Carlo methods confirm that the Markov chains generated with our method exhibit significantly faster convergence speeds, achieve larger effective sample sizes, and produce lower mean square estimation errors at equal computational budget. △ Less

Submitted 19 March, 2020; v1 submitted 23 August, 2019; originally announced August 2019.

Comments: 28 pages, 13 figures. Accepted for publication in SIAM Journal on Imaging Sciences (SIIMS)

Efficient stochastic optimisation by unadjusted Langevin Monte Carlo. Application to maximum marginal likelihood and empirical Bayesian estimation

Authors: Valentin De Bortoli, Alain Durmus, Marcelo Pereyra, Ana F. Vidal

Abstract: Stochastic approximation methods play a central role in maximum likelihood estimation problems involving intractable likelihood functions, such as marginal likelihoods arising in problems with missing or incomplete data, and in parametric empirical Bayesian estimation. Combined with Markov chain Monte Carlo algorithms, these stochastic optimisation methods have been successfully applied to a wide… ▽ More Stochastic approximation methods play a central role in maximum likelihood estimation problems involving intractable likelihood functions, such as marginal likelihoods arising in problems with missing or incomplete data, and in parametric empirical Bayesian estimation. Combined with Markov chain Monte Carlo algorithms, these stochastic optimisation methods have been successfully applied to a wide range of problems in science and industry. However, this strategy scales poorly to large problems because of methodological and theoretical difficulties related to using high-dimensional Markov chain Monte Carlo algorithms within a stochastic approximation scheme. This paper proposes to address these difficulties by using unadjusted Langevin algorithms to construct the stochastic approximation. This leads to a highly efficient stochastic optimisation methodology with favourable convergence properties that can be quantified explicitly and easily checked. The proposed methodology is demonstrated with three experiments, including a challenging application to high-dimensional statistical audio analysis and a sparse Bayesian logistic regression with random effects problem. △ Less

Submitted 30 May, 2020; v1 submitted 28 June, 2019; originally announced June 2019.

Sequential path planning for a formation of mobile robots with split and merge

Authors: M. Estefanía Pereyra, R. Gastón Araguás, Miroslav Kulich

Abstract: An algorithm for robot formation path planning is presented in this paper. Given a map of the working environment, the algorithm finds a path for a formation taking into account possible split of the formation and its consecutive merge. The key part of the solution works on a graph and sequentially employs an extended version of Dijkstra's graph-based algorithm for multiple robots. It is thus dete… ▽ More An algorithm for robot formation path planning is presented in this paper. Given a map of the working environment, the algorithm finds a path for a formation taking into account possible split of the formation and its consecutive merge. The key part of the solution works on a graph and sequentially employs an extended version of Dijkstra's graph-based algorithm for multiple robots. It is thus deterministic, complete, computationally inexpensive, and finds a solution for a fixed source node to another node in the graph. Moreover, the presented solution is general enough to be incorporated into high-level tasks like cooperative surveillance and it can benefit from state-of-the-art formation motion planning approaches, which can be used for evaluation of edges of an input graph. The performed experimental results demonstrate the behavior of the method in complex environments for formations consisting of tens of robots. △ Less

Submitted 22 January, 2019; originally announced January 2019.

Comments: 2017 IEEE Latin American Conference on Computational Intelligence (LA-CCI). arXiv admin note: substantial text overlap with arXiv:1901.07408

Sparse Bayesian mass-mapping with uncertainties: local credible intervals

Authors: Matthew A. Price, Xiaohao Cai, Jason D. McEwen, Marcelo Pereyra, Thomas D. Kitching

Abstract: Until recently mass-mapping techniques for weak gravitational lensing convergence reconstruction have lacked a principled statistical framework upon which to quantify reconstruction uncertainties, without making strong assumptions of Gaussianity. In previous work we presented a sparse hierarchical Bayesian formalism for convergence reconstruction that addresses this shortcoming. Here, we draw on t… ▽ More Until recently mass-mapping techniques for weak gravitational lensing convergence reconstruction have lacked a principled statistical framework upon which to quantify reconstruction uncertainties, without making strong assumptions of Gaussianity. In previous work we presented a sparse hierarchical Bayesian formalism for convergence reconstruction that addresses this shortcoming. Here, we draw on the concept of local credible intervals (cf. Bayesian error bars) as an extension of the uncertainty quantification techniques previously detailed. These uncertainty quantification techniques are benchmarked against those recovered via Px-MALA - a state of the art proximal Markov Chain Monte Carlo (MCMC) algorithm. We find that typically our recovered uncertainties are everywhere conservative, of similar magnitude and highly correlated (Pearson correlation coefficient $\geq 0.85$) with those recovered via Px-MALA. Moreover, we demonstrate an increase in computational efficiency of $\mathcal{O}(10^6)$ when using our sparse Bayesian approach over MCMC techniques. This computational saving is critical for the application of Bayesian uncertainty quantification to large-scale stage IV surveys such as LSST and Euclid. △ Less

Submitted 5 February, 2021; v1 submitted 10 December, 2018; originally announced December 2018.

Dyadic harmonic analysis and weighted inequalities: the sparse revolution

Authors: María Cristina Pereyra

Abstract: We will introduce the basics of dyadic harmonic analysis and how it can be used to obtain weighted estimates for classical Calderón-Zygmund singular integral operators and their commutators. Harmonic analysts have used dyadic models for many years as a first step towards the understanding of more complex continuous operators. In 2000 Stefanie Petermichl discovered a representation formula for the… ▽ More We will introduce the basics of dyadic harmonic analysis and how it can be used to obtain weighted estimates for classical Calderón-Zygmund singular integral operators and their commutators. Harmonic analysts have used dyadic models for many years as a first step towards the understanding of more complex continuous operators. In 2000 Stefanie Petermichl discovered a representation formula for the venerable Hilbert transform as an average (over grids) of dyadic shift operators, allowing her to reduce arguments to finding estimates for these simpler dyadic models. For the next decade the technique used to get sharp weighted inequalities was the Bellman function method introduced by Nazarov, Treil, and Volberg, paired with sharp extrapolation by Dragičević et al. Other methods where introduced by Hytönen, Lerner, Cruz-Uribe, Martell, Pérez, Lacey, Reguera, Sawyer, Uriarte-Tuero, involving stopping time and median oscillation arguments, precursors of the very successful domination by positive sparse operators methodology. The culmination of this work was Tuomas Hytönen's 2012 proof of the $A_2$ conjecture based on a representation formula for any Calderón-Zygmund operator as an average of appropriate dyadic operators. Since then domination by sparse dyadic operators has taken central stage and has found applications well beyond Hytönen's $A_p$ theorem. We will survey this remarkable progression and more in these lecture notes. △ Less

Submitted 3 December, 2018; originally announced December 2018.

Comments: Lecture notes for a mini course taught by the author at the CIMPA 2017 Research School - IX Escuela Santaló: Harmonic Analysis, Geometric Measure Theory and Applications, held in Buenos Aires, Argentina. 65 pages, 3 figures

MSC Class: 42

Quantifying Uncertainty in High Dimensional Inverse Problems by Convex Optimisation

Authors: Xiaohao Cai, Marcelo Pereyra, Jason D. McEwen

Abstract: 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-… ▽ More 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. △ Less

Submitted 5 September, 2019; v1 submitted 4 November, 2018; originally announced November 2018.

Comments: 5 pages, 5 figures

Journal ref: EUSIPCO 2019

Atomic decomposition of product Hardy spaces via wavelet bases on spaces of homogeneous type

Authors: Yongsheng Han, Ji Li, M. Cristina Pereyra, Lesley A. Ward

Abstract: We provide an atomic decomposition of the product Hardy spaces $H^p(\widetilde{X})$ which were recently developed by Han, Li, and Ward in the setting of product spaces of homogeneous type $\widetilde{X} = X_1 \times X_2$. Here each factor $(X_i,d_i,μ_i)$, for $i = 1$, $2$, is a space of homogeneous type in the sense of Coifman and Weiss. These Hardy spaces make use of the orthogonal wavelet base… ▽ More We provide an atomic decomposition of the product Hardy spaces $H^p(\widetilde{X})$ which were recently developed by Han, Li, and Ward in the setting of product spaces of homogeneous type $\widetilde{X} = X_1 \times X_2$. Here each factor $(X_i,d_i,μ_i)$, for $i = 1$, $2$, is a space of homogeneous type in the sense of Coifman and Weiss. These Hardy spaces make use of the orthogonal wavelet bases of Auscher and Hytönen and their underlying reference dyadic grids. However, no additional assumptions on the quasi-metric or on the doubling measure for each factor space are made. To carry out this program, we introduce product $(p,q)$-atoms on $\widetilde{X}$ and product atomic Hardy spaces $H^{p,q}_{\rm at}(\widetilde{X})$. As consequences of the atomic decomposition of $H^p(\widetilde{X})$, we show that for all $q > 1$ the product atomic Hardy spaces coincide with the product Hardy spaces, and we show that the product Hardy spaces are independent of the particular choices of both the wavelet bases and the reference dyadic grids. Likewise, the product Carleson measure spaces ${\rm CMO}^p(\widetilde{X})$, the bounded mean oscillation space ${\rm BMO}(\widetilde{X})$, and the vanishing mean oscillation space ${\rm VMO}(\widetilde{X})$, as defined by Han, Li, and Ward, are also independent of the particular choices of both wavelets and reference dyadic grids. △ Less

Submitted 18 October, 2018; v1 submitted 8 October, 2018; originally announced October 2018.

Comments: Minor typos corrected, reference added. Submitted

MSC Class: 42B35

Scalable Bayesian uncertainty quantification in imaging inverse problems via convex optimization

Authors: Audrey Repetti, Marcelo Pereyra, Yves Wiaux

Abstract: We propose a Bayesian uncertainty quantification method for large-scale imaging inverse problems. Our method applies to all Bayesian models that are log-concave, where maximum-a-posteriori (MAP) estimation is a convex optimization problem. The method is a framework to analyse the confidence in specific structures observed in MAP estimates (e.g., lesions in medical imaging, celestial sources in ast… ▽ More We propose a Bayesian uncertainty quantification method for large-scale imaging inverse problems. Our method applies to all Bayesian models that are log-concave, where maximum-a-posteriori (MAP) estimation is a convex optimization problem. The method is a framework to analyse the confidence in specific structures observed in MAP estimates (e.g., lesions in medical imaging, celestial sources in astronomical imaging), to enable using them as evidence to inform decisions and conclusions. Precisely, following Bayesian decision theory, we seek to assert the structures under scrutiny by performing a Bayesian hypothesis test that proceeds as follows: firstly, it postulates that the structures are not present in the true image, and then seeks to use the data and prior knowledge to reject this null hypothesis with high probability. Computing such tests for imaging problems is generally very difficult because of the high dimensionality involved. A main feature of this work is to leverage probability concentration phenomena and the underlying convex geometry to formulate the Bayesian hypothesis test as a convex problem, that we then efficiently solve by using scalable optimization algorithms. This allows scaling to high-resolution and high-sensitivity imaging problems that are computationally unaffordable for other Bayesian computation approaches. We illustrate our methodology, dubbed BUQO (Bayesian Uncertainty Quantification by Optimization), on a range of challenging Fourier imaging problems arising in astronomy and medicine. △ Less

Submitted 6 November, 2018; v1 submitted 2 March, 2018; originally announced March 2018.

MSC Class: 62F03; 62F15; 49N45; 68U10

Uncertainty quantification for radio interferometric imaging: II. MAP estimation

Authors: Xiaohao Cai, Marcelo Pereyra, Jason D. McEwen

Abstract: 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… ▽ More 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. △ Less

Submitted 11 September, 2018; v1 submitted 13 November, 2017; originally announced November 2017.

Comments: 13 pages, 10 figures, see companion article in this arXiv listing

Journal ref: Monthly Notices of the Royal Astronomical Society, Volume 480, Issue 3, 1 November 2018, Pages 4170--4182

Uncertainty quantification for radio interferometric imaging: I. proximal MCMC methods

Authors: Xiaohao Cai, Marcelo Pereyra, Jason D. McEwen

Abstract: 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 o… ▽ More 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. △ Less

Submitted 11 September, 2018; v1 submitted 13 November, 2017; originally announced November 2017.

Comments: 16 pages, 7 figures, see companion article in this arXiv listing

Journal ref: Monthly Notices of the Royal Astronomical Society, Volume 480, Issue 3, 1 November 2018, Pages 4154--4169

Sampling from a log-concave distribution with compact support with proximal Langevin Monte Carlo

Authors: Nicolas Brosse, Alain Durmus, Éric Moulines, Marcelo Pereyra

Abstract: This paper presents a detailed theoretical analysis of the Langevin Monte Carlo sampling algorithm recently introduced in Durmus et al. (Efficient Bayesian computation by proximal Markov chain Monte Carlo: when Langevin meets Moreau, 2016) when applied to log-concave probability distributions that are restricted to a convex body $\mathsf{K}$. This method relies on a regularisation procedure involv… ▽ More This paper presents a detailed theoretical analysis of the Langevin Monte Carlo sampling algorithm recently introduced in Durmus et al. (Efficient Bayesian computation by proximal Markov chain Monte Carlo: when Langevin meets Moreau, 2016) when applied to log-concave probability distributions that are restricted to a convex body $\mathsf{K}$. This method relies on a regularisation procedure involving the Moreau-Yosida envelope of the indicator function associated with $\mathsf{K}$. Explicit convergence bounds in total variation norm and in Wasserstein distance of order $1$ are established. In particular, we show that the complexity of this algorithm given a first order oracle is polynomial in the dimension of the state space. Finally, some numerical experiments are presented to compare our method with competing MCMC approaches from the literature. △ Less

Submitted 24 May, 2017; originally announced May 2017.

An Atlas of Exotic Variability in IGR J17091-3624: A Comparison with GRS 1915+105

Authors: James Matthew Christopher Court, Diego Altamirano, Margarita Pereyra, Christopher M. Boon, Kazutaka Yamaoka, Tomaso Belloni, Rudy Wijnands, Mayukh Pahari

Abstract: We performed an analysis of all RXTE observations of the Low Mass X-ray Binary and Black Hole Candidate IGR J17091-3624 during the 2011-2013 outburst of the source. By creating lightcurves, hardness-intensity diagrams and power density spectra of each observation, we have created a set of 9 variability `classes' that phenomenologically describe the range of types of variability seen in this object… ▽ More We performed an analysis of all RXTE observations of the Low Mass X-ray Binary and Black Hole Candidate IGR J17091-3624 during the 2011-2013 outburst of the source. By creating lightcurves, hardness-intensity diagrams and power density spectra of each observation, we have created a set of 9 variability `classes' that phenomenologically describe the range of types of variability seen in this object. We compare our set of variability classes to those established by Belloni et al. (2000) to describe the similar behaviour of the LMXB GRS 1915+105, finding that some types of variability seen in IGR J17091-3624 are not represented in data of GRS 1915+105. We also use all available X-ray data of the 2011-2013 outburst of IGR J17091-3624 to analyse its long-term evolution, presenting the first detection of IGR J17091-3624 above 150 keV as well as noting the presence of `re-flares' during the latter stages of the outburst. Using our results we place new constraints on the mass and distance of the object, and find that it accretes at <33% of its Eddington limit. As such, we conclude that Eddington-limited accretion can no longer be considered a sufficient or necessary criterion for GRS 1915+105-like variability to occur in Low Mass X-Ray Binaries. △ Less

Submitted 28 March, 2017; originally announced March 2017.

Comments: 26 Pages, 31 Figures, 8 Tables. Accepted to MNRAS

Efficient Bayesian computation by proximal Markov chain Monte Carlo: when Langevin meets Moreau

Authors: Alain Durmus, Eric Moulines, Marcelo Pereyra

Abstract: Modern imaging methods rely strongly on Bayesian inference techniques to solve challenging imaging problems. Currently, the predominant Bayesian computation approach is convex optimisation, which scales very efficiently to high dimensional image models and delivers accurate point estimation results. However, in order to perform more complex analyses, for example image uncertainty quantification or… ▽ More Modern imaging methods rely strongly on Bayesian inference techniques to solve challenging imaging problems. Currently, the predominant Bayesian computation approach is convex optimisation, which scales very efficiently to high dimensional image models and delivers accurate point estimation results. However, in order to perform more complex analyses, for example image uncertainty quantification or model selection, it is necessary to use more computationally intensive Bayesian computation techniques such as Markov chain Monte Carlo methods. This paper presents a new and highly efficient Markov chain Monte Carlo methodology to perform Bayesian computation for high dimensional models that are log-concave and non-smooth, a class of models that is central in imaging sciences. The methodology is based on a regularised unadjusted Langevin algorithm that exploits tools from convex analysis, namely Moreau-Yoshida envelopes and proximal operators, to construct Markov chains with favourable convergence properties. In addition to scaling efficiently to high dimensions, the method is straightforward to apply to models that are currently solved by using proximal optimisation algorithms. We provide a detailed theoretical analysis of the proposed methodology, including asymptotic and non-asymptotic convergence results with easily verifiable conditions, and explicit bounds on the convergence rates. The proposed methodology is demonstrated with four experiments related to image deconvolution and tomographic reconstruction with total-variation and $\ell_1$ priors, where we conduct a range of challenging Bayesian analyses related to uncertainty quantification, hypothesis testing, and model selection in the absence of ground truth. △ Less

Submitted 22 December, 2016; originally announced December 2016.

Revisiting maximum-a-posteriori estimation in log-concave models

Authors: Marcelo Pereyra

Abstract: Maximum-a-posteriori (MAP) estimation is the main Bayesian estimation methodology in imaging sciences, where high dimensionality is often addressed by using Bayesian models that are log-concave and whose posterior mode can be computed efficiently by convex optimisation. Despite its success and wide adoption, MAP estimation is not theoretically well understood yet. The prevalent view in the communi… ▽ More Maximum-a-posteriori (MAP) estimation is the main Bayesian estimation methodology in imaging sciences, where high dimensionality is often addressed by using Bayesian models that are log-concave and whose posterior mode can be computed efficiently by convex optimisation. Despite its success and wide adoption, MAP estimation is not theoretically well understood yet. The prevalent view in the community is that MAP estimation is not proper Bayesian estimation in a decision-theoretic sense because it does not minimise a meaningful expected loss function (unlike the minimum mean squared error (MMSE) estimator that minimises the mean squared loss). This paper addresses this theoretical gap by presenting a decision-theoretic derivation of MAP estimation in Bayesian models that are log-concave. A main novelty is that our analysis is based on differential geometry, and proceeds as follows. First, we use the underlying convex geometry of the Bayesian model to induce a Riemannian geometry on the parameter space. We then use differential geometry to identify the so-called natural or canonical loss function to perform Bayesian point estimation in that Riemannian manifold. For log-concave models, this canonical loss is the Bregman divergence associated with the negative log posterior density. We then show that the MAP estimator is the only Bayesian estimator that minimises the expected canonical loss, and that the posterior mean or MMSE estimator minimises the dual canonical loss. We also study the question of MAP and MSSE estimation performance in large scales and establish a universal bound on the expected canonical error as a function of dimension, offering new insights into the good performance observed in convex problems. These results provide a new understanding of MAP and MMSE estimation in log-concave settings, and of the multiple roles that convex geometry plays in imaging problems. △ Less

Submitted 18 January, 2019; v1 submitted 19 December, 2016; originally announced December 2016.

Comments: Accepted for publication in SIAM Imaging Sciences

Maximum-a-posteriori estimation with Bayesian confidence regions

Authors: Marcelo Pereyra

Abstract: Solutions to inverse problems that are ill-conditioned or ill-posed may have significant intrinsic uncertainty. Unfortunately, analysing and quantifying this uncertainty is very challenging, particularly in high-dimensional problems. As a result, while most modern mathematical imaging methods produce impressive point estimation results, they are generally unable to quantify the uncertainty in the… ▽ More Solutions to inverse problems that are ill-conditioned or ill-posed may have significant intrinsic uncertainty. Unfortunately, analysing and quantifying this uncertainty is very challenging, particularly in high-dimensional problems. As a result, while most modern mathematical imaging methods produce impressive point estimation results, they are generally unable to quantify the uncertainty in the solutions delivered. This paper presents a new general methodology for approximating Bayesian high-posterior-density credibility regions in inverse problems that are convex and potentially very high-dimensional. The approximations are derived by using recent concentration of measure results related to information theory for log-concave random vectors. A remarkable property of the approximations is that they can be computed very efficiently, even in large-scale problems, by using standard convex optimisation techniques. In particular, they are available as a by-product in problems solved by maximum-a-posteriori estimation. The approximations also have favourable theoretical properties, namely they outer-bound the true high-posterior-density credibility regions, and they are stable with respect to model dimension. The proposed methodology is illustrated on two high-dimensional imaging inverse problems related to tomographic reconstruction and sparse deconvolution, where the approximations are used to perform Bayesian hypothesis tests and explore the uncertainty about the solutions, and where proximal Markov chain Monte Carlo algorithms are used as benchmark to compute exact credible regions and measure the approximation error. △ Less

Submitted 11 July, 2016; v1 submitted 27 February, 2016; originally announced February 2016.

On two weight estimates for dyadic operators

Authors: Oleksandra Beznosova, Daewon Chung, Jean Carlo Moraes, Maria Cristina Pereyra

Abstract: We provide a quantitative two weight estimate for the dyadic paraproduct $π_b$ under certain conditions on a pair of weights $(u;v)$ and $b$ in $Carl_{u,v}$, a new class of functions that we show coincides with BMO when $u = v \in A^d_2$. We obtain quantitative two weight estimates for the dyadic square function and the martingale transforms under the assumption that the maximal function is bounde… ▽ More We provide a quantitative two weight estimate for the dyadic paraproduct $π_b$ under certain conditions on a pair of weights $(u;v)$ and $b$ in $Carl_{u,v}$, a new class of functions that we show coincides with BMO when $u = v \in A^d_2$. We obtain quantitative two weight estimates for the dyadic square function and the martingale transforms under the assumption that the maximal function is bounded from $L_2(u)$ into $L_2(v)$ and $v \in RH^d_1$. Finally we obtain a quantitative two weight estimate from $L_2(u)$ into $L_2(v)$ for the dyadic square function under the assumption that the pair $(u; v)$ is in joint $A^d_2$ and $u^{-1} \in RH^d_1$, this is sharp in the sense that when $u = v$ the conditions reduce to $u \in A^d_2$ and the estimate is the known linear mixed estimate. △ Less

Submitted 5 February, 2016; originally announced February 2016.

The kinematics of the nebular shells around low mass progenitors of PNe with low metallicity

Authors: Margarita Pereyra, José Alberto López, Michael G. Richer

Abstract: We analyze the internal kinematics of 26 Planetary Nebulae (PNe) with low metallicity that appear to derive from progenitor stars of the lowest masses, including the halo PN population. Based upon spatially-resolved, long-slit, echelle spectroscopy drawn from the San Pedro Mártir Kinematic Catalogue of PNe (López et al. 2012), we characterize the kinematics of these PNe measuring their global expa… ▽ More We analyze the internal kinematics of 26 Planetary Nebulae (PNe) with low metallicity that appear to derive from progenitor stars of the lowest masses, including the halo PN population. Based upon spatially-resolved, long-slit, echelle spectroscopy drawn from the San Pedro Mártir Kinematic Catalogue of PNe (López et al. 2012), we characterize the kinematics of these PNe measuring their global expansion velocities based upon the largest sample used to date for this purpose. We find kinematics that follow the trends observed and predicted in other studies, but also find that most of the PNe studied here tend to have expansion velocities less than 20 km/s in all of the emission lines considered. The low expansion velocities that we observe in this sample of low metallicity planetary nebulae with low mass progenitors are most likely a consequence of a weak central star wind driving the kinematics of the nebular shell. This study complements previous results (Pereyra et al. 2013, and references therein) that link the expansion velocities of the PN shells with the characteristics of the central star. △ Less

Submitted 19 January, 2016; originally announced January 2016.

Comments: 26 Pages, 5 Figures, 3 Tables. Accepted for publication in The Astronomical Journal on 14 Dec 2015

Haar bases on quasi-metric measure spaces, and dyadic structure theorems for function spaces on product spaces of homogeneous type

Authors: Anna Kairema, Ji Li, M. Cristina Pereyra, Lesley Ward

Abstract: We give an explicit construction of Haar functions associated to a system of dyadic cubes in a geometrically doubling quasi-metric space equipped with a positive Borel measure, and show that these Haar functions form a basis for $L^p$. Next we focus on spaces $X$ of homogeneous type in the sense of Coifman and Weiss, where we use these Haar functions to define a discrete square function, and hence… ▽ More We give an explicit construction of Haar functions associated to a system of dyadic cubes in a geometrically doubling quasi-metric space equipped with a positive Borel measure, and show that these Haar functions form a basis for $L^p$. Next we focus on spaces $X$ of homogeneous type in the sense of Coifman and Weiss, where we use these Haar functions to define a discrete square function, and hence to define dyadic versions of the function spaces $H^1(X)$ and ${\rm BMO}(X)$. In the setting of product spaces $\widetilde{X} = X_1 \times \cdots \times X_n$ of homogeneous type, we show that the space ${\rm BMO}(\widetilde{X})$ of functions of bounded mean oscillation on $\widetilde{X}$ can be written as the intersection of finitely many dyadic ${\rm BMO}$ spaces on $\widetilde{X}$, and similarly for $A_p(\widetilde{X})$, reverse-Hölder weights on $\widetilde{X}$, and doubling weights on $\widetilde{X}$. We also establish that the Hardy space $H^1(\widetilde{X})$ is a sum of finitely many dyadic Hardy spaces on $\widetilde{X}$, and that the strong maximal function on $\widetilde{X}$ is pointwise comparable to the sum of finitely many dyadic strong maximal functions. These dyadic structure theorems generalize, to product spaces of homogeneous type, the earlier Euclidean analogues for ${\rm BMO}$ and $H^1$ due to Mei and to Li, Pipher and Ward. △ Less

Submitted 12 September, 2015; originally announced September 2015.

Comments: 37 pages, submitted

MSC Class: 42; 43; 46