Search

Equivariant Imaging (EI) regularization has become the de-facto technique for unsupervised training of deep imaging networks, without any need of ground-truth data. Observing that the EI-based unsupervised training paradigm currently has significant computational redundancy leading to inefficiency in high-dimensional applications, we propose a sketched EI regularization which leverages the randomized sketching techniques for acceleration. We then extend our sketched EI regularization to develop an accelerated deep internal learning framework -- Sketched Equivariant Deep Image Prior (Sk-EI-DIP), which can be efficiently applied for single-image and task-adapted reconstruction. Additionally, for network adaptation tasks, we propose a parameter-efficient approach for accelerating both EI-DIP and Sk-EI-DIP via optimizing only the normalization layers. Our numerical study on X-ray CT image reconstruction tasks demonstrate that our approach can achieve order-of-magnitude computational acceleration over standard EI-based counterpart in single-input setting, and network adaptation at test time.

Identifying dynamical system (DS) is a vital task in science and engineering. Traditional methods require numerous calls to the DS solver, rendering likelihood-based or least-squares inference frameworks impractical. For efficient parameter inference, two state-of-the-art techniques are the kernel method for modeling and the "one-step framework" for jointly inferring unknown parameters and hyperparameters. The kernel method is a quick and straightforward technique, but it cannot estimate solutions and their derivatives, which must strictly adhere to physical laws. We propose a model-embedded "one-step" Bayesian framework for joint inference of unknown parameters and hyperparameters by maximizing the marginal likelihood. This approach models the solution and its derivatives using Gaussian process regression (GPR), taking into account smoothness and continuity properties, and treats differential equations as constraints that can be naturally integrated into the Bayesian framework in the linear case. Additionally, we prove the convergence of the model-embedded Gaussian process regression (ME-GPR) for theoretical development. Motivated by Taylor expansion, we introduce a piecewise first-order linearization strategy to handle nonlinear dynamic systems. We derive estimates and confidence intervals, demonstrating that they exhibit low bias and good coverage properties for both simulated models and real data., Comment: 24 pages, 3 figures, 5 tables

Bayesian Experimental Design (BED), which aims to find the optimal experimental conditions for Bayesian inference, is usually posed as to optimize the expected information gain (EIG). The gradient information is often needed for efficient EIG optimization, and as a result the ability to estimate the gradient of EIG is essential for BED problems. The primary goal of this work is to develop methods for estimating the gradient of EIG, which, combined with the stochastic gradient descent algorithms, result in efficient optimization of EIG. Specifically, we first introduce a posterior expected representation of the EIG gradient with respect to the design variables. Based on this, we propose two methods for estimating the EIG gradient, UEEG-MCMC that leverages posterior samples generated through Markov Chain Monte Carlo (MCMC) to estimate the EIG gradient, and BEEG-AP that focuses on achieving high simulation efficiency by repeatedly using parameter samples. Theoretical analysis and numerical studies illustrate that UEEG-MCMC is robust agains the actual EIG value, while BEEG-AP is more efficient when the EIG value to be optimized is small. Moreover, both methods show superior performance compared to several popular benchmarks in our numerical experiments.

Combining the strengths of model-based iterative algorithms and data-driven deep learning solutions, deep unrolling networks (DuNets) have become a popular tool to solve inverse imaging problems. While DuNets have been successfully applied to many linear inverse problems, nonlinear problems tend to impair the performance of the method. Inspired by momentum acceleration techniques that are often used in optimization algorithms, we propose a recurrent momentum acceleration (RMA) framework that uses a long short-term memory recurrent neural network (LSTM-RNN) to simulate the momentum acceleration process. The RMA module leverages the ability of the LSTM-RNN to learn and retain knowledge from the previous gradients. We apply RMA to two popular DuNets -- the learned proximal gradient descent (LPGD) and the learned primal-dual (LPD) methods, resulting in LPGD-RMA and LPD-RMA respectively. We provide experimental results on two nonlinear inverse problems: a nonlinear deconvolution problem, and an electrical impedance tomography problem with limited boundary measurements. In the first experiment we have observed that the improvement due to RMA largely increases with respect to the nonlinearity of the problem. The results of the second example further demonstrate that the RMA schemes can significantly improve the performance of DuNets in strongly ill-posed problems.

