Search

The purpose of this technical note is to summarize the relationship between the marginal variance and correlation length of a Gaussian random field with Mat\'ern covariance and the coefficients of the corresponding partial-differential-equation (PDE)-based precision operator.

Accurate spatiotemporal image reconstruction methods are needed for a wide range of biomedical research areas but face challenges due to data incompleteness and computational burden. Data incompleteness arises from the undersampling often required to increase frame rates and reduce acquisition times, while computational burden emerges due to the memory footprint of high-resolution images with three spatial dimensions and extended time horizons. Neural fields, an emerging class of neural networks that act as continuous representations of spatiotemporal objects, have previously been introduced to solve these dynamic imaging problems by reframing image reconstruction to a problem of estimating network parameters. Neural fields can address the twin challenges of data incompleteness and computational burden by exploiting underlying redundancies in these spatiotemporal objects. This work proposes ProxNF, a novel neural field training approach for spatiotemporal image reconstruction leveraging proximal splitting methods to separate computations involving the imaging operator from updates of the network parameter. Specifically, ProxNF evaluates the (subsampled) gradient of the data-fidelity term in the image domain and uses a fully supervised learning approach to update the neural field parameters. By reducing the memory footprint and the computational cost of evaluating the imaging operator, the proposed ProxNF approach allows for reconstructing large, high-resolution spatiotemporal images. This method is demonstrated in two numerical studies involving virtual dynamic contrast-enhanced photoacoustic computed tomography imaging of an anatomically realistic dynamic numerical mouse phantom and a two-compartment model of tumor perfusion., Comment: 13 pages, 10 figures

The spherical Radon transform (SRT) is an integral transform that maps a function to its integrals over concentric spherical shells centered at specified sensor locations. It has several imaging applications, including synthetic aperture radar and photoacoustic computed tomography. However, computation of the SRT can be expensive. Efficient implementation of SRT on general purpose graphic processing units (GPGPUs) often utilizes non-matched implementation of the adjoint operator, leading to inconsistent gradients in optimization-based image reconstruction methods. This work details an efficient implementation of the SRT and its adjoint for the case of a cylindrical measurement aperture. Exploiting symmetry of the cylindrical geometry, the SRT can then be expressed as the composition of two circular Radon transforms (CRT). Utilizing this formulation then allows for an efficient implementation of the SRT as a discrete-to-discrete operator utilizing sparse matrix representation., Comment: 4 pages

Uncertainty Quantification (UQ) is vital to safety-critical model-based analyses, but the widespread adoption of sophisticated UQ methods is limited by technical complexity. In this paper, we introduce UM-Bridge (the UQ and Modeling Bridge), a high-level abstraction and software protocol that facilitates universal interoperability of UQ software with simulation codes. It breaks down the technical complexity of advanced UQ applications and enables separation of concerns between experts. UM-Bridge democratizes UQ by allowing effective interdisciplinary collaboration, accelerating the development of advanced UQ methods, and making it easy to perform UQ analyses from prototype to High Performance Computing (HPC) scale. In addition, we present a library of ready-to-run UQ benchmark problems, all easily accessible through UM-Bridge. These benchmarks support UQ methodology research, enabling reproducible performance comparisons. We demonstrate UM-Bridge with several scientific applications, harnessing HPC resources even using UQ codes not designed with HPC support., Comment: Fix broken reference, add R. Akers as author for contribution to tritium benchmark

Ultrasound computed tomography (USCT) is actively being developed to quantify acoustic tissue properties such as the speed-of-sound (SOS). Although full-waveform inversion (FWI) is an effective method for accurate SOS reconstruction, it can be computationally challenging for large-scale problems. Deep learning-based image-to-image learned reconstruction (IILR) methods are being investigated as scalable and computationally efficient alternatives. This study investigates the impact of the chosen input modalities on IILR methods for high-resolution SOS reconstruction in USCT. The selected modalities are traveltime tomography (TT) and reflection tomography (RT), which produce a low-resolution SOS map and a reflectivity map, respectively. These modalities have been chosen for their lower computational cost relative to FWI and their capacity to provide complementary information: TT offers a direct -- while low resolution -- SOS measure, while RT reveals tissue boundary information. Systematic analyses were facilitated by employing a stylized USCT imaging system with anatomically realistic numerical breast phantoms. Within this testbed, a supervised convolutional neural network (CNN) was trained to map dual-channel (TT and RT images) to a high-resolution SOS map. Moreover, the CNN was fine-tuned using a weighted reconstruction loss that prioritized tumor regions to address tumor underrepresentation in the training dataset. To understand the benefits of employing dual-channel inputs, single-input CNNs were trained separately using inputs from each modality alone (TT or RT). The methods were assessed quantitatively using normalized root mean squared error and structural similarity index measure for reconstruction accuracy and receiver operating characteristic analysis to assess signal detection-based performance measures.

