Abstract. High-dimensional inverse problems present a challenge for Markov chain Monte Carlo (MCMC)-type sampling schemes. Typically, they rely on finding an efficient proposal distribution, which can be difficult for large-scale problems, even with adaptive approaches. Moreover, the autocorrelations of the samples typically increase with dimension, which leads to the need for long sample chains. We present an alternative method for sampling from posterior distributions in nonlinear inverse problems, when the measurement error and prior are both Gaussian. The approach computes a candidate sample by solving a stochastic optimization problem. In the linear case, these samples are directly from the posterior density, but this is not so in the nonlinear case. We derive the form of the sample density in the nonlinear case, and then show how to use it within both a Metropolis-Hastings and importance sampling framework to obtain samples from the posterior distribution of the parameters. We demonstrate, with various small-and medium-scale problems, that randomize-then-optimize can be efficient compared to standard adaptive MCMC algorithms.

Abstract. The connection between Bayesian statistics and the technique of regularization for inverse problems has been given significant attention in recent years. For example, Bayes' Law is frequently used as motivation for variational regularization methods of Tikhonov type. In this setting, the regularization function corresponds to the negative-log of the prior probability density; the fit-to-data function corresponds to the negative-log of the likelihood; and the regularized solution corresponds to the maximizer of the posterior density function, known as the maximum a posteriori (MAP) estimator of the unknown, which in our case is an image. Much of the work in this direction has focused on the development of techniques for efficient computation of MAP estimators (or regularized solutions). Less explored in the inverse problems community, and of interest to us in this paper, is the problem of sampling from the posterior density. To do this, we use a Markov chain Monte Carlo (MCMC) method which has previously appeared in the Bayesian statistics literature, is straightforward to implement, and provides a means of both estimation and uncertainty quantification for the unknown. Additionally, we show how to use the preconditioned conjugate gradient method to compute image samples in cases where direct methods are not feasible. And finally, the MCMC method provides samples of the noise and prior precision (inverse-variance) parameters, which makes regularization parameter selection unnecessary. We focus on linear models with independent and identically distributed Gaussian noise and define the prior using a Gaussian Markov random field. For our numerical experiments, we consider test-cases from both image deconvolution and computed tomography, and our results show that the approach is effective and surprisingly computationally efficient, even in large-scale cases.

In image processing applications, image intensity is often measured via the counting of incident photons emitted by the object of interest. In such cases, image data-noise is accurately modeled by a Poisson distribution. This motivates the use of Poisson maximum likelihood estimation for image reconstruction. However, when the underlying model equation is ill-posed, regularization is needed. Regularized Poisson likelihood estimation has been studied extensively by the authors, though a problem of high importance remains: the choice of the regularization parameter. We will present three statistically motivated methods for choosing the regularization parameter, and numerical examples will be presented to illustrate their effectiveness.

We consider the problem of solving ill-conditioned linear systems Ax = b subject to the nonnegativity constraint x ≥ 0, and in which the vector b is a realization of a random vector b, i.e. b is noisy. We explore what the statistical literature tells us about solving noisy linear systems; we discuss the effect that a substantial black background in the astronomical object being viewed has on the underlying mathematical and statistical models; and, finally, we present several covariance-based preconditioned iterative methods that incorporate this information. Each of the methods presented can be viewed as an implementation of a preconditioned modified residual-norm steepest descent algorithm with a specific preconditioner, and we show that, in fact, the well-known and often used Richardson-Lucy algorithm is one such method. Ill-conditioning can inhibit the ability to take advantage of a priori statistical knowledge, in which case a more traditional preconditioning approach may be appropriate. We briefly discuss this traditional approach as well. Examples from astronomical imaging are used to illustrate concepts and to test and compare algorithms.

We consider a large-scale convex minimization problem with nonnegativity constraints that arises in astronomical imaging. We develop a cost functional which incorporates the statistics of the noise in the image data and Tikhonov regularization to induce stability. We introduce an efficient hybrid gradient projection-reduced Newton (active set) method. By "reduced Newton" we mean taking Newton steps only in the inactive variables. Due to the large size of our problem, we compute approximate reduced Newton steps using conjugate gradient (CG) iteration. We also introduce a highly effective sparse preconditioner that dramatically speeds up CG convergence. A numerical comparison between our method and other standard large-scale constrained minimization algorithms is presented.