Entropy estimation is of practical importance in information theory and statistical science. Many existing entropy estimators suffer from fast growing estimation bias with respect to dimensionality, rendering them unsuitable for high-dimensional problems. In this work we propose a transform-based method for high-dimensional entropy estimation, which consists of the following two main ingredients. First by modifying the k-NN based entropy estimator, we propose a new estimator which enjoys small estimation bias for samples that are close to a uniform distribution. Second we design a normalizing flow based mapping that pushes samples toward a uniform distribution, and the relation between the entropy of the original samples and the transformed ones is also derived. As a result the entropy of a given set of samples is estimated by first transforming them toward a uniform distribution and then applying the proposed estimator to the transformed samples. The performance of the proposed method is compared against several existing entropy estimators, with both mathematical examples and real-world applications.

Bayesian methods for solving inverse problems are a powerful alternative to classical methods since the Bayesian approach offers the ability to quantify the uncertainty in the solution. In recent years, data-driven techniques for solving inverse problems have also been remarkably successful, due to their superior representation ability. In this work, we incorporate data-based models into a class of Langevin-based sampling algorithms for Bayesian inference in imaging inverse problems. In particular, we introduce NF-ULA (Normalizing Flow-based Unadjusted Langevin algorithm), which involves learning a normalizing flow (NF) as the image prior. We use NF to learn the prior because a tractable closed-form expression for the log prior enables the differentiation of it using autograd libraries. Our algorithm only requires a normalizing flow-based generative network, which can be pre-trained independently of the considered inverse problem and the forward operator. We perform theoretical analysis by investigating the well-posedness and non-asymptotic convergence of the resulting NF-ULA algorithm. The efficacy of the proposed NF-ULA algorithm is demonstrated in various image restoration problems such as image deblurring, image inpainting, and limited-angle X-ray computed tomography (CT) reconstruction. NF-ULA is found to perform better than competing methods for severely ill-posed inverse problems.

The Langevin algorithms are frequently used to sample the posterior distributions in Bayesian inference. In many practical problems, however, the posterior distributions often consist of non-differentiable components, posing challenges for the standard Langevin algorithms, as they require to evaluate the gradient of the energy function in each iteration. To this end, a popular remedy is to utilize the proximity operator, and as a result one needs to solve a proximity subproblem in each iteration. The conventional practice is to solve the subproblems accurately, which can be exceedingly expensive, as the subproblem needs to be solved in each iteration. We propose an approximate primal-dual fixed-point algorithm for solving the subproblem, which only seeks an approximate solution of the subproblem and therefore reduces the computational cost considerably. We provide theoretical analysis of the proposed method and also demonstrate its performance with numerical examples.

In this work we study systems consisting of a group of moving particles. In such systems, often some important parameters are unknown and have to be estimated from observed data. Such parameter estimation problems can often be solved via a Bayesian inference framework. However in many practical problems, only data at the aggregate level is available and as a result the likelihood function is not available, which poses challenge for Bayesian methods. In particular, we consider the situation where the distributions of the particles are observed. We propose a Wasserstein distance based sequential Monte Carlo sampler to solve the problem: the Wasserstein distance is used to measure the similarity between the observed and the simulated particle distributions and the sequential Monte Carlo samplers is used to deal with the sequentially available observations. Two real-world examples are provided to demonstrate the performance of the proposed method.

Solving high-dimensional Bayesian inverse problems (BIPs) with the variational inference (VI) method is promising but still challenging. The main difficulties arise from two aspects. First, VI methods approximate the posterior distribution using a simple and analytic variational distribution, which makes it difficult to estimate complex spatially-varying parameters in practice. Second, VI methods typically rely on gradient-based optimization, which can be computationally expensive or intractable when applied to BIPs involving partial differential equations (PDEs). To address these challenges, we propose a novel approximation method for estimating the high-dimensional posterior distribution. This approach leverages a deep generative model to learn a prior model capable of generating spatially-varying parameters. This enables posterior approximation over the latent variable instead of the complex parameters, thus improving estimation accuracy. Moreover, to accelerate gradient computation, we employ a differentiable physics-constrained surrogate model to replace the adjoint method. The proposed method can be fully implemented in an automatic differentiation manner. Numerical examples demonstrate two types of log-permeability estimation for flow in heterogeneous media. The results show the validity, accuracy, and high efficiency of the proposed method.