Significance: Dynamic photoacoustic computed tomography (PACT) is a valuable technique for monitoring physiological processes. However, current dynamic PACT techniques are often limited to 2D spatial imaging. While volumetric PACT imagers are commercially available, these systems typically employ a rotating gantry in which the tomographic data are sequentially acquired. Because the object varies during the data-acquisition process, the sequential data-acquisition poses challenges to image reconstruction associated with data incompleteness. The proposed method is highly significant in that it will address these challenges and enable volumetric dynamic PACT imaging with existing imagers. Aim: The aim of this study is to develop a spatiotemporal image reconstruction (STIR) method for dynamic PACT that can be applied to commercially available volumetric PACT imagers that employ a sequential scanning strategy. The proposed method aims to overcome the challenges caused by the limited number of tomographic measurements acquired per frame. Approach: A low-rank matrix estimation-based STIR method (LRME-STIR) is proposed to enable dynamic volumetric PACT. The LRME-STIR method leverages the spatiotemporal redundancies to accurately reconstruct a 4D spatiotemporal image. Results: The numerical studies substantiate the LRME-STIR method's efficacy in reconstructing 4D dynamic images from measurements acquired with a rotating gantry. The experimental study demonstrates the method's ability to faithfully recover the flow of a contrast agent at a frame rate of 0.1 s even when only a single tomographic measurement per frame is available. Conclusions: The LRME-STIR method offers a promising solution to the challenges faced by enabling 4D dynamic imaging using commercially available volumetric imagers. By enabling accurate 4D reconstruction, this method has the potential to advance preclinical research.

Ultrasound computed tomography (USCT) is an emerging imaging modality that holds great promise for breast imaging. Full-waveform inversion (FWI)-based image reconstruction methods incorporate accurate wave physics to produce high spatial resolution quantitative images of speed of sound or other acoustic properties of the breast tissues from USCT measurement data. However, the high computational cost of FWI reconstruction represents a significant burden for its widespread application in a clinical setting. The research reported here investigates the use of a convolutional neural network (CNN) to learn a mapping from USCT waveform data to speed of sound estimates. The CNN was trained using a supervised approach with a task-informed loss function aiming at preserving features of the image that are relevant to the detection of lesions. A large set of anatomically and physiologically realistic numerical breast phantoms (NBPs) and corresponding simulated USCT measurements was employed during training. Once trained, the CNN can perform real-time FWI image reconstruction from USCT waveform data. The performance of the proposed method was assessed and compared against FWI using a hold-out sample of 41 NBPs and corresponding USCT data. Accuracy was measured using relative mean square error (RMSE), structural self-similarity index measure (SSIM), and lesion detection performance (DICE score). This numerical experiment demonstrates that a supervised learning model can achieve accuracy comparable to FWI in terms of RMSE and SSIM, and better performance in terms of task performance, while significantly reducing computational time., Comment: 13 pages, 12 figures

Medical imaging systems are often evaluated and optimized via objective, or task-specific, measures of image quality (IQ) that quantify the performance of an observer on a specific clinically-relevant task. The performance of the Bayesian Ideal Observer (IO) sets an upper limit among all observers, numerical or human, and has been advocated for use as a figure-of-merit (FOM) for evaluating and optimizing medical imaging systems. However, the IO test statistic corresponds to the likelihood ratio that is intractable to compute in the majority of cases. A sampling-based method that employs Markov-Chain Monte Carlo (MCMC) techniques was previously proposed to estimate the IO performance. However, current applications of MCMC methods for IO approximation have been limited to a small number of situations where the considered distribution of to-be-imaged objects can be described by a relatively simple stochastic object model (SOM). As such, there remains an important need to extend the domain of applicability of MCMC methods to address a large variety of scenarios where IO-based assessments are needed but the associated SOMs have not been available. In this study, a novel MCMC method that employs a generative adversarial network (GAN)-based SOM, referred to as MCMC-GAN, is described and evaluated. The MCMC-GAN method was quantitatively validated by use of test-cases for which reference solutions were available. The results demonstrate that the MCMC-GAN method can extend the domain of applicability of MCMC methods for conducting IO analyses of medical imaging systems., Comment: Submitted to IEEE Transactions on Medical Imaging

