Bayesian inference provides a systematic and rational framework for
learning complex physical models from complex observational data. In
the resulting Bayesian inverse problem (BIP), we employ observational
data and a model mapping parameters to observations to infer unknown
parameters and quantify their uncertainty, in the form of a posterior
probability distribution. But this begs the question: how do we choose
the data in the first place? The problem of choosing the optimal
training data for a BIP can be formulated as the following optimal
experimental design (OED) problem: optimize the design of a data
acquisition system (sensor locations, what is measured and how often,
etc.) to minimize some measure of the uncertainty in the inferred
parameters.
This constitutes a challenge of the highest order. First, the
parameters to be inferred in the inverse problem are often infinite
dimensional spatially correlated fields, whose discretization yields
high dimensional parameter spaces that can number in the
millions. Second, the physical models to be learned are often complex
PDE systems that are computationally expensive to solve. Third, the
BIP---difficult as it is---is merely an inner problem within the outer
OED problem, to be solved at each optimization iteration. And fourth,
efficient optimization algorithms for the OED problem require
gradients of the OED objective function with respect to the design
parameters, meaning that we must differentiate through the BIP
solution.
An attractive objective function for OED is the expected information
gain (EIG), which can be understood as the expected Kullback--Leibler
divergence between prior and posterior. Naive evaluation of EIG is
intractable for large-scale problems due to the large number of
samples required in double-loop Monte Carlo. Here we invoke an
approximation of the EIG based on the Laplace approximation of the
posterior. Each evaluation of the objective function then requires
computing a sample average approximation (over possible realizations
of the data) of the information gain (IG) between the Laplace
approximation and the Gaussian prior distribution. Since the IG
involves the KL divergence between two Gaussians, it can be expressed
in closed form, in particular in terms of the log-determinant and
trace of the prior-preconditioned posterior covariance operator. A
randomized eigensolver allows us to estimate these invariants at a
cost---measured in number of forward PDE solves---that is independent
of the dimension of the parameter space, resulting in scalable
evaluation of the OED objective. Variational adjoint methods are then
used to efficiently compute the gradient of the OED objective function
with respect to the experimental design variables, thus leading to an
efficient and scalable Laplace-approximation based OED
method. Moreover, the resulting method inherits the fine-grained
parallel efficiency of the forward PDE solver, despite being tightly
coupled. We present numerical results demonstrating the scalability
and efficiency of this approach. In particular, the results
demonstrate that the number of required PDE solves is independent of
the parameter dimension, and only weakly dependent on the data
dimension.
This work is joint with Umberto Villa.
Copyright. All Rights Reserved.