Failure probability estimation problem is an crucial task in engineering. In this work we consider this problem in the situation that the underlying computer models are extremely expensive, which often arises in the practice, and in this setting, reducing the calls of computer model is of essential importance. We formulate the problem of estimating the failure probability with expensive computer models as an sequential experimental design for the limit state (i.e., the failure boundary) and propose a series of efficient adaptive design criteria to solve the design of experiment (DOE). In particular, the proposed method employs the deep neural network (DNN) as the surrogate of limit state function for efficiently reducing the calls of expensive computer experiment. A map from the Gaussian distribution to the posterior approximation of the limit state is learned by the normalizing flows for the ease of experimental design. Three normalizing-flows-based design criteria are proposed in this work for deciding the design locations based on the different assumption of generalization error. The accuracy and performance of the proposed method is demonstrated by both theory and practical examples., Comment: failure probability, normalizing flows, adaptive design of experiment. arXiv admin note: text overlap with arXiv:1509.04613

In many Bayesian inverse problems the goal is to recover a spatially varying random field. Such problems are often computationally challenging especially when the forward model is governed by complex partial differential equations (PDEs). The challenge is particularly severe when the spatial domain is large and the unknown random field needs to be represented by a high-dimensional parameter. In this paper, we present a domain-decomposed method to attack the dimensionality issue and the method decomposes the spatial domain and the parameter domain simultaneously. On each subdomain, a local Karhunen-Lo`eve (KL) expansion is constructed, and a local inversion problem is solved independently in a parallel manner, and more importantly, in a lower-dimensional space. After local posterior samples are generated through conducting Markov chain Monte Carlo (MCMC) simulations on subdomains, a novel projection procedure is developed to effectively reconstruct the global field. In addition, the domain decomposition interface conditions are dealt with an adaptive Gaussian process-based fitting strategy. Numerical examples are provided to demonstrate the performance of the proposed method.

Predicting the behaviors of pedestrian crowds is of critical importance for a variety of real-world problems. Data driven modeling, which aims to learn the mathematical models from observed data, is a promising tool to construct models that can make accurate predictions of such systems. In this work, we present a data-driven modeling approach based on the ODE-Net framework, for constructing continuous-time models of crowd dynamics. We discuss some challenging issues in applying the ODE-Net method to such problems, which are primarily associated with the dimensionality of the underlying crowd system, and we propose to address these issues by incorporating the social-force concept in the ODE-Net framework. Finally application examples are provided to demonstrate the performance of the proposed method.

In many real-world engineering systems, the performance or reliability of the system is characterised by a scalar parameter. The distribution of this performance parameter is important in many uncertainty quantification problems, ranging from risk management to utility optimisation. In practice, this distribution usually cannot be derived analytically and has to be obtained numerically by simulations. To this end, standard Monte Carlo simulations are often used, however, they cannot efficiently reconstruct the tail of the distribution which is essential in many applications. One possible remedy is to use the Multicanonical Monte Carlo method, an adaptive importance sampling scheme. In this method, one draws samples from an importance sampling distribution in a nonstandard form in each iteration, which is usually done via Markov chain Monte Carlo (MCMC). MCMC is inherently serial and therefore struggles with parallelism. In this paper, we present a new approach, which uses the Sequential Monte Carlo sampler to draw from the importance sampling distribution, which is particularly suited for parallel implementation. With both mathematical and practical examples, we demonstrate the competitive performance of the proposed method., Comment: 21 pages, 3 figures

In this paper, we evaluate the performance of the multilevel Monte Carlo method (MLMC) for deterministic and uncertain hyperbolic systems, where randomness is introduced either in the modeling parameters or in the approximation algorithms. MLMC is a well known variance reduction method widely used to accelerate Monte Carlo (MC) sampling. However, we demonstrate in this paper that for hyperbolic systems, whether MLMC can achieve a real boost turns out to be delicate. The computational costs of MLMC and MC depend on the interplay among the accuracy (bias) and the computational cost of the numerical method for a single sample, as well as the variances of the sampled MLMC corrections or MC solutions. We characterize three regimes for the MLMC and MC performances using those parameters, and show that MLMC may not accelerate MC and can even have a higher cost when the variances of MC solutions and MLMC corrections are of the same order. Our studies are carried out by a few prototype hyperbolic systems: a linear scalar equation, the Euler and shallow water equations, and a linear relaxation model, the above statements are proved analytically in some cases, and demonstrated numerically for the cases of the stochastic hyperbolic equations driven by white noise parameters and Glimm's random choice method for deterministic hyperbolic equations.

Elucidating emergent regularities in intriguing crowd dynamics is a fundamental scientific problem arising in multiple fields. In this work, based on the social force model, we simulate the typical scenario of collective escape towards a single exit and reveal the striking analogy of crowd dynamics and crystallisation. With the outflow of the pedestrians, crystalline order emerges in the compact crowd. In this process, the local misalignment and global rearrangement of pedestrians are well rationalized in terms of the characteristic motions of topological defects in the crystal. Exploiting the notions from the physics of crystallisation further reveals the emergence of multiple fast tracks in the collective escape., Comment: 8 pages, 6 figures

We propose an affine-mapping based variational Ensemble Kalman filter for sequential Bayesian filtering problems with generic observation models. Specifically, the proposed method is formulated as to construct an affine mapping from the prior ensemble to the posterior one, and the affine mapping is computed via a variational Bayesian formulation, i.e., by minimizing the Kullback-Leibler divergence between the transformed distribution through the affine mapping and the actual posterior. Some theoretical properties of resulting optimization problem are studied and a gradient descent scheme is proposed to solve the resulting optimization problem. With numerical examples we demonstrate that the method has competitive performance against existing methods.

In this work we consider Bayesian inference problems with intractable likelihood functions. We present a method to compute an approximate of the posterior with a limited number of model simulations. The method features an inverse Gaussian Process regression (IGPR), i.e., one from the output of a simulation model to the input of it. Within the method, we provide an adaptive algorithm with a tempering procedure to construct the approximations of the marginal posterior distributions. With examples we demonstrate that IGPR has a competitive performance compared to some commonly used algorithms, especially in terms of statistical stability and computational efficiency, while the price to pay is that it can only compute a weighted Gaussian approximation of the marginal posteriors.

Quantifying uncertainties in physical or engineering systems often requires a large number of simulations of the underlying computer models that are computationally intensive. Emulators or surrogate models are often used to accelerate the computation in such problems, and in this regard the Gaussian Process (GP) emulator is a popular choice for its ability to quantify the approximation error in the emulator itself. However, a major limitation of the GP emulator is that it can not handle problems of very high dimensions, which is often addressed with dimension reduction techniques. In this work we hope to address an issue that the models of interest are so complex that they admit different low dimensional structures in different parameter regimes. Building upon the active subspace method for dimension reduction, we propose a clustered active subspace method which identifies the local low-dimensional structures as well as the parameter regimes they are in (represented as clusters), and then construct low dimensional and local GP emulators within the clusters. Specifically we design a clustering method based on the gradient information to identify these clusters, and a local GP construction procedure to construct the GP emulator within a local cluster. With numerical examples, we demonstrate that the proposed method is effective when the underlying models are of complex low-dimensional structures.

Many real-world problems require one to estimate parameters of interest, in a Bayesian framework, from data that are collected sequentially in time. Conventional methods for sampling from posterior distributions, such as {Markov Chain Monte Carlo} can not efficiently address such problems as they do not take advantage of the data's sequential structure. To this end, sequential methods which seek to update the posterior distribution whenever a new collection of data become available are often used to solve these types of problems. Two popular choices of sequential method are the Ensemble Kalman filter (EnKF) and the sequential Monte Carlo sampler (SMCS). While EnKF only computes a Gaussian approximation of the posterior distribution, SMCS can draw samples directly from the posterior. Its performance, however, depends critically upon the kernels that are used. In this work, we present a method that constructs the kernels of SMCS using an EnKF formulation, and we demonstrate the performance of the method with numerical examples.

By facilitating the generation of samples from arbitrary probability distributions, Markov Chain Monte Carlo (MCMC) is, arguably, \emph{the} tool for the evaluation of Bayesian inference problems that yield non-standard posterior distributions. In recent years, however, it has become apparent that Sequential Monte Carlo (SMC) samplers have the potential to outperform MCMC in a number of ways. SMC samplers are better suited to highly parallel computing architectures and also feature various tuning parameters that are not available to MCMC. One such parameter - the `L-kernel' - is a user-defined probability distribution that can be used to influence the efficiency of the sampler. In the current paper, the authors explain how to derive an expression for the L-kernel that minimises the variance of the estimates realised by an SMC sampler. Various approximation methods are then proposed to aid implementation of the proposed L-kernel. The improved performance of the resulting algorithm is demonstrated in multiple scenarios. For the examples shown in the current paper, the use of an approximately optimum L-kernel has reduced the variance of the SMC estimates by up to 99 % while also reducing the number of times that resampling was required by between 65 % and 70 %. Python code and code tests accompanying this manuscript are available through the Github repository \url{https://github.com/plgreenLIRU/SMC_approx_optL}., Comment: 29 pages, 14 figures