We discuss solution algorithms for calibrating a tumor growth model using imaging data posed as a deterministic inverse problem. The forward model consists of a nonlinear and time-dependent reaction-diffusion partial differential equation (PDE) with unknown parameters (diffusivity and proliferation rate) being spatial fields. We use a dimension-independent globalized, inexact Newton Conjugate Gradient algorithm to solve the PDE-constrained optimization. The required gradient and Hessian actions are also presented using the adjoint method and Lagrangian formalism.

Ultrasound computed tomography (USCT) is an emerging medical imaging modality that holds great promise for improving human health. Full-waveform inversion (FWI)-based image reconstruction methods account for the relevant wave physics to produce high spatial resolution images of the acoustic properties of the breast tissues. A practical USCT design employs a circular ring-array comprised of elevation-focused ultrasonic transducers, and volumentric imaging is achieved by translating the ring-array orthogonally to the imaging plane. In commonly deployed slice-by-slice (SBS) reconstruction approaches, the three-dimensional (3D) volume is reconstructed by stacking together two-dimensional (2D) images reconstructed for each position of the ring-array. A limitation of the SBS reconstruction approach is that it does not account for 3D wave propagation physics and the focusing properties of the transducers, which can result in significant image artifacts and inaccuracies. To perform 3D image reconstruction when elevation-focused transducers are employed, a numerical description of the focusing properties of the transducers should be included in the forward model. To address this, a 3D computational model of an elevation-focused transducer is developed to enable 3D FWI-based reconstruction methods to be deployed in ring-array-based USCT. The focusing is achieved by applying a spatially varying temporal delay to the ultrasound pulse (emitter mode) and recorded signal (receiver mode). The proposed numerical transducer model is quantitatively validated and employed in computer-simulation studies that demonstrate its use in image reconstruction for ring-array USCT

Reliably predicting the future spread of brain tumors using imaging data and on a subject-specific basis requires quantifying uncertainties in data, biophysical models of tumor growth, and spatial heterogeneity of tumor and host tissue. This work introduces a Bayesian framework to calibrate the spatial distribution of the parameters within a tumor growth model to quantitative magnetic resonance imaging (MRI) data and demonstrates its implementation in a pre-clinical model of glioma. The framework leverages an atlas-based brain segmentation of grey and white matter to establish subject-specific priors and tunable spatial dependencies of the model parameters in each region. Using this framework, the tumor-specific parameters are calibrated from quantitative MRI measurements early in the course of tumor development in four rats and used to predict the spatial development of the tumor at later times. The results suggest that the tumor model, calibrated by animal-specific imaging data at one time point, can accurately predict tumor shapes with a Dice coefficient > 0.89. However, the reliability of the predicted volume and shape of tumors strongly relies on the number of earlier imaging time points used for calibrating the model. This study demonstrates, for the first time, the ability to determine the uncertainty in the inferred tissue heterogeneity and the model predicted tumor shape.

We propose derivative-informed neural operators (DINOs), a general family of neural networks to approximate operators as infinite-dimensional mappings from input function spaces to output function spaces or quantities of interest. After discretizations both inputs and outputs are high-dimensional. We aim to approximate not only the operators with improved accuracy but also their derivatives (Jacobians) with respect to the input function-valued parameter to empower derivative-based algorithms in many applications, e.g., Bayesian inverse problems, optimization under parameter uncertainty, and optimal experimental design. The major difficulties include the computational cost of generating derivative training data and the high dimensionality of the problem leading to large training cost. To address these challenges, we exploit the intrinsic low-dimensionality of the derivatives and develop algorithms for compressing derivative information and efficiently imposing it in neural operator training yielding derivative-informed neural operators. We demonstrate that these advances can significantly reduce the costs of both data generation and training for large classes of problems (e.g., nonlinear steady state parametric PDE maps), making the costs marginal or comparable to the costs without using derivatives, and in particular independent of the discretization dimension of the input and output functions. Moreover, we show that the proposed DINO achieves significantly higher accuracy than neural operators trained without derivative information, for both function approximation and derivative approximation (e.g., Gauss-Newton Hessian), especially when the training data are limited.

Dynamic imaging is essential for analyzing various biological systems and behaviors but faces two main challenges: data incompleteness and computational burden. For many imaging systems, high frame rates and short acquisition times require severe undersampling, which leads to data incompleteness. Multiple images may then be compatible with the data, thus requiring special techniques (regularization) to ensure the uniqueness of the reconstruction. Computational and memory requirements are particularly burdensome for three-dimensional dynamic imaging applications requiring high resolution in both space and time. Exploiting redundancies in the object's spatiotemporal features is key to addressing both challenges. This contribution investigates neural fields, or implicit neural representations, to model the sought-after dynamic object. Neural fields are a particular class of neural networks that represent the dynamic object as a continuous function of space and time, thus avoiding the burden of storing a full resolution image at each time frame. Neural field representation thus reduces the image reconstruction problem to estimating the network parameters via a nonlinear optimization problem (training). Once trained, the neural field can be evaluated at arbitrary locations in space and time, allowing for high-resolution rendering of the object. Key advantages of the proposed approach are that neural fields automatically learn and exploit redundancies in the sought-after object to both regularize the reconstruction and significantly reduce memory storage requirements. The feasibility of the proposed framework is illustrated with an application to dynamic image reconstruction from severely undersampled circular Radon transform data., Comment: 13 pages, 10 figures

Tomographic imaging is in general an ill-posed inverse problem. Typically, a single regularized image estimate of the sought-after object is obtained from tomographic measurements. However, there may be multiple objects that are all consistent with the same measurement data. The ability to generate such alternate solutions is important because it may enable new assessments of imaging systems. In principle, this can be achieved by means of posterior sampling methods. In recent years, deep neural networks have been employed for posterior sampling with promising results. However, such methods are not yet for use with large-scale tomographic imaging applications. On the other hand, empirical sampling methods may be computationally feasible for large-scale imaging systems and enable uncertainty quantification for practical applications. Empirical sampling involves solving a regularized inverse problem within a stochastic optimization framework to obtain alternate data-consistent solutions. In this work, a new empirical sampling method is proposed that computes multiple solutions of a tomographic inverse problem that are consistent with the same acquired measurement data. The method operates by repeatedly solving an optimization problem in the latent space of a style-based generative adversarial network (StyleGAN), and was inspired by the Photo Upsampling via Latent Space Exploration (PULSE) method that was developed for super-resolution tasks. The proposed method is demonstrated and analyzed via numerical studies that involve two stylized tomographic imaging modalities. These studies establish the ability of the method to perform efficient empirical sampling and uncertainty quantification., Comment: Submitted to IEEE Transactions on Medical Imaging

Bayesian inference provides a systematic framework for integration of data with mathematical models to quantify the uncertainty in the solution of the inverse problem. However, the solution of Bayesian inverse problems governed by complex forward models described by partial differential equations (PDEs) remains prohibitive with black-box Markov chain Monte Carlo (MCMC) methods. We present hIPPYlib-MUQ, an extensible and scalable software framework that contains implementations of state-of-the art algorithms aimed to overcome the challenges of high-dimensional, PDE-constrained Bayesian inverse problems. These algorithms accelerate MCMC sampling by exploiting the geometry and intrinsic low-dimensionality of parameter space via derivative information and low rank approximation. The software integrates two complementary open-source software packages, hIPPYlib and MUQ. hIPPYlib solves PDE-constrained inverse problems using automatically-generated adjoint-based derivatives, but it lacks full Bayesian capabilities. MUQ provides a spectrum of powerful Bayesian inversion models and algorithms, but expects forward models to come equipped with gradients and Hessians to permit large-scale solution. By combining these two libraries, we created a robust, scalable, and efficient software framework that realizes the benefits of each and allows us to tackle complex large-scale Bayesian inverse problems. To illustrate the capabilities of hIPPYlib-MUQ, we present a comparison of a number of MCMC methods on several inverse problems. These include problems with linear and nonlinear PDEs, various noise models, and different parameter dimensions. The results demonstrate that large ($\sim 50\times$) speedups over conventional black box and gradient-based MCMC algorithms can be obtained by exploiting Hessian information (from the log posterior), underscoring the power of the integrated hIPPYlib-MUQ framework.

This work develops a multiphase thermomechanical model of porous silica aerogel and implements an uncertainty analysis framework consisting of the Sobol methods for global sensitivity analyses and Bayesian inference using a set of experimental data of silica aerogel. A notable feature of this work is implementing a new noise model within the Bayesian inversion to account for data uncertainty and modeling error. The hyper-parameters in the likelihood balance data misfit and prior contribution to the parameter posteriors and prevent their biased estimation. The results indicate that the uncertainty in solid conductivity and elasticity are the most influential parameters affecting the model output variance. Also, the Bayesian inference shows that despite the microstructural randomness in the thermal measurements, the model captures the data with 2% error. However, the model is inadequate in simulating the stress-strain measurements resulting in significant uncertainty in the computational prediction of a building insulation component., Comment: Jingye Tan and Pedram Maleki contributed equally to this work