Data collection is a critical step in statistical inference and data science, and the goal of statistical experimental design (ED) is to find the data collection setup that can provide most information for the inference. In this work we consider a special type of ED problems where the likelihoods are not available in a closed form. In this case, the popular information-theoretic Kullback-Leibler divergence (KLD) based design criterion can not be used directly, as it requires to evaluate the likelihood function. To address the issue, we derive a new utility function, which is a lower bound of the original KLD utility. This lower bound is expressed in terms of the summation of two or more entropies in the data space, and thus can be evaluated efficiently via entropy estimation methods. We provide several numerical examples to demonstrate the performance of the proposed method., Comment: To appear in AISTATS 2020

Global optimization finds applications in a wide range of real world problems. The multi-start methods are a popular class of global optimization techniques, which are based on the ideas of conducting local searches at multiple starting points. In this work we propose a new multi-start algorithm where the starting points are determined in a Bayesian optimization framework. Specifically, the method can be understood as to construct a new function by conducting local searches of the original objective function, where the new function attains the same global optima as the original one. Bayesian optimization is then applied to find the global optima of the new local search defined function.

The performance of Hamiltonian Monte Carlo (HMC) sampler depends critically on some algorithm parameters such as the total integration time and the numerical integration stepsize. The parameter tuning is particularly challenging when the mass matrix of the HMC sampler is adapted. We propose in this work a Kolmogorov-Sinai entropy (KSE) based design criterion to optimize these algorithm parameters, which can avoid some potential issues in the often used jumping-distance based measures. For near-Gaussian distributions, we are able to derive the optimal algorithm parameters with respect to the KSE criterion analytically. As a byproduct the KSE criterion also provides a theoretical justification for the need to adapt the mass matrix in HMC sampler. Based on the results, we propose an adaptive HMC algorithm, and we then demonstrate the performance of the proposed algorithm with numerical examples.

Failure probability evaluation for complex physical and engineering systems governed by partial differential equations (PDEs) are computationally intensive, especially when high-dimensional random parameters are involved. Since standard numerical schemes for solving these complex PDEs are expensive, traditional Monte Carlo methods which require repeatedly solving PDEs are infeasible. Alternative approaches which are typically the surrogate based methods suffer from the so-called ``curse of dimensionality'', which limits their application to problems with high-dimensional parameters. For this purpose, we develop a novel hierarchical neural hybrid (HNH) method to efficiently compute failure probabilities of these challenging high-dimensional problems. Especially, multifidelity surrogates are constructed based on neural networks with different levels of layers, such that expensive highfidelity surrogates are adapted only when the parameters are in the suspicious domain. The efficiency of our new HNH method is theoretically analyzed and is demonstrated with numerical experiments. From numerical results, we show that to achieve an accuracy in estimating the rare failure probability (e.g., $10^{-5}$), the traditional Monte Carlo method needs to solve PDEs more than a million times, while our HNH only requires solving them a few thousand times.

We consider supervised dimension reduction problems, namely to identify a low dimensional projection of the predictors $\-x$ which can retain the statistical relationship between $\-x$ and the response variable $y$. We follow the idea of the sliced inverse regression (SIR) and the sliced average variance estimation (SAVE) type of methods, which is to use the statistical information of the conditional distribution $\pi(\-x|y)$ to identify the dimension reduction (DR) space. In particular we focus on the task of computing this conditional distribution without slicing the data. We propose a Bayesian framework to compute the conditional distribution where the likelihood function is obtained using the Gaussian process regression model. The conditional distribution $\pi(\-x|y)$ can then be computed directly via Monte Carlo sampling. We then can perform DR by considering certain moment functions (e.g. the first or the second moment) of the samples of the posterior distribution. With numerical examples, we demonstrate that the proposed method is especially effective for small data problems.