This paper presents the use of element-based algebraic multigrid (AMGe) hierarchies, implemented in the ParELAG (Parallel Element Agglomeration Algebraic Multigrid Upscaling and Solvers) library, to produce multilevel preconditioners and solvers for H(curl) and H(div) formulations. ParELAG constructs hierarchies of compatible nested spaces, forming an exact de Rham sequence on each level. This allows the application of hybrid smoothers on all levels and AMS (Auxiliary-space Maxwell Solver) or ADS (Auxiliary-space Divergence Solver) on the coarsest levels, obtaining complete multigrid cycles. Numerical results are presented, showing the parallel performance of the proposed methods. As a part of the exposition, this paper demonstrates some of the capabilities of ParELAG and outlines some of the components and procedures within the library.

Ultrasound computed tomography (USCT) is an emerging imaging modality for breast imaging that can produce quantitative images that depict the acoustic properties of tissues. Computer-simulation studies, also known as virtual imaging trials, provide researchers with an economical and convenient route to systematically explore imaging system designs and image reconstruction methods. When simulating an imaging technology intended for clinical use, it is essential to employ realistic numerical phantoms that can facilitate the objective, or task-based, assessment of image quality. Moreover, when computing objective image quality measures, an ensemble of such phantoms should be employed that display the variability in anatomy and object properties that is representative of the to-be-imaged patient cohort. Such stochastic phantoms for clinically relevant applications of USCT are currently lacking. In this work, a methodology for producing realistic three-dimensional (3D) numerical breast phantoms for enabling clinically relevant computer-simulation studies of USCT breast imaging is presented. By extending and adapting an existing stochastic 3D breast phantom for use withUSCT, methods for creating ensembles of numerical acoustic breast phantoms are established. These breast phantoms will possess clinically relevant variations in breast size, composition, acoustic properties, tumor locations, and tissue textures. To demonstrate the use of the phantoms in virtual USCT studies, two brief case studies are presented that address the development and assessment of image reconstruction procedures. Examples of breast phantoms produced by use of the proposed methods and a collection of 52 sets of simulated USCT measurement data have been made open source for use in image reconstruction development

The Alternating Direction Method of Multipliers (ADMM) provides a natural way of solving inverse problems with multiple partial differential equations (PDE) forward models and nonsmooth regularization. ADMM allows splitting these large-scale inverse problems into smaller, simpler sub-problems, for which computationally efficient solvers are available. In particular, we apply large-scale second-order optimization methods to solve the fully-decoupled Tikhonov regularized inverse problems stemming from each PDE forward model. We use fast proximal methods to handle the nonsmooth regularization term. In this work, we discuss several adaptations (such as the choice of the consensus norm) needed to maintain consistency with the underlining infinite-dimensional problem. We present two imaging applications inspired by electrical impedance tomography and quantitative photoacoustic tomography to demonstrate the proposed method's effectiveness.

Uncertainty quantification of groundwater (GW) aquifer parameters is critical for efficient management and sustainable extraction of GW resources. These uncertainties are introduced by the data, model, and prior information on the parameters. Here we develop a Bayesian inversion framework that uses Interferometric Synthetic Aperture Radar (InSAR) surface deformation data to infer the laterally heterogeneous permeability of a transient linear poroelastic model of a confined GW aquifer. The Bayesian solution of this inverse problem takes the form of a posterior probability density of the permeability. Exploring this posterior using classical Markov chain Monte Carlo (MCMC) methods is computationally prohibitive due to the large dimension of the discretized permeability field and the expense of solving the poroelastic forward problem. However, in many partial differential equation (PDE)-based Bayesian inversion problems, the data are only informative in a few directions in parameter space. For the poroelasticity problem, we prove this property theoretically for a one-dimensional problem and demonstrate it numerically for a three-dimensional aquifer model. We design a generalized preconditioned Crank--Nicolson (gpCN) MCMC method that exploits this intrinsic low dimensionality by using a low-rank based Laplace approximation of the posterior as a proposal, which we build scalably. The feasibility of our approach is demonstrated through a real GW aquifer test in Nevada. The inherently two dimensional nature of InSAR surface deformation data informs a sufficient number of modes of the permeability field to allow detection of major structures within the aquifer, significantly reducing the uncertainty in the pressure and the displacement quantities of interest.