We provide a complete framework for performing infinite-dimensional Bayesian inference and uncertainty quantification for image reconstruction with Poisson data. In particular, we address the following issues to make the Bayesian framework applicable in practice. We first introduce a positivity-preserving reparametrization, and we prove that under the reparametrization and a hybrid prior, the posterior distribution is well-posed in the infinite dimensional setting. Second we provide a dimension-independent MCMC algorithm, based on the preconditioned Crank-Nicolson Langevin method, in which we use a primal-dual scheme to compute the offset direction. Third we give a method combining the model discrepancy method and maximum likelihood estimation to determine the regularization parameter in the hybrid prior. Finally we propose to use the obtained posterior distribution to detect artifacts in a recovered image. We provide an example to demonstrate the effectiveness of the proposed method.

Bayesian inference methods have been widely applied in inverse problems, {largely due to their ability to characterize the uncertainty associated with the estimation results.} {In the Bayesian framework} the prior distribution of the unknown plays an essential role in the Bayesian inference, {and a good prior distribution can significantly improve the inference results.} In this paper, we extend the total~variation-Gaussian (TG) prior in \cite{Z.Yao2016}, and propose a hybrid prior distribution which combines the nonlocal total variation regularization and the Gaussian (NLTG) distribution. The advantage of the new prior is two-fold. The proposed prior models both texture and geometric structures present in images through the NLTV. The Gaussian reference measure also provides a flexibility of incorporating structure information from a reference image. Some theoretical properties are established for the NLTG prior. The proposed prior is applied to limited-angle tomography reconstruction problem with difficulties of severe data missing. We compute both MAP and CM estimates through two efficient methods and the numerical experiments validate the advantages and feasibility of the proposed NLTG prior.

Importance sampling (IS) is an important technique to reduce the estimation variance in Monte Carlo simulations. In many practical problems, however, the use of IS method may result in unbounded variance, and thus fail to provide reliable estimates. To address the issue, we propose a method which can prevent the risk of unbounded variance; the proposed method performs the standard IS for the integral of interest in a region only in which the IS weight is bounded and use the result as an approximation to the original integral. It can be verified that the resulting estimator has a finite variance. Moreover, we also provide a normality test based method to identify the region with bounded IS weight (termed as the safe region) from the samples drawn from the standard IS distribution. With numerical examples, we demonstrate that the proposed method can yield rather reliable estimate when the standard IS fails, and it also outperforms the defensive IS, a popular method to prevent unbounded variance.

In Bayesian inverse problems sampling the posterior distribution is often a challenging task when the underlying models are computationally intensive. To this end, surrogates or reduced models are often used to accelerate the computation. However, in many practical problems, the parameter of interest can be of high dimensionality, which renders standard model reduction techniques infeasible. In this paper, we present an approach that employs the ANOVA decomposition method to reduce the model with respect to the unknown parameters, and the reduced basis method to reduce the model with respect to the physical parameters. Moreover, we provide an adaptive scheme within the MCMC iterations, to perform the ANOVA decomposition with respect to the posterior distribution. With numerical examples, we demonstrate that the proposed model reduction method can significantly reduce the computational cost of Bayesian inverse problems, without sacrificing much accuracy.

Particle filtering (PF) is an often used method to estimate the states of dynamical systems. A major limitation of the standard PF method is that the dimensionality of the state space increases as the time proceeds and eventually may cause degeneracy of the algorithm. A possible approach to alleviate the degeneracy issue is to compute the marginal posterior distribution at each time step, which leads to the so-called marginal PF method. A key issue in the marginal PF method is to construct a good sampling distribution in the marginal space. When the posterior distribution is close to Gaussian, the Ensemble Kalman filter (EnKF) method can usually provide a good sampling distribution; however the EnKF approximation may fail completely when the posterior is strongly non-Gaussian. In this work we propose a defensive marginal PF (DMPF) algorithm which constructs a sampling distribution in the marginal space by combining the standard PF and the EnKF approximation using a multiple importance sampling (MIS) scheme. An important feature of the proposed algorithm is that it can automatically adjust the relative weight of the PF and the EnKF components in the MIS scheme in each step, according to how non-Gaussian the posterior is. With numerical examples we demonstrate that the proposed method can perform well regardless of whether the posteriors can be well approximated by Gaussian.

We study Bayesian inference methods for solving linear inverse problems, focusing on hierarchical formulations where the prior or the likelihood function depend on unspecified hyperparameters. In practice, these hyperparameters are often determined via an empirical Bayesian method that maximizes the marginal likelihood function, i.e., the probability density of the data conditional on the hyperparameters. Evaluating the marginal likelihood, however, is computationally challenging for large-scale problems. In this work, we present a method to approximately evaluate marginal likelihood functions, based on a low-rank approximation of the update from the prior covariance to the posterior covariance. We show that this approximation is optimal in a minimax sense. Moreover, we provide an efficient algorithm to implement the proposed method, based on a combination of the randomized SVD and a spectral approximation method to compute square roots of the prior covariance matrix. Several numerical examples demonstrate good performance of the proposed method.

We consider Bayesian inference problems with computationally intensive likelihood functions. We propose a Gaussian process (GP) based method to approximate the joint distribution of the unknown parameters and the data. In particular, we write the joint density approximately as a product of an approximate posterior density and an exponentiated GP surrogate. We then provide an adaptive algorithm to construct such an approximation, where an active learning method is used to choose the design points. With numerical examples, we illustrate that the proposed method has competitive performance against existing approaches for Bayesian computation.

Estimating failure probabilities of engineering systems is an important problem in many engineering fields. In this work we consider such problems where the failure probability is extremely small (e.g $\leq10^{-10}$). In this case, standard Monte Carlo methods are not feasible due to the extraordinarily large number of samples required. To address these problems, we propose an algorithm that combines the main ideas of two very powerful failure probability estimation approaches: the subset simulation (SS) and the multicanonical Monte Carlo (MMC) methods. Unlike the standard MMC which samples in the entire domain of the input parameter in each iteration, the proposed subset MMC algorithm adaptively performs MMC simulations in a subset of the state space and thus improves the sampling efficiency. With numerical examples we demonstrate that the proposed method is significantly more efficient than both of the SS and the MMC methods. Moreover, the proposed algorithm can reconstruct the complete distribution function of the parameter of interest and thus can provide more information than just the failure probabilities of the systems.

The preconditioned Crank-Nicolson (pCN) method is a Markov Chain Monte Carlo (MCMC) scheme, specifically designed to perform Bayesian inferences in function spaces. Unlike many standard MCMC algorithms, the pCN method can preserve the sampling efficiency under the mesh refinement, a property referred to as being dimension independent. In this work we consider an adaptive strategy to further improve the efficiency of pCN. In particular we develop a hybrid adaptive MCMC method: the algorithm performs an adaptive Metropolis scheme in a chosen finite dimensional subspace, and a standard pCN algorithm in the complement space of the chosen subspace. We show that the proposed algorithm satisfies certain important ergodicity conditions. Finally with numerical examples we demonstrate that the proposed method has competitive performance with existing adaptive algorithms., Comment: arXiv admin note: text overlap with arXiv:1511.05838