Multiscale models of materials, consisting of upscaling discrete simulations to continuum models, are unique in their capability to simulate complex materials behavior. The fundamental limitation in multiscale models is the presence of uncertainty in the computational predictions delivered by them. In this work, a sequential multiscale model has been developed, incorporating discrete dislocation dynamics (DDD) simulations and a strain gradient plasticity (SGP) model to predict the size effect in plastic deformations of metallic micro-pillars. The DDD simulations include uniaxial compression of micro-pillars with different sizes and over a wide range of initial dislocation densities and spatial distributions of dislocations. An SGP model is employed at the continuum level that accounts for the size-dependency of flow stress and hardening rate. Sequences of uncertainty analyses have been performed to assess the predictive capability of the multiscale model. The variance-based global sensitivity analysis determines the effect of parameter uncertainty on the SGP model prediction. The multiscale model is then constructed by calibrating the continuum model using the data furnished by the DDD simulations. A Bayesian calibration method is implemented to quantify the uncertainty due to microstructural randomness in discrete dislocation simulations (density and spatial distributions of dislocations) on the macroscopic continuum model prediction (size effect in plastic deformation). The outcomes of this study indicate that the discrete-continuum multiscale model can accurately simulate the plastic deformation of micro-pillars, despite the significant uncertainty in the DDD results. Additionally, depending on the macroscopic features represented by the DDD simulations, the SGP model can reliably predict the size effect in plasticity responses of the micropillars with below 10% of error, Comment: Submitted to the International Journal of Plasticity

Many-query problems, arising from uncertainty quantification, Bayesian inversion, Bayesian optimal experimental design, and optimization under uncertainty-require numerous evaluations of a parameter-to-output map. These evaluations become prohibitive if this parametric map is high-dimensional and involves expensive solution of partial differential equations (PDEs). To tackle this challenge, we propose to construct surrogates for high-dimensional PDE-governed parametric maps in the form of projected neural networks that parsimoniously capture the geometry and intrinsic low-dimensionality of these maps. Specifically, we compute Jacobians of these PDE-based maps, and project the high-dimensional parameters onto a low-dimensional derivative-informed active subspace; we also project the possibly high-dimensional outputs onto their principal subspace. This exploits the fact that many high-dimensional PDE-governed parametric maps can be well-approximated in low-dimensional parameter and output subspace. We use the projection basis vectors in the active subspace as well as the principal output subspace to construct the weights for the first and last layers of the neural network, respectively. This frees us to train the weights in only the low-dimensional layers of the neural network. The architecture of the resulting neural network captures to first order, the low-dimensional structure and geometry of the parametric map. We demonstrate that the proposed projected neural network achieves greater generalization accuracy than a full neural network, especially in the limited training data regime afforded by expensive PDE-based parametric maps. Moreover, we show that the number of degrees of freedom of the inner layers of the projected network is independent of the parameter and output dimensions, and high accuracy can be achieved with weight dimension independent of the discretization dimension.

In this work we develop a new hierarchical multilevel approach to generate Gaussian random field realizations in an algorithmically scalable manner that is well-suited to incorporate into multilevel Markov chain Monte Carlo (MCMC) algorithms. This approach builds off of other partial differential equation (PDE) approaches for generating Gaussian random field realizations; in particular, a single field realization may be formed by solving a reaction-diffusion PDE with a spatial white noise source function as the righthand side. While these approaches have been explored to accelerate forward uncertainty quantification tasks, e.g. multilevel Monte Carlo, the previous constructions are not directly applicable to multilevel MCMC frameworks which build fine scale random fields in a hierarchical fashion from coarse scale random fields. Our new hierarchical multilevel method relies on a hierarchical decomposition of the white noise source function in $L^2$ which allows us to form Gaussian random field realizations across multiple levels of discretization in a way that fits into multilevel MCMC algorithmic frameworks. After presenting our main theoretical results and numerical scaling results to showcase the utility of this new hierarchical PDE method for generating Gaussian random field realizations, this method is tested on a four-level MCMC algorithm to explore its feasibility.

Non-smooth regularization is widely used in image reconstruction to eliminate the noise while preserving subtle image structures. In this work, we investigate the use of proximal Newton (PN) method to solve an optimization problem with a smooth data-fidelity term and total variation (TV) regularization arising from image reconstruction applications. Specifically, we consider a nonlinear Poisson-modeled single-energy X-ray computed tomography reconstruction problem with the data-fidelity term given by the I-divergence. The PN algorithm is compared to state-of-the-art first-order proximal algorithms, such as the well-established fast iterative shrinkage and thresholding algorithm (FISTA), both in terms of the number of iterations and time to solutions. We discuss the key factors that influence the performance of PN, including the strength of regularization, the stopping criterion for both sub-problem and main-problem, and the use of exact or approximated Hessian operators.

We present an extensible software framework, hIPPYlib, for solution of large-scale deterministic and Bayesian inverse problems governed by partial differential equations (PDEs) with infinite-dimensional parameter fields (which are high-dimensional after discretization). hIPPYlib overcomes the prohibitive nature of Bayesian inversion for this class of problems by implementing state-of-the-art scalable algorithms for PDE-based inverse problems that exploit the structure of the underlying operators, notably the Hessian of the log-posterior. The key property of the algorithms implemented in hIPPYlib is that the solution of the deterministic and linearized Bayesian inverse problem is computed at a cost, measured in linearized forward PDE solves, that is independent of the parameter dimension. The mean of the posterior is approximated by the MAP point, which is found by minimizing the negative log-posterior. This deterministic nonlinear least-squares optimization problem is solved with an inexact matrix-free Newton-CG method. The posterior covariance is approximated by the inverse of the Hessian of the negative log posterior evaluated at the MAP point. This Gaussian approximation is exact when the parameter-to-observable map is linear; otherwise, its logarithm agrees to two derivatives with the log-posterior at the MAP point, and thus it can serve as a proposal for Hessian-based MCMC methods. The construction of the posterior covariance is made tractable by invoking a low-rank approximation of the Hessian of the log-likelihood. Scalable tools for sample generation are also implemented. hIPPYlib makes all of these advanced algorithms easily accessible to domain scientists and provides an environment that expedites the development of new algorithms. hIPPYlib is also a teaching tool to educate researchers and practitioners who are new to inverse problems and the Bayesian inference framework.

This paper develops a multifidelity method that enables estimation of failure probabilities for expensive-to-evaluate models via information fusion and importance sampling. The presented general fusion method combines multiple probability estimators with the goal of variance reduction. We use low-fidelity models to derive biasing densities for importance sampling and then fuse the importance sampling estimators such that the fused multifidelity estimator is unbiased and has mean-squared error lower than or equal to that of any of the importance sampling estimators alone. By fusing all available estimators, the method circumvents the challenging problem of selecting the best biasing density and using only that density for sampling. A rigorous analysis shows that the fused estimator is optimal in the sense that it has minimal variance amongst all possible combinations of the estimators. The asymptotic behavior of the proposed method is demonstrated on a convection-diffusion-reaction partial differential equation model for which $10^5$ samples can be afforded. To illustrate the proposed method at scale, we consider a model of a free plane jet and quantify how uncertainties at the flow inlet propagate to a quantity of interest related to turbulent mixing. Compared to an importance sampling estimator that uses the high-fidelity model alone, our multifidelity estimator reduces the required CPU time by 65\% while achieving a similar coefficient of variation.

In this work we develop a scalable computational framework for the solution of PDE-constrained optimal control under high-dimensional uncertainty. Specifically, we consider a mean-variance formulation of the control objective and employ a Taylor expansion with respect to the uncertain parameter either to directly approximate the control objective or as a control variate for variance reduction. The expressions for the mean and variance of the Taylor approximation are known analytically, although their evaluation requires efficient computation of the trace of the (preconditioned) Hessian of the control objective. We propose to estimate this trace by solving a generalized eigenvalue problem using a randomized algorithm that only requires the action of the Hessian on a small number of random directions. Then, the computational work does not depend on the nominal dimension of the uncertain parameter, but depends only on the effective dimension, thus ensuring scalability to high-dimensional problems. Moreover, to increase the estimation accuracy of the mean and variance of the control objective by the Taylor approximation, we use it as a control variate for variance reduction, which results in considerable computational savings (several orders of magnitude) compared to a plain Monte Carlo method. We demonstrate the accuracy, efficiency, and scalability of the proposed computational method for two examples with high-dimensional uncertain parameters: subsurface flow in a porous medium modeled as an elliptic PDE, and turbulent jet flow modeled by the Reynolds-averaged Navier--Stokes equations coupled with a nonlinear advection-diffusion equation characterizing model uncertainty. In particular, for the latter more challenging example we show scalability of our algorithm up to one million parameters resulting from discretization of the uncertain parameter field.