Many scientific and engineering problems require to perform Bayesian inferences for unknowns of infinite dimension. In such problems, many standard Markov Chain Monte Carlo (MCMC) algorithms become arbitrary slow under the mesh refinement, which is referred to as being dimension dependent. To this end, a family of dimensional independent MCMC algorithms, known as the preconditioned Crank-Nicolson (pCN) methods, were proposed to sample the infinite dimensional parameters. In this work we develop an adaptive version of the pCN algorithm, where the covariance operator of the proposal distribution is adjusted based on sampling history to improve the simulation efficiency. We show that the proposed algorithm satisfies an important ergodicity condition under some mild assumptions. Finally we provide numerical examples to demonstrate the performance of the proposed method.

In this note, we present a derivative-free trust-region (TR) algorithm for reliability based optimization (RBO) problems. The proposed algorithm consists of solving a set of subproblems, in which simple surrogate models of the reliability constraints are constructed and used in solving the subproblems. Taking advantage of the special structure of the RBO problems, we employ a sample reweighting method to evaluate the failure probabilities, which constructs the surrogate for the reliability constraints by performing only a single full reliability evaluation in each iteration. With numerical experiments, we illustrate that the proposed algorithm is competitive against existing methods.

Many scientific and engineering problems require to perform Bayesian inferences in function spaces, in which the unknowns are of infinite dimension. In such problems, choosing an appropriate prior distribution is an important task. In particular we consider problems where the function to infer is subject to sharp jumps which render the commonly used Gaussian measures unsuitable. On the other hand, the so-called total variation (TV) prior can only be defined in a finite dimensional setting, and does not lead to a well-defined posterior measure in function spaces. In this work we present a TV-Gaussian (TG) prior to address such problems, where the TV term is used to detect sharp jumps of the function, and the Gaussian distribution is used as a reference measure so that it results in a well-defined posterior measure in the function space. We also present an efficient Markov Chain Monte Carlo (MCMC) algorithm to draw samples from the posterior distribution of the TG prior. With numerical examples we demonstrate the performance of the TG prior and the efficiency of the proposed MCMC algorithm.

An important task of uncertainty quantification is to identify {the probability of} undesired events, in particular, system failures, caused by various sources of uncertainties. In this work we consider the construction of Gaussian {process} surrogates for failure detection and failure probability estimation. In particular, we consider the situation that the underlying computer models are extremely expensive, and in this setting, determining the sampling points in the state space is of essential importance. We formulate the problem as an optimal experimental design for Bayesian inferences of the limit state (i.e., the failure boundary) and propose an efficient numerical scheme to solve the resulting optimization problem. In particular, the proposed limit-state inference method is capable of determining multiple sampling points at a time, and thus it is well suited for problems where multiple computer simulations can be performed in parallel. The accuracy and performance of the proposed method is demonstrated by both academic and practical examples.

In this work we consider a class of uncertainty quantification problems where the system performance or reliability is characterized by a scalar parameter $y$. The performance parameter $y$ is random due to the presence of various sources of uncertainty in the system, and our goal is to estimate the probability density function (PDF) of $y$. We propose to use the multicanonical Monte Carlo (MMC) method, a special type of adaptive importance sampling algorithm, to compute the PDF of interest. Moreover, we develop an adaptive algorithm to construct local Gaussian process surrogates to further accelerate the MMC iterations. With numerical examples we demonstrate that the proposed method can achieve several orders of magnitudes of speedup over the standard Monte Carlo method.

Many scientific and engineering problems require to perform Bayesian inferences in function spaces, in which the unknowns are of infinite dimension. In such problems, many standard Markov Chain Monte Carlo (MCMC) algorithms become arbitrary slow under the mesh refinement, which is referred to as being dimension dependent. In this work we develop an independence sampler based MCMC method for the infinite dimensional Bayesian inferences. We represent the proposal distribution as a mixture of a finite number of specially parametrized Gaussian measures. We show that under the chosen parametrization, the resulting MCMC algorithm is dimension independent. We also design an efficient adaptive algorithm to adjust the parameter values of the mixtures from the previous samples. Finally we provide numerical examples to demonstrate the efficiency and robustness of the proposed method, even for problems with multimodal posterior distributions.