This work describes a domain embedding technique between two non-matching meshes used for generating realizations of spatially correlated random fields with applications to large-scale sampling-based uncertainty quantification. The goal is to apply the multilevel Monte Carlo (MLMC) method for the quantification of output uncertainties of PDEs with random input coefficients on general, unstructured computational domains. We propose a highly scalable, hierarchical sampling method to generate realizations of a Gaussian random field on a given unstructured mesh by solving a reaction-diffusion PDE with a stochastic right-hand side. The stochastic PDE is discretized using the mixed finite element method on an embedded domain with a structured mesh, and then the solution is projected onto the unstructured mesh. This work describes implementation details on how to efficiently transfer data from the structured and unstructured meshes at coarse levels, assuming this can be done efficiently on the finest level. We investigate the efficiency and parallel scalability of the technique for the scalable generation of Gaussian random fields in three dimensions. An application of the MLMC method is presented for quantifying uncertainties of subsurface flow problems. We demonstrate the scalability of the sampling method with non-matching mesh embedding, coupled with a parallel forward model problem solver, for large-scale 3D MLMC simulations with up to $1.9\cdot 10^9$ unknowns., Comment: To appear in Numerical Linear Algebra with Applications

Numerical solution of the system of partial differential equations arising from the Reynolds-Averaged Navier-Stokes (RANS) equations with $k-\epsilon$ turbulence model presents several challenges due to the advection dominated nature of the problem and the presence of highly nonlinear reaction terms. State-of-the-art software for the numerical solution of the RANS equations address these challenges by introducing non-differentiable perturbations in the model to ensure numerical stability. However, this approach leads to difficulties in the formulation of the higher-order forward/adjoint problems, which are needed for scalable Hessian-based uncertainty quantification (UQ) methods. In this note, we present the construction of a UQ-ready flow solver, i.e., one that is smoothly differentiable and provides not only forward solver capabilities but also adjoint and higher derivatives capabilities.

In this work we propose and analyze a Hessian-based adaptive sparse quadrature to compute infinite-dimensional integrals with respect to the posterior distribution in the context of Bayesian inverse problems with Gaussian prior. Due to the concentration of the posterior distribution in the domain of the prior distribution, a prior-based parametrization and sparse quadrature may fail to capture the posterior distribution and lead to erroneous evaluation results. By using a parametrization based on the Hessian of the negative log-posterior, the adaptive sparse quadrature can effectively allocate the quadrature points according to the posterior distribution. A dimension-independent convergence rate of the proposed method is established under certain assumptions on the Gaussian prior and the integrands. Dimension-independent and faster convergence than $O(N^{-1/2})$ is demonstrated for a linear as well as a nonlinear inverse problem whose posterior distribution can be effectively approximated by a Gaussian distribution at the MAP point.

We propose an alternative method to generate samples of a spatially correlated random field with applications to large-scale problems for forward propagation of uncertainty. A classical approach for generating these samples is the Karhunen-Lo\`{e}ve (KL) decomposition. However, the KL expansion requires solving a dense eigenvalue problem and is therefore computationally infeasible for large-scale problems. Sampling methods based on stochastic partial differential equations provide a highly scalable way to sample Gaussian fields, but the resulting parametrization is mesh dependent. We propose a multilevel decomposition of the stochastic field to allow for scalable, hierarchical sampling based on solving a mixed finite element formulation of a stochastic reaction-diffusion equation with a random, white noise source function. Numerical experiments are presented to demonstrate the scalability of the sampling method as well as numerical results of multilevel Monte Carlo simulations for a subsurface porous media flow application using the proposed sampling method., Comment: 22 pages, 9 figures, To appear in SIAM Journal on Scientific Computing Special Session: 2016 Copper Mountain Conference

Current state of the art preconditioners for the reduced Hessian and the Karush-Kuhn-Tucker (KKT) operator for large scale inverse problems are typically based on approximating the reduced Hessian with the regularization operator. However, the quality of this approximation degrades with increasingly informative observations or data. Thus the best case scenario from a scientific standpoint (fully informative data) is the worse case scenario from a computational perspective. In this paper we present an augmented Lagrangian-type preconditioner based on a block diagonal approximation of the augmented upper left block of the KKT operator. The preconditioner requires solvers for two linear subproblems that arise in the augmented KKT operator, which we expect to be much easier to precondition than the reduced Hessian. Analysis of the spectrum of the preconditioned KKT operator indicates that the preconditioner is effective when the regularization is chosen appropriately. In particular, it is effective when the regularization does not over-penalize highly informed parameter modes and does not under-penalize uninformed modes. Finally, we present a numerical study for a large data/low noise Poisson source inversion problem, demonstrating the effectiveness of the preconditioner. In this example, three MINRES iterations on the KKT system with our preconditioner results in a reconstruction with better accuracy than 50 iterations of CG on the reduced Hessian system with regularization preconditioning., Comment: 30 pages, 4 figures, 1 table. Accepted for publication in SIAM Journal on Scientific Computing (SISC)