# Bayesian Level Set Approach for Inverse Problems with Piecewise Constant
Reconstructions

William Reese1, Arvind K. Saibaba1, Jonghyun Lee2 1 Department of Mathematics,
North Carolina State University, Raleigh NC 27695, USA  2 Department of Civil
and Environmental Engineering and Water Resources Research Center Honolulu
Hawaii 96822, USA

###### Abstract

There are several challenges associated with inverse problems in which we seek
to reconstruct a piecewise constant field, and which we model using multiple
level sets. Adopting a Bayesian viewpoint, we impose prior distributions on
both the level set functions that determine the piecewise constant regions as
well as the parameters that determine their magnitudes. We develop a Gauss-
Newton approach with a backtracking line search to efficiently compute the
maximum a priori (MAP) estimate as a solution to the inverse problem. We use
the Gauss-Newton Laplace approximation to construct a Gaussian approximation
of the posterior distribution and use preconditioned Krylov subspace methods
to sample from the resulting approximation. To visualize the uncertainty
associated with the parameter reconstructions we compute the approximate
posterior variance using a matrix-free Monte Carlo diagonal estimator, which
we develop in this paper. We will demonstrate the benefits of our approach and
solvers on synthetic test problems (photoacoustic and hydraulic tomography,
respectively a linear and nonlinear inverse problem) as well as an application
to X-ray imaging with real data.

††: Inverse Problems

Keywords: Bayesian Inverse Problems, Level set method, Uncertainty
Quantification, Matrix-free methods, Iterative solvers.

##  1 Introduction

Inverse problems allow us to estimate physical quantities that may not be
directly observable [2, 37, 19, 21]. Two examples that we consider in this
paper are X-ray tomography and hydraulic tomography. X-ray tomography is a non
invasive technique of to image the inside of a body. X-ray sources are placed
around the object and then radiation passes through the sources to
corresponding detectors. The attenuation of the X-rays is used to reconstruct
the internal details of the body [21]. In hydraulic tomography a series of
pumping tests is performed with multiple pumping and observation wells to
determine complex subsurface structures such as spatially distributed
hydraulic conductivity of an aquifer [24]. In both of these examples the
unknown quantity (an image or hydraulic conductivity field) may be represented
as a piecewise constant spatial field while the negligible intra-variability
can be ignored. Beyond these two examples, inverse problems have applications
in many areas of science and engineering such as biomedicine, geosciences, and
material science.

In deterministic inverse problems, the unknown parameter field is uniquely
determined such that the estimated field can reproduce the observations within
a reasonable error level. In addition measurements may be sparse or limited
and the recovery of the unknown parameter field from the measurements may be
an ill-posed problem. Such is the case, for example, in limited angle X-ray
tomography or subsurface structure characterization in hydraulic tomography.
Uncertainty Quantification (UQ) allows us to address these different
challenges. In deterministic inverse problems the ill-posedness is often
alleviated by solving a regularized optimization problem in which the
regularization term (e.g., Tikhonov regularization) often penalizes undesired
behavior in the reconstructions. Alternatively, the Bayesian approach assumes
the unknown field is random and characterizes the prior beliefs as a prior
distribution. The prior distribution is combined with the likelihood
distribution using the Bayes rule to produce the posterior distribution which
serves as the solution to the Bayesian inverse problem. In this paper, we are
primarily interested in inverse problems in which the unknown fields are
assumed to be piecewise constant. For example in hydraulic tomography, one
would like to identify a subsurface permeability field with discrete
geological zones or facies, such as a regional aquifer formed with deposited
layers of sediment, that are relatively homogeneous with each zone but with
abrupt changes at the boundaries between zones. Widely used methods Tikhonov
regularization or Gaussian-priors require a large amount of data and have high
computational costs associated with the data size to provide sharply changing
parameter estimation. This motivates advanced, scalable techniques to solve
inverse problems with piecewise constant fields.

The literature on inverse problems with piecewise constant unknowns is vast
and so we only mention a few key references. Perhaps the most common is the
use of total variation (TV) [19, 37, 21, 2]. TV regularization is favored over
Tikhonov regularization when the unknown parameter is discontinuous because it
preserves sharp interfaces as opposed to smoothing them. A difficulty of using
the TV regularization is the fact that it is non differentiable; however,
there are many ways this can be handled (e.g., using smoothed TV
approximations). The TV approach can also be used in the Bayesian setting by
incorporating the TV regularization as a prior information (this can
correspond to the Laplace prior). However, in [22], the authors demonstrate
that the conditional mean estimates for the TV prior are not edge preserving
through fine discretizations of the model space. There are other approaches
for promoting edge preserving behavior such as Besov priors [4], level set
approaches [7, 20, 14], and Cauchy difference priors [27]. In this paper, we
develop a multiple level set based approach for enforcing piecewise constant
reconstructions. The multiple level set approach utilizes a combination of
level set functions, that determine the shapes, and constants, representing
the magnitudes, to model the unknown parameter field. The advantage of level
set functions is their ability to model the discontinuous parameter field in
terms of a differentiable function. This changes the inversion problem from
estimating a piecewise constant field into a few (possibly smooth) level set
functions that determine the regions and the magnitudes of the field in those
regions.

##### Overview of the main contributions.

This paper proposes a Bayesian level set approach for solving inverse problems
in which the unknown parameters are represented as piecewise constant fields.
There are three main contributions of this paper:

  1. 1.

We use multiple level sets to enforce piecewise constant reconstructions. This
approach allows the flexibility to incorporate prior information for the level
set functions $\boldsymbol{\Phi}$ which control the shape and the parameters
governing the magnitudes of the representations $\mathbf{c}$. We discuss
several possible choices for prior distributions.

  2. 2.

Assuming a Gaussian distribution for the measurement errors and the prior
distribution, we derive the posterior distribution, which, in general is non-
Gaussian because of the nonlinear parameterization. We develop an inexact
Gauss-Newton approach to obtain the maximum a posteriori (MAP) estimate. The
algorithm is efficient in that it can take advantage of the adjoint-based
method for computing derivatives and only incurs a small additional cost due
to the level set parameterization.

  3. 3.

To quantify the uncertainty in the reconstructed parameters, we approximate
the posterior distribution using Laplace’s approximation based on the Gauss-
Newton Hessian; this is known as linearized Bayesian inference (see e.g.,
[36]). Using this approximate distribution we develop matrix-free (here and
henceforth, by matrix-free, we mean that the approaches only rely on matrix-
vector products and do not require us to form the operators explicitly)
iterative methods to generate samples using the preconditioned Lanczos
approach and estimate the posterior variance using two different methods. As
with the MAP estimate, the posterior samples also have the piecewise constant
structure. We also develop a matrix-free estimator for the diagonals of the
inverse of a matrix (LanczosMC) that combines the Lanczos and Monte Carlo
approaches, and which may be of independent interest beyond this paper.

Finally, we demonstrate the benefits of the Bayesian level set approach and
the computational benefits of the proposed algorithms on a range of
applications such as Photoacoustic Tomography, X-ray Tomography using a real-
world data set, and hydraulic tomography (which is a nonlinear inverse
problem).

##### Related approaches.

We position our work in relation to other works on level methods in inverse
problems. Surveys of level set methods in deterministic inverse problems are
provided in [6, 11]. Another deterministic approach that uses level sets is
parametric level set approach (PaLS) [1]. In PaLS, the level set function(s)
are represented a weighted sum of radial basis functions and the unknown
parameters are the weights. The advantage of this approach over our multiple
level set approach is the use of a few parameters (e.g., on the order of
$10-100$) rather than a pixel or voxel-wise approach that can require millions
of unknowns. Our work is closely related to [35], who also use a multiple
level set approach and a Gauss-Newton type approach for solving the
optimization problem. However, quantifying the reconstruction uncertainty is
still an active area of research.

Several authors have utilized level sets to perform Bayesian inversion for
piecewise constant parameters. Cardiff and Kitanidis [7] demonstrate the use
of level set functions for facies detection in a hydraulic tomography problem
and provide uncertainty estimates for the facies’ boundaries. This approach
relies on computing of the velocity field for the level set functions to
update the interfaces. In [20] the authors formulate a Bayesian level set
inverse problem in terms of a single level set function, but assume that the
number of levels and their magnitudes are known a priori. The authors
formulate the Bayesian inverse problem in an infinite dimensional setting and
utilize Monte Carlo Markov Chain (MCMC) methods to sample from the posterior
distribution. This approach does not require computing the velocity field for
the level set function evolution as it is updated through the sampling. In
[13] the authors extend the previous approach by incorporating length and
amplitude scales as parameters in addition to the level set function in the
prior. In contrast to these approaches, which use MCMC techniques to explore
the posterior distribution, our approach uses matrix-free techniques that are
scalable to large problem sizes (both in the number of unknown parameters and
measurements).

##### Overview of the paper.

We conclude this section with a brief overview of this paper. In Section 2, we
review the basics of inverse problems and the level set approach. In Section
3, we present the Bayesian level set approach, discuss several options for
prior information, and use an inexact Gauss-Newton approach for computing the
maximum a posteriori (MAP) estimate. In Section 4, we derive an expression for
the approximate posterior distribution and develop matrix-free techniques to
generate samples from this distribution and the estimate the posterior
variance. In Section 5, we apply this framework and the solvers we have
developed to synthetic test problems from photoacoustic and hydraulic
tomography, and real-world problem from X-ray tomography. In A, we derive a
matrix-free technique to estimate the diagonals of the inverse of a matrix and
validate on several different matrices; we dissociate it from the rest of the
paper since it may be of interest beyond this paper. Finally, we offer some
concluding remarks and future directions in Section 6.

##  2 Background

In this section, we give background information on the Bayesian inverse
problem and level set technique for representing piecewise constant functions.

###  2.1 Bayesian Inverse Problem

Let $\mathbf{f}\in\mathbb{R}^{n}\rightarrow\mathbb{R}^{n_{\rm obs}}$ represent
the forward operator or the parameter-to-observable map, let
$\mathbf{m}\in\mathbb{R}^{n}$ denote the unknown parameters to be recovered.
If we assume that the measurement are corrupted by additive Gaussian noise
$\boldsymbol{\varepsilon}$, then the following equation relates the
measurements $\mathbf{d}\in\mathbb{R}^{n_{\rm obs}},$ to the unknown
parameters

$\mathbf{d}=\mathbf{f}(\mathbf{m})+\boldsymbol{\varepsilon}\;\;\;\;\boldsymbol{\varepsilon}\sim\mathcal{N}(\mathbf{0},\boldsymbol{\Gamma}_{\text{noise}}).$
(1)

Here, we represent the measurement error $\boldsymbol{\varepsilon}$ as a
Gaussian distribution with zero mean, and a positive definite covariance
matrix $\boldsymbol{\Gamma}_{\text{noise}}$. The inverse problem can be stated
as follows: given measurements $\mathbf{d}$ use (1) to recover the unknown
parameters $\mathbf{m}$.

Fundamental limitations in data acquisition, means that the recovery of
spatial parameters $\mathbf{m}$ is mathematically an ill-posed problem: there
may be no solution, the solution may not be unique, or it may depend
sensitively on the data $\mathbf{d}$. The Bayesian approach specifies a prior
probability distribution $\pi_{\mathrm{pr}}(\mathbf{m})$ that describes the
prior information, or “belief,” about the parameter $\mathbf{m}$. Then Bayes’
rule is applied to derive the posterior probability distribution, defined as
the distribution of the parameter $\mathbf{m}$ conditioned on the data
$\mathbf{d}$ and is represented by
$\pi_{\text{post}}(\mathbf{m}\mid\mathbf{d})$. The posterior distribution
takes the form

$\pi_{\text{post}}(\mathbf{m}\mid\mathbf{d})\propto\>\pi_{\text{like}}(\mathbf{d}\mid\mathbf{m})\pi_{\text{pr}}(\mathbf{m})\propto\exp\left(-\frac{1}{2}\|\mathbf{d}-\mathbf{f}(\mathbf{m})\|_{\boldsymbol{\Gamma}_{\text{noise}}^{-1}}^{2}\right)\pi_{\text{pr}}(\mathbf{m}).$
(2)

A standard approach is to compute the so-called maximum a posteriori (MAP)
estimate. That is, the MAP point $\mathbf{m}_{\text{MAP}}$ maximizes the
posterior distribution $\pi_{\text{post}}(\mathbf{m}\mid\mathbf{d})$. If the
prior distribution is of the form
$\pi_{\text{prior}}(\mathbf{m})\propto\exp(-\mathcal{R}(\mathbf{m}))$, the MAP
estimate is the solution to the variational problem

$\min_{\mathbf{m}\in\mathbb{R}^{n}}\>\frac{1}{2}\|\mathbf{d}-\mathbf{f}(\mathbf{m})\|_{\boldsymbol{\Gamma}_{\text{noise}}^{-1}}^{2}+\mathcal{R}(\mathbf{m}).$
(3)

The MAP estimate is easy to interpret as the “most likely” value of the
parameter $\mathbf{m}$, given the data $\mathbf{d}$. The MAP estimate can be
computed using optimization methods and the specific choice of the solver
(either a general purpose or a special purpose) depends on the choice of the
prior distribution [34].

###  2.2 Level Set Representation

To explain the idea behind the level set method, consider a piecewise constant
function $m$ that only takes two positive values $c_{p}$ and $c_{b}$, in the
domain $\Omega\subset{}^{d}$; here the dimension takes values $d=1,2,3$. The
function can be represented mathematically as follows. Let $A\subset\Omega$ be
a closed subset and let $\partial\Omega$ denote its boundary. Let $\chi_{A}$
denote an indicator function, which takes two values $1$ and $0$, depending on
whether the point $\mathbf{r}\in A$ or outside of $A$, respectively. Then, the
function can be written in terms of the indicator function as

$m(\mathbf{r})=c_{p}\chi_{A}(\mathbf{r})+c_{b}[1-\chi_{A}(\mathbf{r})]\qquad\forall\mathbf{r}\in\Omega.$

Then, in order to completely specify the function $m$, we only need knowledge
of the region $A$. The level set approach is shape-based and represents the
boundary $\partial A$ as the level set of a function
$\phi:\Omega\rightarrow\real$. Specifically, if we take the zero level set,
then $\phi(\mathbf{r})$ is related to $A$ and $\partial A$ as

$\begin{array}[]{ll}\phi(\mathbf{r})>0&\forall\mathbf{r}\in A\\\
\phi(\mathbf{r})=0&\forall\mathbf{r}\in\partial A\\\
\phi(\mathbf{r})<0&\forall\mathbf{r}\in\Omega\backslash A.\end{array}$

Alternatively, we can represent the function $m(\mathbf{r})$ in terms of the
Heaviside function $H(x)=\frac{1}{2}(1+\text{sign}(x))$ as

$m(\mathbf{r})=c_{p}H(\phi(\mathbf{r}))+c_{b}[1-H(\phi(\mathbf{r}))].$

By this formulation, identifying the unknown region $A$, is tantamount to
estimating the continuous function $\phi(\mathbf{r})$ and the coefficients
$c_{p}$ and $c_{b}$. To represent $m(\mathbf{r})$ when it has 4 regions we can
use 2 level set functions and write

$\displaystyle m(\mathbf{r})$
$\displaystyle=c_{1}(1-H(\boldsymbol{\phi}_{1}))(1-H(\boldsymbol{\phi}_{2}))+c_{2}H(\boldsymbol{\phi}_{1})(1-H(\boldsymbol{\phi}_{2}))$
$\displaystyle+c_{3}(1-H(\boldsymbol{\phi}_{1}))H(\boldsymbol{\phi}_{2})+c_{4}H(\boldsymbol{\phi}_{1})H(\boldsymbol{\phi}_{2})$

We illustrate this representation for two level set functions in Figure 1.

Figure 1: Left: piecewise constant field with 4 regions, Right: level set
functions $\phi_{1}$ and $\phi_{2}$ that are used to compute $m$.

We now explain how to extend this approach to recover piecewise constant
functions with more than $2$ regions. This follows the discussion in [35].
Suppose that there are $N_{\text{R}}$ regions to be recovered. We choose the
number of level set functions $N_{\text{LS}}$ to be such that

$2^{N_{\text{LS}}-1}<N_{\text{R}}\leq 2^{N_{\text{LS}}}.$ (4)

Denote the corresponding level set functions $\phi_{i}(\mathbf{r})$ for $1\leq
i\leq N_{\text{LS}}$. Then, the unknown function $f(\mathbf{r})$ can be
expressed in terms of the unknown level set functions as

$m(\mathbf{r})=\sum_{i_{1},\dots,i_{N_{\text{LS}}}=\\{0,1\\}}c_{i_{1},\dots,i_{N_{\text{LS}}}}\widehat{H}_{i_{1}}(\phi_{1}(\mathbf{r}))\cdots\widehat{H}_{i_{N_{\text{LS}}}}(\phi_{N_{\text{LS}}}(\mathbf{r})),$
(5)

where $c_{i_{1},\dots,i_{N_{\text{LS}}}}$ are the magnitudes of the piecewise
constant regions, and summation involves $2^{N_{\text{LS}}}$ terms.
Furthermore, the modified level set functions $\widehat{H}_{i}(\cdot)$ have
the form

$\widehat{H}_{i}(x)=\left\\{\begin{array}[]{ll}H(x)&i=0\\\
1-H(x)&i=1.\end{array}\right.$

In this formulation, recovering the unknown function $m(\mathbf{r})$ now
involves recovering the level set functions $\phi_{i}(\mathbf{r})$ and the
coefficients $c_{i_{1},\dots,i_{N_{\text{LS}}}}$. In practice, the Heaviside
function $H(x)=\frac{1}{2}(1+\text{sign}(x))$ is not differentiable; to
address this issue we replace it with the mollified Heaviside function
$H_{\epsilon}$ as suggested in [1]:

$H_{\epsilon}(x)=\left\\{\begin{array}[]{ll}1,&x>\epsilon,\\\
\frac{1}{2}+\frac{x}{2\epsilon}+\frac{1}{2\epsilon}\sin(\frac{x\pi}{\epsilon}),&|x|\leq\epsilon,\\\
0,&x<-\epsilon,\end{array}\right.$ (6)

where $\epsilon$ is a positive parameter used to control the resolution of
$H_{\epsilon}$ at the point $x=0$. This choice of the Heaviside function
ensures sparsity in the Jacobian which results in lower computational costs.

Figure 2: (left) Heaviside function $H(x)$ and its approximation
$H_{\epsilon}(x)$, (right) derivative of Heaviside function
$\delta_{\epsilon}(x)$.

##### Discretization

The unknown coefficients $c_{i_{1},\dots,i_{N_{\text{LS}}}}$ can be rearranged
in column major order into a vector $\mathbf{c}$ of the size
$2^{N_{\text{LS}}}$. The domain $\Omega$ is discretized into $N$ points, and
the vectors
$\boldsymbol{\phi}^{(1)},\dots,\boldsymbol{\phi}^{(N_{\text{LS}})}$ represent
the discretized versions with respect to some basis. Concatenating the vectors
$\boldsymbol{\phi}^{(i)}$ into the vector

$\boldsymbol{\Phi}=\begin{bmatrix}\boldsymbol{\phi}^{(1)}\\\ \vdots\\\
\boldsymbol{\phi}^{(N_{\text{LS}})}\end{bmatrix}\in{}^{n{N_{\text{LS}}}}$

and utilizing (5) we can write the unknown vector $\mathbf{m}$ as

$\mathbf{m}=\mathbf{m}(\boldsymbol{\Phi},\mathbf{c})$ (7)

i.e., $\mathbf{m}$ is a vector valued function
$\mathbf{m}:{}^{(nN_{\text{LS}})}\times{}^{2^{N_{\text{LS}}}}\longrightarrow\mathbb{R}^{n}$.
Now that we have a way to present piecewise constant images we now address the
problem of estimating the parameters $\boldsymbol{\Phi}$ and $\mathbf{c}$,
which describe the image, from a set of measurements $\mathbf{d}$. To this
end, we adopt the Bayesian approach described in [21, 34].

##  3 Bayesian Level Set Formulation and MAP Estimate

With the parameterization of $\mathbf{m}$ as defined in (7), we can rewrite
the measurement equation as

$\mathbf{d}=\mathbf{f}(\mathbf{m}(\boldsymbol{\Phi},\mathbf{c}))+\boldsymbol{\varepsilon}.$
(8)

As in Section 2.1, we assume that $\boldsymbol{\varepsilon}$ is Gaussian with
a known mean and covariance
$\boldsymbol{\varepsilon}\sim\mathcal{N}(\mathbf{0},\boldsymbol{\Gamma}_{\text{noise}})$.
This equation defines the likelihood

$\pi(\mathbf{d}|\boldsymbol{\Phi},\mathbf{c})\propto\exp\left(-\frac{1}{2}\left\lVert\mathbf{f}(\mathbf{m}(\boldsymbol{\Phi},\mathbf{c}))-\mathbf{d}\right\rVert_{\boldsymbol{\Gamma}_{\text{noise}}^{-1}}^{2}\right).$

We assume that our unknown parameters $\boldsymbol{\Phi}$ and $\mathbf{c}$ are
independent random vectors and our goal is to estimate these unknown
parameters, given measurements $\mathbf{d}$. The posterior distribution
combines the likelihood and the prior information through Bayes’ theorem.
However, the specific form of the posterior distribution depends on the
specific choice of the prior distribution. We briefly discuss various possible
choices for the prior distributions.

###  3.1 Choice of prior distributions

Since our unknowns consist of the coefficients $\mathbf{c}$ and the level set
functions $\boldsymbol{\Phi}$ it is reasonable to assume that they are
independently distributed random vectors. Thus we will discuss possible
choices of prior for $\mathbf{c}$ and $\boldsymbol{\Phi}$ separately.

####  3.1.1 Prior for the level set functions $\boldsymbol{\Phi}$

Following the level set formulation it is apparent that we want each of the
level set functions $\boldsymbol{\phi}_{i}$ to be “smooth” for $1\leq i\leq
N_{\text{LS}}$. One approach to enforce smoothness in prior distributions is
to use Gaussian distribution, which we adopt in this paper; it takes the form

$\boldsymbol{\phi}\sim\mathcal{N}(\boldsymbol{\mu}_{\boldsymbol{\phi}},\lambda_{\boldsymbol{\phi}}^{-2}\boldsymbol{\Gamma}_{\boldsymbol{\phi}})$
(9)

where $\lambda_{\boldsymbol{\phi}}^{-2}$ is the precision and
$\boldsymbol{\Gamma}_{\boldsymbol{\phi}}\in\mathbb{R}^{n\times n}$ is the
covariance matrix and $\boldsymbol{\mu}_{\boldsymbol{\phi}}$ is the mean. To
assign a prior for $\boldsymbol{\Phi}$ we assume that each level set function
is independently distributed and each level set function has the distribution
given by (9). In other words we assume
$\boldsymbol{\phi}_{i}\sim\mathcal{N}(\boldsymbol{\mu}_{\boldsymbol{\phi}},\lambda_{\boldsymbol{\Phi}}^{-2}\boldsymbol{\Gamma}_{\boldsymbol{\phi}})$,
are independent for $1\leq i\leq N_{\text{LS}}.$ Thus, we can write the prior
for $\boldsymbol{\Phi}$ as

$\boldsymbol{\Phi}\sim\mathcal{N}(\boldsymbol{\mu}_{\boldsymbol{\Phi}},\lambda_{\boldsymbol{\Phi}}^{-2}\boldsymbol{\Gamma}_{\boldsymbol{\Phi}})$
(10)

where the mean
$\boldsymbol{\mu}_{\boldsymbol{\Phi}}=\mbox{\leavevmode\hbox{\small
1\kern-3.3pt\normalsize 1}}\otimes\boldsymbol{\mu}_{\boldsymbol{\phi}}$ (here
11= $[{1,\dots,1}]^{T}\in\mathbb{R}^{N_{\text{LS}}}$), and the covariance
matrix $\boldsymbol{\Gamma}_{\boldsymbol{\Phi}}$ takes the form
$\boldsymbol{\Gamma}_{\boldsymbol{\Phi}}=\mathbf{I}_{N_{\text{LS}}}\otimes\boldsymbol{\Gamma}_{\boldsymbol{\phi}}.$
The symbol $\otimes$ represents the Kronecker product. We assume that each
level set function is independent and Gaussian with mean
$\boldsymbol{\mu}_{\boldsymbol{\phi}}\in\mathbb{R}^{n}$ and covariance
$\boldsymbol{\Gamma}_{\boldsymbol{\phi}}\in\mathbb{R}^{n\times n}$, so that we
assume the same regularization parameter
$\lambda_{\boldsymbol{\Phi}}=\lambda_{\boldsymbol{\phi}}$ for each level set.
This reduces the number of hyperparameters to be estimated.

There are several possible choices for the covariance matrix
$\boldsymbol{\Gamma}_{\boldsymbol{\phi}}$. If we assume that the underlying
field $\phi$ is a Gaussian random field, then the covariance matrix
$\boldsymbol{\Gamma}_{\boldsymbol{\phi}}$ can be chosen based on a kernel
function. There are several choices here, e.g., Matérn class,
$\gamma$-exponential, etc [31]. Another approach is to choose priors from the
Whittle-Matérn class as was done in [23]. Our choice of covariance matrix
$\boldsymbol{\Gamma}_{\boldsymbol{\phi}}$ will depend on the test problem. Yet
another choice, which we adopt in this paper, is to use the Gauss-Markov
Random field (GMRF) to specify the covariance matrix [2, 32]. We choose the
precision matrices (inverse of the covariance matrix
$\boldsymbol{\Gamma}_{\boldsymbol{\phi}}$) as the discretization of a partial
differential operator. This has close connections with the Gaussian random
field approach [26].

####  3.1.2 Prior for level set magnitudes $\mathbf{c}$

The elements of $\mathbf{c}$ represent the magnitudes of the values taken by
$\mathbf{m}$. When $\mathbf{m}$ represents a piecewise constant image, the
elements $c_{i}$, $1\leq i\leq 2^{N_{\text{LS}}}$ assign intensity of pixels
belonging to different regions. Therefore, imposing a prior distribution for
$\mathbf{c}$ controls the distribution of the intensity of pixels. We briefly
mention a few choices for prior distributions. A convenient representation for
the prior information for the intensities is the multivariate Gaussian
distribution, i.e.,

$\mathbf{c}\sim\mathcal{N}(\boldsymbol{\mu}_{\mathbf{c}},\boldsymbol{\Gamma}_{\mathbf{c}})$
(11)

where $\boldsymbol{\mu}_{\mathbf{c}}\in{}^{2^{N_{\text{LS}}}}$ and
$\boldsymbol{\Gamma}_{\mathbf{c}}\in{}^{{2^{N_{\text{LS}}}}\times{2^{N_{\text{LS}}}}}$
are the mean and covariance respectively. This is the approach we adopt in
this paper.

Alternatively, if one requires that $\mathbf{c}$ have positive elements, a
lognormal or gamma prior may be more appropriate. In the lognormal case, we
assume that
$\mathbf{y}=\log\mathbf{c}\sim\mathcal{N}(\boldsymbol{\mu}_{c},\boldsymbol{\Gamma}_{\mathbf{c}})$
(where the log is computed elementwise). Finally, if we assume the pixel
intensities are independent and identically distributed according to the gamma
distribution with its shape and (inverse) scale parameters $\alpha>0$ and
$\beta>0$. That is, $c_{i}\sim\text{Gamma}(\alpha,\beta),1\leq i\leq
N_{\text{LS}}.$ This implies that the pdf of the prior distribution satisfies

$\pi(\mathbf{c})\propto\prod_{i=1}^{2^{N_{\text{LS}}}}\Bigg{(}c_{i}^{\alpha-1}\text{exp}(-\beta
c_{i})\Bigg{)}.$

###  3.2 Efficient solver for computing the MAP estimate

For the rest of the section, we assume that the prior distribution for
$\boldsymbol{\Phi}$ and $\mathbf{c}$ are Gaussian. That is, we assume the
following model

$\boldsymbol{\Phi}\sim\mathcal{N}(\boldsymbol{\mu}_{\boldsymbol{\phi}},\lambda^{-2}_{\boldsymbol{\Phi}}\boldsymbol{\Gamma}_{\boldsymbol{\phi}})\qquad\mathbf{c}\sim\mathcal{N}(\boldsymbol{\mu}_{\mathbf{c}},\lambda^{-2}_{\mathbf{c}}\boldsymbol{\Gamma}_{\mathbf{c}}).$
(12)

With this choice of prior distribution, the posterior distribution takes the
form

$\pi_{\text{post}}(\boldsymbol{\Phi},\mathbf{c}\mid\mathbf{d})\propto\exp\Big{(}-\frac{1}{2}\left\lVert\mathbf{f}(\mathbf{m}(\boldsymbol{\Phi},\mathbf{c}))-\mathbf{d}\right\rVert_{\boldsymbol{\Gamma}_{\text{noise}}^{-1}}^{2}-\frac{\lambda_{\boldsymbol{\Phi}}^{2}}{2}\left\lVert\boldsymbol{\Phi}-\boldsymbol{\mu}_{\boldsymbol{\Phi}}\right\rVert_{\boldsymbol{\Gamma}_{\boldsymbol{\Phi}}^{-1}}^{2}-\frac{\lambda^{2}_{\mathbf{c}}}{2}\left\lVert\mathbf{c}-\boldsymbol{\mu}_{\mathbf{c}}\right\rVert_{\boldsymbol{\Gamma}_{\mathbf{c}}{}^{-1}}^{2}\Big{)}.$
(13)

It is worth noting that due to the parameterization, the posterior
distribution is non-Gaussian, even though the measurement noise and the prior
distributions are Gaussian. This is due to the nonlinearity in the forward
operator as well as the nonlinearity in the level set representation
$\mathbf{m}(\boldsymbol{\Phi},\mathbf{c})$.

The MAP point is a maximizer of the posterior distribution; alternatively, it
can be obtained as the minimizer of the negative logarithm of the posterior
distribution. More specifically, the MAP point can be found by solving the
optimization problem

$\min_{\boldsymbol{\Phi},\mathbf{c}}\>\left[-\log\pi_{\text{post}}(\boldsymbol{\Phi},\mathbf{c}\mid\mathbf{d})\right]=\min_{\boldsymbol{\Phi},\mathbf{c}}\textrm{F}(\boldsymbol{\Phi},\mathbf{c})$
(14)

where

$\textrm{F}(\boldsymbol{\Phi},\mathbf{c})=\frac{1}{2}\left\lVert\mathbf{f}(\mathbf{m}(\boldsymbol{\Phi},\mathbf{c}))-\mathbf{d}\right\rVert_{\boldsymbol{\Gamma}_{\text{noise}}^{-1}}^{2}+\frac{\lambda^{2}_{\boldsymbol{\Phi}}}{2}\left\lVert\boldsymbol{\Phi}-\boldsymbol{\mu}_{\boldsymbol{\Phi}}\right\rVert_{\boldsymbol{\Gamma}_{\boldsymbol{\Phi}}^{-1}}^{2}+\frac{\lambda^{2}_{\mathbf{c}}}{2}\left\lVert\mathbf{c}-\boldsymbol{\mu}_{\mathbf{c}}\right\rVert_{\boldsymbol{\Gamma}_{\mathbf{c}}{}^{-1}}^{2}.$
(15)

Since (14) is a nonlinear least squares problem we solve this efficiently
using an inexact Gauss-Newton approach.

For simplicity, we denote the concatenated vector
$\mathbf{x}=\begin{bmatrix}\boldsymbol{\Phi}^{T}&\mathbf{c}^{T}\end{bmatrix}^{T}.$
Using the chain rule, the Jacobian of $\mathbf{f}$ with respect to
$\mathbf{x}$ is

$\mathbf{J}=\frac{\partial\mathbf{f}}{\partial\mathbf{m}}\frac{\partial\mathbf{m}}{\partial\mathbf{x}}.$
(16)

The term $\frac{\partial\mathbf{f}}{\partial\mathbf{m}}$ is the part of the
Jacobian that comes from the forward operator and is problem dependent. It is
important to note that this matrix need not be computed explicitly but we only
need to be able to compute its action, i.e., compute matrix-vector products
with this matrix. The Jacobian matrix of $\mathbf{m}$ with respect to
$\mathbf{x}$ has the following form

$\frac{\partial\mathbf{m}}{\partial\mathbf{x}}=\begin{bmatrix}\frac{\partial\mathbf{m}}{\partial\boldsymbol{\phi}_{1}}&\dots&\frac{\partial\mathbf{m}}{\partial\boldsymbol{\phi}_{N_{\text{LS}}}}&\frac{\partial\mathbf{m}}{\partial\mathbf{c}}\end{bmatrix}$
(17)

where the matrices

$\frac{\partial\mathbf{m}}{\partial\boldsymbol{\phi}_{i}}\in\mathbb{R}^{n\times
n}\qquad 1\leq i\leq
N_{\text{LS}}\qquad\text{and}\quad\frac{\partial\mathbf{m}}{\partial\mathbf{c}}\in\mathbb{R}^{n\times
2^{N_{\text{LS}}}}.$ (18)

The details are given in [35].

It is worth noting that the matrix
$\frac{\partial\mathbf{m}}{\partial\mathbf{x}}$ has a lot of sparsity that can
be exploited for computational advantages. The matrices
$\frac{\partial\mathbf{m}}{\partial\boldsymbol{\phi}_{i}}$ are diagonal and
can be stored in sparse format. There may be additional sparsity in the
diagonal entries and in ${\frac{\partial\mathbf{m}}{\partial\mathbf{c}}}$
depending on the support of the derivative of the mollified Heaviside
function. In Figure 3, we illustrate the structure of the Jacobian.

Figure 3: Visualization of the structure of the Jacobian.

To compute the MAP estimate, we solve the nonlinear least squares problem (14)
with an inexact Gauss-Newton approach [29, Chapter 10.2]. Given the current
estimate $\mathbf{x}_{k}$, the gradient of the objective function is given by

$\mathbf{g}_{k}=\mathbf{J}^{T}_{k}\boldsymbol{\Gamma}_{\text{noise}}^{-1}\left(\mathbf{f}(\mathbf{m}_{k})-\mathbf{d}\right)+\boldsymbol{\Gamma}_{\text{prior}}^{-1}(\mathbf{x}_{k}-\boldsymbol{\mu}_{\text{prior}})$
(19)

where $\mathbf{m}_{k}=\mathbf{m}(\boldsymbol{\Phi}_{k},\mathbf{c}_{k})$ and

$\boldsymbol{\mu}_{\text{prior}}=\begin{bmatrix}\boldsymbol{\mu}_{\boldsymbol{\Phi}}\\\
\boldsymbol{\mu}_{\mathbf{c}}\end{bmatrix}\qquad\boldsymbol{\Gamma}_{\text{prior}}=\begin{pmatrix}\lambda^{-2}_{\boldsymbol{\Phi}}\boldsymbol{\Gamma}_{\boldsymbol{\Phi}}&\\\
&\lambda^{-2}_{\mathbf{c}}\boldsymbol{\Gamma}_{\mathbf{c}}\end{pmatrix}.$ (20)

The Gauss-Newton step updates the current iterate as
$\mathbf{x}_{k+1}=\mathbf{x}_{k}+\alpha_{k}\delta\mathbf{x}_{k}$ where
$\alpha_{k}$ is a step size and the search direction $\delta\mathbf{x}_{k}$ is
obtained by solving the linear system of equations

$\mathbf{H}_{k}\delta\mathbf{x}_{k}=-\mathbf{g}_{k}\qquad\text{where}\qquad\mathbf{H}_{k}=\mathbf{J}_{k}^{T}\boldsymbol{\Gamma}_{\text{noise}}^{-1}\mathbf{J}_{k}+\boldsymbol{\Gamma}_{\text{prior}}^{-1}.$
(21)

The matrix $\mathbf{H}_{k}$ is known as the Gauss-Newton Hessian at step $k$.
The system of equations is solved using the conjugate gradient (CG) iterative
method. The stopping criterion for the CG iterations is
$\|\mathbf{H}_{k}\delta\mathbf{x}_{k}+\mathbf{g}_{k}\|_{2}\leq\eta_{k}\|\mathbf{g}_{k}\|_{2}$
where the residual tolerance is taken to be a decreasing sequence of the form

$\eta_{k}=\min\left(0.5,\sqrt{\frac{\|\mathbf{g}_{k}\|_{2}}{\|\mathbf{g}_{0}\|_{2}}}\right).$
(22)

The step length $\alpha_{k}$ is chosen using a backtracking line search
procedure in order to satisfy the strong Wolfe Conditions with parameters
$c_{1}=10^{-4}$, $c_{2}=0.9$ [29, Section 3.1]. To stop the inexact Gauss-
Newton iterations we use a combination of two different stopping criteria. The
first criterion is based on the Morozov discrepancy principle [37, Chapter 7]
and terminates the iterations if the data misfit falls below a specified
amount, whereas the second criterion is based on the first order optimality
condition and terminates the iterations if the gradient is smaller than a
user-defined tolerance.

###  3.3 Computational cost

To assess computational costs of each Gauss-Newton iteration we consider
storage of matrices as well as the number of matrix-vector products required.
In terms of storage, the additional storage required is $N_{\text{LS}}$
diagonal $n\times n$ matrices
${\frac{\partial\mathbf{m}}{\partial\boldsymbol{\phi}_{i}}}$ with sparse
diagonals, and one dense $n\times 2^{N_{\text{LS}}}$ matrix
${\frac{\partial\mathbf{m}}{\partial\mathbf{c}}}$; the total cost of storage
is $n(N_{\text{LS}}+2^{N_{\text{LS}}})$. Note that the matrix $\mathbf{J}_{k}$
need not be formed explicitly.

Let $N_{\text{CG}}$ represent the number of CG iterations in one Gauss-Newton
iteration. Let $T_{\text{der}}$ be the cost of applying the derivative
${\frac{\partial\mathbf{f}}{\partial\mathbf{m}}}$ or its adjoint, and let
$T_{\text{prior}}$ be the cost of applying the precision matrix
$\boldsymbol{\Gamma}_{\text{prior}}^{-1}$. Then, the cost of applying the
Jacobian $\mathbf{J}_{k}$ or its adjoint is

$T_{\text{jac}}=T_{\text{der}}+\mathcal{O}(n(N_{\text{LS}}+2^{N_{\text{LS}}}))\>\text{flops}.$

Since each Gauss-Newton iteration requires solving the system (21) we need to
compute the gradient $\mathbf{g}_{k}$ and matrix vector products with
$\mathbf{H}_{k}$. The gradient involves a matrix-vector product with the
adjoint of $\mathbf{J}_{k}$ and $\boldsymbol{\Gamma}_{\text{prior}}^{-1}$.
Each matrix vector product with $\mathbf{H}_{k}$ requires a matrix vector
product with $\boldsymbol{\Gamma}_{\text{prior}}^{-1}$, one with
$\mathbf{J}_{k}$ and its adjoint. Therefore, the total cost of a Gauss-Newton
step (not including step length computation) is

$T_{\text{GN}}=N_{\text{CG}}(2T_{\text{jac}}+T_{\text{prior}})+T_{\text{jac}}+T_{\text{prior}}\>\text{flops}.$

An important point to make here is that the use of the level set approach has
only a marginal increase in the computational cost per iteration and the
memory footprint. Furthermore, the proposed approach can be easily
incorporated into standard solvers for inverse problems, since it reutilizes
the ability to apply Jacobians and their adjoints.

##  4 Linearized Bayesian Inference

A common approach for uncertainty quantification is to generate independent
samples from the posterior distribution and use these samples along with Monte
Carlo estimates to generate statistics of the quantities of interest. Due to
the nonlinearity of the Heaviside function, the posterior distribution is, in
general, non-Gaussian even if the measurement operator
$\mathbf{f}(\mathbf{m})=\mathbf{A}\mathbf{m}$ is linear and the prior
distributions are Gaussian. This implies that the problem of generating
samples from the posterior distribution can be challenging. The prevalent
approach for sampling is some version of Markov Chain Monte Carlo approach.
However, the performance of MCMC approaches suffers in high dimensional spaces
and good proposal distributions are necessary for successful implementation of
MCMC methods. In this paper, we settle for linearized Bayesian inference
(following the terminology in [5, 36]) and approximate the posterior
distribution using the Laplace’s approximation (with the Gauss-Newton Hessian
instead of the full Hessian) . We develop matrix-free techniques for
generating samples from the approximate distribution and to estimate the
posterior variance.

A word on the notation: We let
$\mathbf{x}=\begin{bmatrix}\boldsymbol{\Phi}^{T}&\mathbf{c}^{T}\end{bmatrix}^{T}$
as the concatenated unknown vector and consider the mapping
$\mathbf{m}(\mathbf{x})$ as the mapping from the unknown vector to the pixel
space. With this notation, the posterior distribution can be written as
$\pi_{\mathrm{post}}(\mathbf{x}|\mathbf{d})$.

###  4.1 Sampling from the posterior distribution

We first derive an expression for the approximate posterior distribution by
using the Laplace’s approximation to the posterior distribution. This gives us
a Gaussian distribution centered about the MAP estimate
$\boldsymbol{\mu}_{\mathbf{x}}=[\boldsymbol{\Phi}_{\text{MAP}}^{T},\mathbf{c}_{\text{MAP}}^{T}]^{T}$
and with the covariance $\boldsymbol{\Gamma}_{\mathbf{x}}$, given by the
expression

$\boldsymbol{\Gamma}_{\mathbf{x}}=\>(\widehat{\mathbf{J}}^{T}\boldsymbol{\Gamma}_{\text{noise}}^{-1}\widehat{\mathbf{J}}+\boldsymbol{\Gamma}_{\text{prior}}^{-1})^{-1},$

where $\widehat{\mathbf{J}}$ is the Jacobian of $\mathbf{f}(\cdot)$ computed
at $\boldsymbol{\mu}_{\mathbf{x}}$. In other words, the covariance matrix is
taken to be the inverse of the Gauss-Newton Hessian computed at the MAP point,
i.e., $\boldsymbol{\Gamma}_{\mathbf{x}}=\mathbf{H}_{\text{GN}}^{-1}$. Note
that we are using the Gauss-Newton Hessian rather than the full Newton Hessian
since the former is guaranteed to be positive definite in the context of our
problem. We denote this approximate posterior distribution
$\mathbf{x}|\mathbf{d}\sim\mathcal{N}(\boldsymbol{\mu}_{\mathbf{x}},\boldsymbol{\Gamma}_{\mathbf{x}})$
with the resulting pdf denoted
$\widehat{\pi}_{\mathbf{x}}(\mathbf{x}|\mathbf{d})$.

An important point to note here is that because of the size of the problems
involved, forming the Gauss-Newton Hessian explicitly is not possible. As a
result, explicit computation and storage of $\boldsymbol{\Gamma}_{\mathbf{x}}$
(or an appropriate factorization for it) is not feasible for large-scale
problems since it is likely to be dense; for the same reason, computing the
factor $\mathbf{S}$ is also computationally infeasible. Therefore, to generate
samples we use the preconditioned Lanczos approach outlined in [9]. The key
advantage of this approach is that it only requires matrix-vector products
with $\boldsymbol{\Gamma}_{\mathbf{x}}$ or its inverse.

To explain this approach, consider a factorization of the form
$\boldsymbol{\Gamma}_{\mathbf{x}}=\mathbf{S}\mathbf{S}^{T}$. Note that we will
not compute $\mathbf{S}$ explicitly but use it only for illustration purpose.
Given this factorization, we can generate samples from this distribution using
the formula

$\mathbf{x}=\boldsymbol{\mu}_{\mathbf{x}}+\mathbf{S}\boldsymbol{\eta}\qquad\boldsymbol{\eta}\sim\mathcal{N}(\mathbf{0},\mathbf{I})$
(23)

Suppose we have the preconditioner $\mathbf{G}$ such that
$\mathbf{G}\mathbf{G}^{T}\approx\boldsymbol{\Gamma}_{\mathbf{x}}^{-1}$ and let

$\mathbf{S}=\mathbf{G}^{-T}(\mathbf{G}^{-1}\mathbf{H}_{\text{GN}}\mathbf{G}^{-T})^{-1/2},$

where $\mathbf{H}_{\text{GN}}=\boldsymbol{\Gamma}_{\mathbf{x}}^{-1}$ is the
Gauss-Newton Hessian evaluated at the MAP point. It is easy to verify that

$\mathbf{S}\mathbf{S}^{T}=\mathbf{G}^{-T}(\mathbf{G}^{-1}\mathbf{H}_{\text{GN}}\mathbf{G}^{-T})^{-1/2}(\mathbf{G}^{-1}\mathbf{H}_{\text{GN}}\mathbf{G}^{-T})^{-1/2}\mathbf{G}^{-1}=\boldsymbol{\Gamma}_{\mathbf{x}}.$

With this definition of $\mathbf{S}$, the use of (23) requires the action of
the matrix $(\mathbf{G}^{-1}\mathbf{H}_{\text{GN}}\mathbf{G}^{-T})^{-1/2}$ on
the vector $\boldsymbol{\eta}$. This is accomplished using the Lanczos
algorithm and it only requires matrix-vector products with the matrix
$\mathbf{G}^{-1}\mathbf{H}_{\text{GN}}\mathbf{G}^{-T}$. See [9, Section 2];
the details are omitted here. The goal of the preconditioner here is to reduce
the number of iterations required by the Lanczos algorithm. In practice,
throughout this paper, we choose the preconditioner $\mathbf{G}$ to be the
Cholesky factorization of
$\boldsymbol{\Gamma}_{\text{prior}}^{-1}=\mathbf{G}\mathbf{G}^{T}$. Such a
factorization is easy to compute since
$\boldsymbol{\Gamma}_{\text{prior}}^{-1}$ is block-diagonal with sparse
diagonal blocks.

Let $N_{\text{Lanc}}$ be the number of Lanczos iterations and let
$T_{\text{prec}}$ be the cost of applying a preconditioner. Similar to the
analysis in Section 3.3, the cost of generating a sample is

$T_{\text{sample}}=N_{\text{Lanc}}(2T_{\text{jac}}+T_{\text{prior}}+2T_{\text{prec}})+\mathcal{O}(nN_{\text{LS}}N_{\text{Lanc}})\>\text{flops}.$

Since we use full reorthogonalization in the Lanczos process, there is an
additional cost of $\mathcal{O}((nN_{\text{LS}})N_{\text{Lanc}}^{2})$ flops.

###  4.2 Approximating the posterior variance

A useful way to visualize the uncertainty associated with the reconstruction
is to plot the posterior variance. Since the posterior distribution is non-
Gaussian, computing the conditional variance is not straightforward. Here, we
derive two methods for approximating the posterior variance based on the
Gauss-Newton-Laplace approximation. First, we derive an approximate posterior
distribution for $\mathbf{m}$, conditioned on the data $\mathbf{d}$. Since the
mapping $\mathbf{m}(\mathbf{x})$ is differentiable, using a first-order Taylor
expansion about the posterior mean $\boldsymbol{\mu}_{\mathbf{x}}$

$\mathbf{m}(\mathbf{x})\approx\mathbf{m}(\boldsymbol{\mu}_{\mathbf{x}})+\mathbf{J}_{\mathbf{m}}(\mathbf{x}-\boldsymbol{\mu}_{\mathbf{x}}),$

where $\mathbf{J}_{\mathbf{m}}=\frac{\partial\mathbf{m}}{\partial\mathbf{x}}$
is the Jacobian of $\mathbf{m}$ with respect to $\mathbf{x}$ computed at the
MAP estimate. Let
$\boldsymbol{\mu}_{\mathbf{m}}=\mathbf{m}(\boldsymbol{\mu}_{\mathbf{x}})-\mathbf{J}_{\mathbf{m}}\boldsymbol{\mu}_{\mathbf{x}}$,
then the approximate posterior distribution is

$\mathbf{m}|\mathbf{d}\sim\mathcal{N}(\boldsymbol{\mu}_{\mathbf{m}},\boldsymbol{\Gamma}_{\mathbf{m}})\qquad\boldsymbol{\Gamma}_{\mathbf{m}}=\mathbf{J}_{\mathbf{m}}\boldsymbol{\Gamma}_{\mathbf{x}}\mathbf{J}_{\mathbf{m}}^{\top},$

with the corresponding pdf denoted
$\widehat{\pi}_{\mathbf{m}}(\mathbf{m}|\mathbf{d})$. An approximation of the
posterior variance can, therefore, be obtained from the diagonals of the
matrix $\boldsymbol{\Gamma}_{\mathbf{m}}$. However, forming
$\boldsymbol{\Gamma}_{\mathbf{m}}$ is computationally infeasible for the size
of the problems, since we cannot compute $\boldsymbol{\Gamma}_{\mathbf{x}}$
explicitly. We now propose two different matrix-free methods for estimating
the posterior variance. An advantage of both approaches is that they are
embarrassingly parallel across the random vectors.

####  4.2.1 Method 1: Using LanczosMC diagonal estimators

In the appendix, we propose a new method LanczosMC for estimating the
diagonals of the inverse of a matrix. This method is a hybrid method which
combines the advantages of the Lanczos method for estimating the diagonals and
the Monte Carlo diagonal estimator.

A detailed algorithm, explanation of computational costs, as well as examples
on test matrices will be give in A. Here we only give a high-level overview.
Our approach first approximates $\mathbf{H}_{\text{GN}}^{-1}$ using a low-rank
approximation computed using the preconditioned Lanczos approach so that
$\mathbf{H}_{\text{GN}}^{-1}\approx\mathbf{W}_{k}\mathbf{W}_{k}^{T}$. Recall
that $\boldsymbol{\Gamma}_{\mathbf{x}}=\mathbf{H}_{\text{GN}}^{-1}$, which we
can write as

$\mathbf{H}_{\text{GN}}^{-1}=\mathbf{W}_{k}\mathbf{W}_{k}^{T}+\left(\mathbf{H}_{\text{GN}}^{-1}-\mathbf{W}_{k}\mathbf{W}_{k}^{T}\right).$

The definition of $\boldsymbol{\Gamma}_{\mathbf{m}}$ implies that

$\mathop{\mathrm{diag}}(\boldsymbol{\Gamma}_{\mathbf{m}})=\mathop{\mathrm{diag}}(\mathbf{J}_{\mathbf{m}}\mathbf{W}_{k}\mathbf{W}_{k}^{T}\mathbf{J}_{\mathbf{m}}^{T})+\mathop{\mathrm{diag}}(\boldsymbol{\Gamma}_{\mathbf{m}}-\mathbf{J}_{\mathbf{m}}\mathbf{W}_{k}\mathbf{W}_{k}^{T}\mathbf{J}_{\mathbf{m}}^{T}).$
(24)

The first summand is easy to calculate because it is a low-rank outer product.
The second summand can be estimated with a Monte Carlo diagonal estimator,
similar to the one described in [3]. To estimate
$\mathop{\mathrm{diag}}(\boldsymbol{\Gamma}_{\mathbf{m}}-\mathbf{J}_{\mathbf{m}}\mathbf{W}_{k}\mathbf{W}_{k}^{T}\mathbf{J}_{\mathbf{m}}^{T})$,
we draw $N$ independent random vectors $\\{\mathbf{z}_{j}\\}_{j=1}^{N}$ with
the Rademacher distribution (i.e., entries $\pm 1$ with equal probability). We
then estimate

$\mathop{\mathrm{diag}}(\boldsymbol{\Gamma}_{\mathbf{m}}-\mathbf{J}_{\mathbf{m}}\mathbf{W}_{k}\mathbf{W}_{k}^{T}\mathbf{J}_{\mathbf{m}}^{T})\approx\frac{1}{N}\sum_{j=1}^{N}\mathbf{y}_{j}\odot\mathbf{z}_{j},$

where $\odot$ denotes the elementwise product and

$\mathbf{y}_{j}=\mathbf{J}_{\mathbf{m}}(\mathbf{H}_{\text{GN}}^{-1}(\mathbf{J}_{\mathbf{m}}^{T}\mathbf{z}_{j}))-\mathbf{J}_{\mathbf{m}}\mathbf{W}_{k}\mathbf{W}_{k}^{T}(\mathbf{J}_{\mathbf{m}}^{T}\mathbf{z}_{j}),\qquad
1\leq j\leq N.$

Here to compute
$\mathbf{H}_{\text{GN}}^{-1}(\mathbf{J}_{\mathbf{m}}^{T}\mathbf{z}_{j})$ we
use the CG approach similar to that used in Section 3.2. Note that the only
source of error in (24) is in the second summand because the diagonal of a
low-rank outer product in the first summand can be computed analytically. In
A, we show through a series of numerical experiments that the LanczosMC
approach has lower variance compared to directly applying the Monte Carlo
diagonal estimator to $\boldsymbol{\Gamma}_{\mathbf{m}}$. The cost of solving
a linear system is comparable to that of computing the Gauss-Newton step; we
denote it by $T_{\rm solve}$. The matrix $\mathbf{J}_{\mathbf{m}}$ is very
sparse, so the cost of applying $\mathbf{J}_{\mathbf{m}}$ is negligible
compared to the other costs. Therefore, the cost of using this approach is

$T_{\text{var1}}=k(2T_{\text{jac}}+2T_{\text{prec}}+T_{\text{prior}})+NT_{\rm
solve}+\mathcal{O}(nN_{\text{LS}}(Nk+k^{2}))\>\text{flops}.$

####  4.2.2 Method 2: Using approximate posterior samples

An alternative method to estimate the posterior variance is to use the samples
generated from the approximate posterior distribution with pdf
$\widehat{\pi}_{\mathbf{x}}(\mathbf{x}|\mathbf{d})$. Suppose we have
independent samples $\\{\mathbf{x}_{j}\\}_{j=1}^{N}$ from the distribution
with the pdf $\widehat{\pi}_{\mathbf{x}}(\mathbf{x}|\mathbf{d})$ which can be
computed using the methods proposed in Section 4.1, then define the sample
average
$\widehat{\mathbf{m}}=\frac{1}{N}\sum_{j=1}^{N}\mathbf{m}(\mathbf{x}_{j})$.
Using the multivariate delta method [25, Section 1.8, Theorem 8.22], it can be
shown that

$\sqrt{N}(\widehat{\mathbf{m}}-\mathbf{m}(\boldsymbol{\mu}_{\mathbf{x}}))\longrightarrow\mathcal{N}(\mathbf{0},\mathbf{J}_{\mathbf{m}}\boldsymbol{\Gamma}_{\mathbf{x}}\mathbf{J}_{\mathbf{m}}^{T}),$

as $N\rightarrow\infty$. Here the convergence is in the sense of
distributions. As before, $\mathbf{J}_{\mathbf{m}}$ is the Jacobian computed
at the MAP estimate. Therefore, we can approximate the posterior variance as
the diagonals of the sample covariance operator

$\mathop{\mathrm{diag}}(\boldsymbol{\Gamma}_{\mathbf{m}})\approx\frac{1}{N-1}\mathop{\mathrm{diag}}\left(\sum_{j=1}^{N}(\mathbf{m}(\mathbf{x}_{j})-\widehat{\mathbf{m}})(\mathbf{m}(\mathbf{x}_{j})-\widehat{\mathbf{m}}))^{\top}\right).$

The cost of this approach is simply $T_{\text{var2}}=NT_{\text{sample}}$
flops, where $T_{\text{sample}}$ is the cost of computing a sample using the
approach described in Section 4.1.

##  5 Numerical Experiments

In this section, we demonstrate the Bayesian level set approach and the
performance of our solvers on a set of test problems.

###  5.1 Description of the test problems and parameters

We describe the various test problems and images that we use in this paper.
The first test application (Section 5.2) comes from Photoacoustic Tomography.
The linear forward operator matrix $\mathbf{A}$ is of size $23168\times 16384$
and is generated using the PRspherical function from the IR Tools toolbox
[15]. We use two different test images from IR tools: the first one is called
‘Three Phases’, which has three regions, and the other is called ‘Grains’. The
second application comes from X-ray tomography and uses an open-source
tomographic dataset corresponding to a walnut image [18]. A brief summary of
the problem sizes and the number of level set functions used are summarized
below:

Image | Three Phases | Grains | Walnut
---|---|---|---
Application | Photoacoustic | X-ray
$n_{\rm obs}$ | $23168$ | $23168$ | $9840$
$n$ | $128\times 128$ | $128\times 128$ | $82\times 82$
$N_{\text{LS}}$ | $2$ | $3$ | $2$

##### Level set parameters

In the experiments, we choose the scale parameter $\epsilon=10^{-2}$. We
experimented with different values of $\epsilon$ and observed that with a
smaller $\epsilon$, the images were a lot sharper but consequently the
optimization problem became harder. The number of level sets used depended on
the test image used; in applications, a domain expert may have some prior
knowledge on the maximum number of regions $N_{R}$ possible, and one can use
$N_{\text{LS}}\sim\log_{2}(N_{R})$ using the formula in (4). In the X-ray
dataset, the true image represents a walnut and even though we do not know the
exact number of regions, our approach does a reasonable job estimating the
image.

##### Prior information

In all the experiments, the prior distributions for $\boldsymbol{\Phi}$ is a
Gaussian distribution with zero mean and covariance defined by
$\boldsymbol{\Gamma}_{\boldsymbol{\Phi}}^{-1}=\mathbf{I}_{N_{\text{LS}}}\otimes\lambda_{\boldsymbol{\Phi}}^{2}\left(\alpha\mathbf{L}+\gamma\mathbf{I}\right)^{2}$,
where $\mathbf{L}$ is the finite difference discretization of $(-\Delta)$ with
homogeneous Neumann boundary conditions. The choice of the precision matrix
for each individual level set is related to [4, 2]. The parameters $\alpha$
and $\gamma$ control the correlation length and smoothness of the samples
drawn from the prior and we take $\alpha=0.01$ and $\gamma=0.1$. We also take
the prior information for $\mathbf{c}$ as a Gaussian distribution with zero
mean and covariance matrix
$\boldsymbol{\Gamma}_{\mathbf{c}}=\lambda_{\mathbf{c}}^{2}\mathbf{I}$. We
discuss how to estimate $\lambda_{\boldsymbol{\Phi}}^{2}$ in the experiments
below.

###  5.2 Photoacoustic Tomography

The first set of numerical experiments is based on a test problem from the
`IRTools` package [15] that models photoacoustic tomography (PAT). In all the
experiments in this subsection, we used a regular grid of size $n=128^{2}$.
The matrix $\mathbf{A}$ has the size $23168\times 16384$. We use two different
test images: ‘Three Phases’ and ‘Grains’. The noise level is set to be $2\%$
for both the test images. For the ‘Three Phases’ problem, we take
$N_{\text{LS}}=2$, which is capable of handling $2^{N_{\text{LS}}}=4$ regions.
On the other hand for the ‘Grains’ test problem, which has multiple regions,
we pick $N_{\text{LS}}=3$ which is capable of $8$ distinct regions. For both
problems take $\lambda_{\mathbf{c}}^{2}=10\lambda_{\boldsymbol{\Phi}}^{2}$.
The parameter $\lambda_{\boldsymbol{\Phi}}^{2}$ for each test problem is
determined using an L-curve type approach. In this approach, we solve the
inverse problem for several values of $\lambda_{\boldsymbol{\Phi}}^{2}$ and
for each solution, plot the data misfit
$\left\lVert\mathbf{d}-\mathbf{A}\mathbf{m}(\mathbf{x})\right\rVert_{\boldsymbol{\Gamma}_{\text{noise}}^{-1}}^{2}$
and the regularization term (without the scaling factor
$\lambda_{\boldsymbol{\Phi}}^{2}$). We also plot the estimate of noise, and we
choose $\lambda_{\boldsymbol{\Phi}}^{2}$ at the first instance at which the
solution is smaller than the estimate of the noise. An alternative approach
would be to use the continuation approach described in [17].

Figure 4: These images correspond to the ‘Three Phases’ test problem. In the
top row, we plot the true image (left) and the MAP estimate (right). In the
bottom row, we plot the level set functions $\boldsymbol{\phi}_{1}$ and
$\boldsymbol{\phi}_{2}$ that determine the MAP estimate.

A summary of the relative error and the solver statistics for both test
problems is provided in Table 1. For the ‘Three Phases’ test problems, we
provide the true image and the MAP estimate in Figure 4. In the bottom row of
that same figure, we also provide the level set functions
$\boldsymbol{\phi}_{1}$ and $\boldsymbol{\phi}_{2}$. In Figure 5, we plot
samples from the approximate posterior distribution
$\widehat{\pi}_{\mathbf{m}}(\mathbf{m}|\mathbf{d})$, and on the right panel,
we plot the approximate posterior variance computed using the LanczosMC
approach (Method 1). We use $N=1000$ Monte Carlo samples and the size of the
Lanczos basis was $k=200$. To solve linear systems with $\mathbf{H}_{\rm GN}$,
we used preconditioned GMRES (with no restart) instead of CG, even though the
matrix is symmetric positive definite. This is because CG did not converge
with this preconditioner in a reasonable number of iterations. The number of
GMRES iterations, on average, were $273.3$. For generating the samples and for
the LanczosMC approach, we used the Cholesky factorization of the
$\boldsymbol{\Gamma}_{\text{prior}}^{-1}$ as the preconditioner. The Lanczos
solver for generating the samples took $409$ iterations, on average.

Figure 5: These results correspond to the ‘Three Phases’ test problem: (left)
Samples drawn from the approximate posterior distribution
$\widehat{\pi}_{\mathbf{m}}(\mathbf{m}|\mathbf{d})$, and (right) approximate
posterior variance computed using $N=1000$ Monte Carlo samples.

For the ‘Grains’ test problems, we provide the true image and the MAP estimate
in Figure 6; we also plot the results from the L-curve for the grains example.
In Figure 7, we plot samples from the approximate posterior distribution
$\widehat{\pi}_{\mathbf{m}}(\mathbf{m}|\mathbf{d})$, and on the right panel,
we plot the approximate posterior variance computed using $N=1000$ Monte Carlo
samples and the size of the Lanczos basis was $k=200$; the number of
preconditioned GMRES iteration for solving systems with $\mathbf{H}_{\rm GN}$
was $182.3$. As before, for generating the samples, we used the Cholesky
factorization of $\boldsymbol{\Gamma}_{\text{prior}}^{-1}$ as the
preconditioner and the Lanczos solver took $235.3$ iterations, on average.

Image | $\lambda^{2}_{\boldsymbol{\Phi}}$ | relative error | $N_{\text{LS}}$ | # GN | # CG
---|---|---|---|---|---
Three Phases | $10^{2}/2^{6}$ | $13.9\%$ | $2$ | $8$ | $38$
Grains | $10^{3}/2^{5}$ | $9.78\%$ | $3$ | $18$ | $1847$
Table 1: Summary of the errors and the solver performance corresponding to
‘Three Phases’ and ‘Grains’ test problems.  Figure 6: The image on the left
displays the L-curve corresponding to noise levels $2\%$. The center image
represents the true image and the image on the right displays the MAP
estimate. The relative errors and a summary of the solver performance is given
in Table 1.

Figure 7: These results correspond to the ‘Grains’ test problem corresponding
to $2\%$ noise level: (top) Samples drawn from the Laplace’s approximation to
the posterior distribution
$\widehat{\pi}_{\mathbf{m}}(\mathbf{m}|\mathbf{d})$, and (bottom) approximate
posterior variance.

####  5.2.1 Limited Angle Tomography

For the next test problem, we consider a limited angle PAT problem.
Measurements for the original forward operator $\mathbf{A}$ for PAT consists
of line integrals along concentric circles whose centers span the original
image [15]. For the limited angle case, we cut the number of angles in half,
while maintaining the same angular range. As a consequence the new forward
operator is under-determined and $n_{\rm obs}=11584<n$. The test image we use
is ‘Grains’. We set the noise level at $2\%$, and since the image has multiple
regions we set $N_{\text{LS}}=3$ which can handle up to $2^{N_{\text{LS}}}=8$
regions. As before we take
$\boldsymbol{\Gamma}_{\mathbf{c}}=\lambda^{2}_{\mathbf{c}}\mathbf{I}$ with
$\lambda^{2}_{\mathbf{c}}=10\lambda_{\boldsymbol{\Phi}}^{2}$. In Figure 8 we
plot the true image and the corresponding map estimate.

Figure 8: These images correspond to the ‘Grains’ limited angle test problem:
the true image (left) and the MAP estimate (right).

In Table 2 we summarize the error and solver statistics for the limited angle
‘Grains’ example. In Figure 9 we plot the approximate variance and samples
from the posterior distribution for this example. The Lanczos solver for
approximate samples took $203.8$ iterations, on average, and the number of
GMRES iterations was $166.9$.

noise | $\lambda^{2}_{\boldsymbol{\Phi}}$ | relative error | $N_{\text{LS}}$ | # GN | # CG
---|---|---|---|---|---
$2\%$ | $10^{3}/2^{6}$ | $11.2\%$ | $3$ | $6$ | $35$
Table 2: Summary of the errors and the solver performance corresponding to
limited angle ‘Grains’ test problem.

Figure 9: These results correspond to the ‘Grains’ limited angle test problem
with $2\%$ noise level and the limited angles: (left) Samples drawn from the
approximate posterior distribution
$\widehat{\pi}_{\mathbf{m}}(\mathbf{m}|\mathbf{d})$, and (right) approximate
posterior variance computed using $N=1000$ Monte Carlo samples.

Looking at Tables 1 and 2 we notice that the number of CG iterations for the
‘Grains’ example is much larger than for the ‘Three Phases example. One
explanation for the ‘Grains’ example having several more CG iterations could
be a combination of its larger regularization parameter and the fact that it
has a much smaller relative error than both ’Grains’ limited angle and ‘Three
phases.’ Since the regularization parameter was larger it may have taken many
more optimization iterations to converge.

####  5.2.2 Comparing optimization solvers

To demonstrate the robustness of our inexact Gauss-Newton approach we will
compare it to MATLAB’s Trust Region algorithm and the LBFGS algorithm from the
Poblano toolbox [12]. By default, the LBFGS solver has a limited memory
parameter of $5$. The default stopping criterion for the LBFGS algorithm are:
a maximum of 100 optimization iterations, a maximum of 100 function
evaluations, a scaled gradient norm less than or equal to $10^{-5}$, and a
relative function change less than or equal to $10^{-6}$. See [12] for more
details. As discussed in Subsection 3.2 our solver (inexact Gauss-Newton) uses
both the discrepancy principle and a relative reduction in the gradient as
stopping criterion. For the gradient stopping criterion we require that the
initial gradient reduces by three orders of magnitude. For all of the
optimizers we will set the maximum number of optimization iterations to $50$.
We set the maximum number of CG iterations per optimization iteration to $200$
and set the tolerance for the relative residuals to $10^{-6}$. Note that we
provide the Gauss-Newton Hessian in lieu of the exact Hessian for the Trust
Region solver. We will examine the performance of computing the map estimate
on the ‘Grains’ test problem. For this example we use the same noise and
$\lambda_{\boldsymbol{\Phi}}^{2}$ as reported in Table 1. To compare we report
the number of objective function evaluations, optimization iterations, the
total number of CG iterations and report the relative error. Note that CG
iterations are not reported for the LBFGS solver because the search direction
in the LBFGS algorithm is calculated using a Hessian inverse approximation
formula [12]. The results are given in Table 3.

Solver | #F evals | #Opt iters | #CG iters | Rel error
---|---|---|---|---
LBFGS | 120 | 50(∗) | N/A | 11.4 %
Trust Region | 52 | 50(∗) | 8200 | 9.68 %
Inexact GN | 71 | 18 | 1847 | 9.78%
Table 3: Table comparing performance of different optimization solvers for the
Grains example. The symbol $(*)$ indicates that the optimizer did not converge
within the $50$ iteration limit.

Table 3 demonstrates the benefits of utilizing an inexact Gauss-Newton
approach in favor of other popular methods such as LBFGS. LBFGS completed
$120$ function evaluations as opposed to Trust Region and both inexact Gauss-
Newton methods which were $52$ and $71$ respectively. The Trust Region method
utilized less function evaluations than the inexact Gauss-Newton methods but
had almost four times the number of CG iterations. One explanation for this
may be that since inexact Gauss-Newton only requires CG to be solved within a
tolerance specified by the forcing sequence as in (22). On the other hand,
inexact Gauss-Newton requires more function evaluations at the linesearch
stage.

###  5.3 X-ray tomography

In this set of numerical experiments, we use the real-data from [18]. The
dataset consists of the logarithm of the x-ray sinogram of a single two-
dimensional slice of a walnut and the forward operator $\mathbf{A}$ modeling
the x-ray transform for this particular set. We use all 120 projections that
are provided in the dataset; therefore, the size of the matrix $\mathbf{A}$ is
$9840\times 6724$. From the ground truth image, while the image is not exactly
piecewise constant, we use $N_{\text{LS}}=2$ which yields
$2^{N_{\text{LS}}}=4$ regions. The parameter $\lambda_{\boldsymbol{\Phi}}^{2}$
is chosen using an L-curve type analysis; although we do not know the exact
noise level we use the curvature of the L-curve to guide the selection of the
regularization parameter. We take the noise covariance matrix to be
$\boldsymbol{\Gamma}_{\mathrm{noise}}=\sigma^{2}\mathbf{I}$ where
$\sigma=0.035$ and $\lambda_{\boldsymbol{\Phi}}=10^{3}/2^{5}$. We also choose
$\lambda_{\mathbf{c}}^{2}=50\lambda_{\boldsymbol{\Phi}}^{2}$.

Figure 10: (left) Ground truth obtained using the filtered backprojection
technique. The image is provided at a much higher resolution and is meant for
visual comparison only. (right) MAP estimate obtained using $N_{\text{LS}}=2$
level set functions.

To compute the MAP estimate, we use the inexact Gauss-Newton solver which
converged in $18$ iterations and required $683$ CG iterations overall. Since
prior knowledge of the magnitude of the noise is not available, the
optimization routine uses a stopping criterion based on the reduction of the
relative magnitude by four orders of magnitude. The MAP estimate has been
displayed alongside the ground truth image obtained using a filtered
backprojection technique; note that this image is at a much higher resolution
than what is being used to solve the inverse problem. From a visual
perspective, it is clear that the MAP estimate captures the main outline of
the walnut shape and we do not provide a quantitative comparison for the
reasons described above. In addition to the MAP estimate, we provide samples
from the approximate posterior
$\widehat{\pi}_{\mathbf{m}}(\mathbf{m}|\mathbf{d})$ in the left panel of
Figure 11. Finally, we provide the approximate posterior variance, computed
using $N=1000$ Monte Carlo samples and Lanczos basis $k=100$, in the right
panel of the same figure. The number of GMRES iterations, on average, was
$195.8$. The Lanczos solver for the approximate samples took $248.7$
iterations, on average. The image shows a higher variance around the edges of
the walnut and at the ends of the images.

Figure 11: These results correspond to the ‘Walnut’ test problem: (left)
Samples drawn from the approximate posterior distribution
$\widehat{\pi}_{\mathbf{m}}(\mathbf{m}|\mathbf{d})$, and (right) approximate
posterior variance.

###  5.4 Hydraulic Tomography

In this example, we present a nonlinear hydraulic tomography problem adapted
from [24]. A synthetic vertical confined aquifer consisting of two geologic
facies is considered and the subsurface flow system is governed by

$-\nabla\cdot\left(K\nabla
h_{i}\right)=q_{i}\delta\left(\mathbf{s}_{i}\right)\qquad 1\leq i\leq n_{s}$
(25)

where $K$ is the hydraulic conductivity [$m/s$], $h$ is the hydraulic head
[$m$], $q_{i}$ is the 2-D pumping rate for a pumping test at a well location
$\mathbf{s}_{i}$ (marked with a circle in Figure 12) for $1\leq i\leq n_{s}$,
$n_{s}$ is the number of sources, and $\delta(\mathbf{x})$ is the Dirac delta
function. A constant-head boundary condition on the left, right, and upper
boundaries as well as no-flux condition on the bottom is imposed. The
measurements are steady state head changes due to sequential pump tests at
each well location. For each test water is extracted from one of the $20$ well
locations at the rate of $2$ L/s and the change in the head is measured at all
other 19 locations resulting in $380$ drawdown measurements ($20$ pump tests
$\times$ $19$ measurements per test). To ensure positive hydraulic
conductivity values during the inversion, log-transformed hydraulic
conductivity $\log K$ is estimated and 1% Gaussian noise ($\approx$ std(error)
= $0.12$ m) was added to the measurements. The domain was discretized into
$100$ by $100$ grids so that $n=10^{4}$. The governing equation in (25) is
solved by a finite volume method.

Figure 12: Ground truth log hydraulic conductivity field ($\log_{10}K$) for a
hydraulic tomography example consisting of gravel ($\log_{10}K=-5$ m/s) and
sand ($\log_{10}K=-3$ m/s). White circles are the pumping/observation well
locations.

To compute the MAP estimate, the inexact Gauss-Newton solver is used and it
converged in $17$ iterations with $1006$ CG iterations in total. Because of
log-transformation, the prior covariance for $\boldsymbol{\Phi}$ is
$\boldsymbol{\Gamma}_{\boldsymbol{\Phi}}^{-1}=\mathbf{I}_{N_{\text{LS}}}\otimes\lambda_{\boldsymbol{\Phi}}^{2}\left(\alpha\mathbf{L}+\gamma\mathbf{I}\right)^{2}$
and for $\mathbf{c}$ as a Gaussian distribution with mean $\log_{10}K=-4$ and
covariance matrix
$\boldsymbol{\Gamma}_{\mathbf{c}}=\lambda_{\mathbf{c}}^{2}\mathbf{I}$ and
$\lambda_{\mathbf{c}}^{2}=\lambda_{\boldsymbol{\Phi}}^{2}$. The MAP estimate
is displayed along with the true $\log_{10}K$ field in Figure 12. The estimate
identifies the main circle structures with different values of $K$. Similar to
the previous examples we plot samples from the approximate posterior
distribution as well as the approximate posterior variance (using LanczosMC
with the same settings as before) in Figure 13. The Lanczos solver took $44.2$
iterations, on average. To compute the approximate posterior variance, we used
preconditioned CG instead of GMRES in contrast to the previous applications;
the total number of iterations, on average, was $207.9$.

Figure 13: These results correspond to the hydraulic tomography test problem:
(left) Samples from the approximate posterior distribution
$\widehat{\pi}_{\mathbf{m}}(\mathbf{m}|\mathbf{d})$ and (right) approximate
posterior variance. White circles are the pumping/observation well locations.

In Figure 13 we see that the magnitude of the uncertainty is high at the
interfaces of the regions especially in the areas far away from the well
locations as expected. The high uncertainty regions are well matched with the
interfaces of the geologic structures. In the interior of largest high
permeable circle there is a decrease in uncertainty. The largest circle
contains more measurements than the other ones which would imply that the
uncertainty would be smaller there.

##  6 Conclusions

In this paper we have developed a Bayesian level set approach to inverse
problems where the unknown is a piecewise constant field. This was
accomplished by assigning prior distributions to level set functions which
define different piecewise constant regions and the constants that define
their magnitudes. Assuming an additive Gaussian noise model and Gaussian prior
we derived the (potentially) non-Gaussian posterior distribution and developed
an inexact Gauss-Newton approach to compute the MAP estimate that leveraged
efficient adjoint-based methods for computing derivatives and sparsity in the
Jacobian, of the piecewise constant model. We employed Laplace’s approximation
to approximate the posterior and developed matrix-free methods for sampling
from the posterior and estimating the posterior variance. In addition we
developed a matrix-free estimator for the diagonals of the inverse of a
matrix. Our algorithms were performed on computational examples in
photoacoustic, hydraulic, and X-ray tomography. In each example we obtained an
accurate reconstruction and gained qualitative information about the
uncertainty in the reconstruction using samples from the approximate posterior
and the approximate posterior variance. It is important to emphasize that the
approximated uncertainty estimates are to be interpreted qualitatively.
Although we have proposed efficient and matrix-free methods to Bayesian
inverse problems with piecewise constant unknowns there are a number of
possible future directions for our work.

An immediate extension of this work would be to implement a Gaussian prior
with the Matérn covariance. The Matérn covariance is a stationary covariance
kernel that allows for control of the smoothness and variance of prior samples
[16]. As a result, the Matérn covariance provides more flexibility than the
inverse elliptical differential operator we have implemented in this work.
This flexibility comes at the cost of sparsity in the prior precision matrix
$\boldsymbol{\Gamma}^{-1}_{\text{prior}}$ which we have leveraged in our
algorithms. The new prior covariance will be a large dense matrix and
explicitly forming or factorizing the precision matrix will be computationally
infeasible and introduces new challenges in our approaches for MAP estimation
and uncertainty quantification. To tackle these challenges we would build on
the methods in [33]. A topic that we have not addressed in this paper is the
validation of linearized Bayesian inference. One way to accomplish this would
be to compare samples generated from the method in 4.1 to samples generated
from a Gauss-Newton Hessian based MCMC sampler [30].

##  7 Acknowledgements

This work was supported, in part, by the National Science Foundation through
the award DMS-1845406 (W.R. and A.K.S.) and DGE-1746939 (W.R.).

##  8 References

## References

  * [1] A. Aghasi, M. Kilmer, and E. L. Miller.  Parametric level set methods for inverse problems.  SIAM Journal on Imaging Sciences, 4(2):618–650, 2011.
  * [2] J. M. Bardsley.  Computational uncertainty quantification for inverse problems, volume 19 of Computational Science & Engineering.  Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2018.
  * [3] C. Bekas, E. Kokiopoulou, and Y. Saad.  An estimator for the diagonal of a matrix.  Applied numerical mathematics, 57(11-12):1214–1229, 2007.
  * [4] T. Bui-Thanh and O. Ghattas.  A scalable algorithm for MAP estimators in Bayesian inverse problems with Besov priors.  Inverse Problems & Imaging, 9(1):27, 2015.
  * [5] T. Bui-Thanh, O. Ghattas, J. Martin, and G. Stadler.  A computational framework for infinite-dimensional Bayesian inverse problems part I: The linearized case, with application to global seismic inversion.  SIAM Journal on Scientific Computing, 35(6):A2494–A2523, 2013.
  * [6] M. Burger and S. J. Osher.  A survey on level set methods for inverse problems and optimal design.  European journal of applied mathematics, 16(2):263–301, 2005.
  * [7] M. Cardiff and P. Kitanidis.  Bayesian inversion for facies detection: An extensible level set framework.  Water Resources Research, 45, 2009.
  * [8] G. Chantas, N. Galatsanos, A. Likas, and M. Saunders.  Variational Bayesian image restoration based on a product of t-distributions image prior.  IEEE Transactions on Image Processing, 17(10):1795–1805, 2008.
  * [9] E. Chow and Y. Saad.  Preconditioned Krylov subspace methods for sampling multivariate Gaussian distributions.  SIAM Journal on Scientific Computing, 36(2):A588–A608, 2014.
  * [10] T. A. Davis and Y. Hu.  The university of Florida sparse matrix collection.  ACM Transactions on Mathematical Software, 38(1):1–25, 2011.
  * [11] O. Dorn and D. Lesselier.  Level set methods for inverse scattering.  Inverse Problems, 22(4):R67, 2006.
  * [12] D. M. Dunlavy, T. G. Kolda, and E. Acar.  Poblano v1.0: A matlab toolbox for gradient-based optimization.  Number SAND2010-1422, March 2010.
  * [13] M. M. Dunlop, M. A. Iglesias, and A. M. Stuart.  Hierarchical Bayesian level set inversion.  Statistics and Computing, 27(6):1555–1584, 2017.
  * [14] M. M. Dunlop and A. M. Stuart.  MAP estimators for piecewise continuous inversion.  Inverse Problems, 32(10):105003, 2016.
  * [15] S. Gazzola, P. C. Hansen, and J. G. Nagy.  IR tools: a MATLAB package of iterative regularization methods and large-scale test problems.  Numerical Algorithms, 81(3):773–811, 2019.
  * [16] M. G. Genton.  Classes of kernels for machine learning: A statistics perspective.  J. Mach. Learn. Res., 2:299–312, mar 2002.
  * [17] E. Haber, U. M. Ascher, and D. Oldenburg.  On optimization techniques for solving nonlinear inverse problems.  Inverse problems, 16(5):1263, 2000.
  * [18] K. Hämäläinen, L. Harhanen, A. Kallonen, A. Kujanpää, E. Niemi, and S. Siltanen.  Tomographic x-ray data of a walnut.  arXiv preprint arXiv:1502.04064, 2015.
  * [19] P. C. Hansen.  Discrete inverse problems, volume 7 of Fundamentals of Algorithms.  Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2010.  Insight and algorithms.
  * [20] M. A. Iglesias, Y. Lu, and A. M. Stuart.  A Bayesian level set method for geometric inverse problems.  Interfaces and Free Boundaries, 18, 2016.
  * [21] J. Kaipio and E. Somersalo.  Statistical and computational inverse problems, volume 160 of Applied Mathematical Sciences.  Springer-Verlag, New York, 2005.
  * [22] M. Lassas and S. Siltanen.  Can one use total variation prior for edge preserving Bayesian inversion?  Inverse Problems, 20(5):1537–1563, 2004.
  * [23] S. L. Lassi Roininen, Janne M. J. Huttunen.  Whittle-Matérn priors for Bayesian statistical inversion with applications in electrical impedance tomography.  Inverse Problems & Imaging, 8(2):561–586, 2014.
  * [24] J. Lee and P. Kitandis.  Bayesian inversion with total variation prior for discrete geologic structure identification.  Water Resources Research, 49:7658–7669, 2013.
  * [25] E. L. Lehmann and G. Casella.  Theory of point estimation.  Springer Science & Business Media, 2006.
  * [26] F. Lindgren, H. Rue, and J. Lindström.  An explicit link between Gaussian fields and Gaussian Markov random fields: the stochastic partial differential equation approach.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 73(4):423–498, 2011.
  * [27] M. Markkanen, L. Roininen, J. M. Huttunen, and S. Lasanen.  Cauchy difference priors for edge-preserving Bayesian inversion.  Journal of Inverse and Ill-posed Problems, 27(2):225–240, 2019\.
  * [28] R. A. Meyer, C. Musco, C. Musco, and D. P. Woodruff.  Hutch++: Optimal stochastic trace estimation.  In Symposium on Simplicity in Algorithms (SOSA), pages 142–155. SIAM, 2021.
  * [29] J. Nocedal and S. J. Wright.  Numerical Optimization.  Springer, 1 edition, 1999.
  * [30] N. Petra, J. Martin, G. Stadler, and O. Ghattas.  A computational framework for infinite-dimensional Bayesian inverse problems, part ii: Stochastic Newton MCMC with application to ice sheet flow inverse problems.  SIAM Journal on Scientific Computing, 36(4):A1525–A1555, 2014.
  * [31] C. E. Rasmussen.  Gaussian processes in machine learning.  In Summer school on machine learning, pages 63–71. Springer, 2003\.
  * [32] H. Rue and L. Held.  Gaussian Markov Random Fields: Theory and Applications, volume 104 of Monographs on Statistics and Applied Probability.  Chapman & Hall, London, 2005.
  * [33] A. Saibaba, J. Chung, and K. Petroske.  Efficient Krylov subspace methods for uncertainty quantification in large Bayesian linear inverse problems.  Numerical Linear Algebra with Applications, 27, 08 2020.
  * [34] A. M. Stuart.  Inverse problems: a Bayesian perspective.  Acta Numerica, 19:451–559, 2010.
  * [35] K. van den Doel, U. M. Ascher, and A. Leitao.  Multiple level sets for piecewise constant surface reconstruction in highly ill-posed problems.  Journal of Scientific Computing, 43(1):44–66, 2010.
  * [36] U. Villa, N. Petra, and O. Ghattas.  hIPPYlib: An Extensible Software Framework for Large-Scale Inverse Problems Governed by PDEs: Part I: Deterministic Inversion and Linearized Bayesian Inference.  ACM Transactions on Mathematical Software (TOMS), 47(2):1–34, 2021\.
  * [37] C. R. Vogel.  Computational methods for inverse problems, volume 23 of Frontiers in Applied Mathematics.  Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2002.  With a foreword by H. T. Banks.
  * [38] L. Wu, J. Laeuchli, V. Kalantzis, A. Stathopoulos, and E. Gallopoulos.  Estimating the trace of the matrix inverse by interpolating from the diagonal of an approximate inverse.  Journal of Computational Physics, 326:828–844, 2016.

##  Appendix A Estimating the Diagonal of $\mathbf{A}^{-1}$

In this appendix, we give the details of the new estimator (LanczosMC) for the
diagonals of the inverse of a positive definite matrix. This algorithm is
matrix-free and combines two existing diagonal estimators: based on the the
Lanczos approach [8] and the Monte Carlo approach [3].

###  A.1 Lanczos Diagonal Estimator

Let $\mathbf{A}\in\mathbb{R}^{n\times n}$ be a symmetric positive definite
(SPD) matrix. To estimate the diagonal of $\mathbf{A}^{-1}$ we make use of a
method introduced by [8, Section IV], which uses the Lanczos iteration. We
present this in a slightly different form using the preconditioned Lanczos
iteration. Suppose we have a preconditioner $\mathbf{G}\in\mathbb{R}^{n\times
n}$ such that $\mathbf{A}\approx\mathbf{G}\mathbf{G}^{T}$, $\mathbf{G}$ is
easy to invert, and the preconditioned matrix
$\mathbf{G}^{-1}\mathbf{A}\mathbf{G}^{-T}$ has a smaller condition number than
$\mathbf{A}$. Then after $k<n$ steps of the Lanczos iteration on the
preconditioned matrix, assuming no breakdown, we have

$\mathbf{G}^{-1}\mathbf{A}\mathbf{G}^{-T}\mathbf{V}_{k}=\mathbf{V}_{k}\mathbf{T}_{k}+\beta_{k+1}\mathbf{v}_{k+1}\mathbf{e}_{k}^{T},$

where $\mathbf{V}_{k}\in\mathbb{R}^{n\times k}$ has orthonormal columns and
$\mathbf{T}_{k}\in\mathbb{R}^{k\times k}$ is a symmetric, tridiagonal matrix.
Since $\mathbf{A}$ is SPD, $\mathbf{T}_{k}$ is also SPD due to Cauchy
interlacing theorem; therefore, we can compute the Cholesky factorization
$\mathbf{T}_{k}=\mathbf{L}_{k}\mathbf{L}_{k}^{T}$. Let us define
$\mathbf{W}_{k}=\mathbf{G}^{-T}\mathbf{V}_{k}\mathbf{L}_{k}^{-T}$. After $n$
iterations of the Lanczos process, we have
$\mathbf{V}_{n}^{T}\mathbf{G}^{-1}\mathbf{A}\mathbf{G}^{-T}\mathbf{V}_{n}=\mathbf{T}_{n}$,
so that

$\mathbf{A}^{-1}=\mathbf{G}^{-T}\mathbf{V}_{n}\mathbf{T}_{n}^{-1}\mathbf{V}_{n}^{T}\mathbf{G}^{-1}=\mathbf{W}_{n}\mathbf{W}_{n}^{T}.$

Suppose we define $\mathbf{B}_{k}=\mathbf{W}_{k}\mathbf{W}_{k}^{T}$, then we
can approximate $\mathbf{A}^{-1}\approx\mathbf{B}_{k}$. We can estimate the
diagonals of $\mathbf{A}^{-1}$ using the relation

$\mathbf{e}_{i}^{T}\mathbf{A}^{-1}\mathbf{e}_{i}\approx\mathbf{e}_{i}^{T}\mathbf{B}_{k}\mathbf{e}_{i}=\mathbf{e}_{i}^{T}\mathbf{W}_{k}\mathbf{W}_{k}^{T}\mathbf{e}_{i}=\sum_{j=1}^{k}w_{ij}^{2}.$

Note that this is a monotonically increasing estimate for the $i-$th diagonal
since $\mathbf{B}_{k}=\mathbf{B}_{k-1}+\mathbf{w}_{k}\mathbf{w}_{k}^{T}$. Note
that we can easily adapt this algorithm to estimate the diagonals of
$\mathbf{J}\mathbf{A}^{-1}\mathbf{J}^{T}$ by instead working with
$\mathbf{J}\mathbf{B}_{k}\mathbf{J}^{T}$. In practice, the vectors which form
the columns of $\mathbf{V}_{k}$ lose orthogonality in finite precision, so we
reorthogonalize them against previous vectors. A potential downside of this
approach is that a large number of iterations $k$ may be required to estimate
the diagonals accurately, even with the use of a preconditioner. Although this
estimator does not require storing $\mathbf{W}_{k}$, our approach in A.3
requires the storage of $\mathbf{W}_{k}$ and this motivates the need to use a
small $k$.

###  A.2 Monte Carlo Diagonal Estimator

Another method to estimate $\mathop{\mathrm{diag}}(\mathbf{A}^{-1})$ is using
the matrix-free Monte Carlo diagonal estimator proposed in [3]. We describe
this method to estimate the diagonals of $\mathbf{A}^{-1}$, but this approach
is applicable to any matrix $\mathbf{B}$. In this approach, we draw a set of
random vectors $\mathbf{z}_{\ell}$ drawn from the Rademacher distribution,
i.e., each vector has entries $\pm 1$ with equal probability. Then, we have
the following approximation

$\mathbf{d}^{N}_{\text{MC}}=\frac{1}{N}\left(\sum_{\ell=1}^{N}\mathbf{y}_{\ell}\odot\mathbf{z}_{\ell}\right)\qquad\mathbf{y}_{\ell}=\mathbf{A}^{-1}\mathbf{z}_{\ell},\quad
1\leq\ell\leq N,$ (26)

where $\odot$ denotes elementwise multiplication. The vectors
$\mathbf{y}_{\ell}$ are computed in a matrix-free fashion with a CG solver
applied to the system $\mathbf{A}\mathbf{y}_{\ell}=\mathbf{z}_{\ell}$ for
$1\leq\ell\leq N$; solving these linear systems dominates the computational
cost of the MC estimator. It can be shown that this estimate is unbiased in
the sense that
$\mathbb{E}[\mathbf{d}^{N}_{\text{MC}}]=\mathop{\mathrm{diag}}(\mathbf{A}^{-1})$,
where $\mathbb{E}$ denotes the expectation. However, the convergence of this
estimator is slow in the sense that a large number of vectors $N$ are required
for an accurate estimate of the diagonals of $\mathbf{A}^{-1}$. The authors in
[3] used a probing technique for estimating the diagonals of a matrix, that
takes into account special properties of the matrix such as bandedness and
decay of the magnitudes of the off-diagonal elements of the matrix.

###  A.3 LanczosMC Estimator

Our approach for approximating $\mathop{\mathrm{diag}}(\mathbf{A}^{-1})$
combines the ideas from Lanczos and Monte Carlo diagonal estimator. First, we
approximate $\mathbf{A}^{-1}\approx\mathbf{W}_{k}\mathbf{W}_{k}^{T}$ using the
Lanczos method as described in A.1. Then we decompose $\mathbf{A}^{-1}$ as

$\mathbf{A}^{-1}=\mathbf{W}_{k}\mathbf{W}_{k}^{T}+\left(\mathbf{A}^{-1}-\mathbf{W}_{k}\mathbf{W}_{k}^{T}\right).$
(27)

This equation implies that
$\mathop{\mathrm{diag}}(\mathbf{A}^{-1})=\mathop{\mathrm{diag}}(\mathbf{W}_{k}\mathbf{W}_{k}^{T})+\mathop{\mathrm{diag}}\left(\mathbf{A}^{-1}-\mathbf{W}_{k}\mathbf{W}_{k}^{T}\right)$.
Once the preconditioned Lanczos vectors are computed, the first summand can be
calculated efficiently as it is a low-rank outer product. The second summand
can be estimated via the MC method applied to
$\mathbf{A}^{-1}-\mathbf{W}_{k}\mathbf{W}_{k}^{T}$ rather than
$\mathbf{A}^{-1}$, as described in A.2. The rationale behind this approach is
that when $k$ is sufficiently large, we expect the approximation
$\mathbf{A}^{-1}\approx\mathbf{W}_{k}\mathbf{W}_{k}^{T}$ to be reasonably
accurate. However, numerical evidence suggests that the variance of the
diagonal estimator applied to
$\mathbf{A}^{-1}-\mathbf{W}_{k}\mathbf{W}_{k}^{T}$ is much smaller than
$\mathbf{A}^{-1}$ resulting in more accurate estimators. Therefore, the use of
Lanczos diagonal estimator can be interpreted as a variance reduction
technique. This idea has been explored by [38, 28] in the context of trace
estimators. At this point, we lack a complete theoretical understanding of
this but numerical experiments in A.4 provide convincing numerical evidence in
support of our observation. We summarize this new method for estimating
$\mathop{\mathrm{diag}}(\mathbf{A}^{-1})$ in Algorithm 1.

Related to the diagonal estimator, we can also estimate the trace of the
matrix $\mathbf{A}^{-1}$ as

$\text{trace}(\mathbf{A}^{-1})\approx\text{trace}(\mathbf{W}_{k}^{T}\mathbf{W}_{k})+\frac{1}{N}\sum_{\ell=1}^{N}\mathbf{y}_{\ell}^{T}\mathbf{z}_{\ell},$

where we have used the relation (27) and the cyclic property of trace.

0: SPD matrix $\mathbf{A}$, preconditioner $\mathbf{G}$, integers $k$ denoting
the number of iterations of the Lanczos approach and $N$ denoting the number
of Monte Carlo samples.

1: return An estimator $\mathbf{d}^{k,N}_{\text{LanczosMC}}$ for
$\mathop{\mathrm{diag}}(\mathbf{A}^{-1})$.

2: Run $k$ steps of the Lanczos approach on
$\mathbf{G}^{-1}\mathbf{A}\mathbf{G}^{-T}$ to compute $\mathbf{V}_{k}$ and
$\mathbf{T}_{k}$.

3: Compute $\mathbf{W}_{k}=\mathbf{G}^{-T}\mathbf{B}_{k}\mathbf{L}_{k}^{-T}$
and
$\mathbf{d}^{k}_{\text{Lanczos}}=\mathop{\mathrm{diag}}(\mathbf{W}_{k}\mathbf{W}_{k}^{T}).$

4: Draw independent vectors $\\{\mathbf{z}_{\ell}\\}_{\ell=1}^{N}$ from the
Radamacher distribution and compute
$\mathbf{y}_{\ell}=\left(\mathbf{A}^{-1}-\mathbf{W}_{k}\mathbf{W}_{k}^{T}\right)\mathbf{z}_{\ell}$
for $1\leq\ell\leq N$.

5: Compute the Monte Carlo estimate
$\mathbf{d}^{N}_{\text{MC}}=\frac{1}{N}\left(\sum_{\ell=1}^{N}\mathbf{y}_{\ell}\odot\mathbf{z}_{\ell}\right)$
where $\odot$ denotes elementwise multiplication.

6:
$\mathbf{d}^{k,N}_{\text{LanczosMC}}=\mathbf{d}^{k}_{\text{Lanczos}}+\mathbf{d}^{N}_{\text{MC}}$

Algorithm 1 LanczosMC approach for estimating
$\mathop{\mathrm{diag}}(\mathbf{A}^{-1})$

##### Computational Costs

To analyze the costs of Algorithm 1, we make the following assumptions. We
assume that the cost of applying the preconditioner $\mathbf{G}^{-1}$ and
$\mathbf{G}^{-T}$ is the same and is denoted by $T_{\text{prec}}$.
Furthermore, we denote the cost of applying $\mathbf{A}$ and $\mathbf{A}^{-1}$
to a vector as $T_{\text{apply}}$ and $T_{\text{solve}}$ respectively. If an
iterative method, such as CG, is used to apply $\mathbf{A}^{-1}$, then the
computational cost is
$T_{\text{solve}}=N_{\text{iter}}(T_{\text{apply}}+2T_{\text{prec}})+\mathcal{O}(nN_{\text{iter}})$
flops; here $N_{\text{iter}}$ is the number of iterations taken by CG which
may be different than $k$.

To calculate $\mathbf{d}^{k}_{\text{Lanczos}}$, we first perform $k$ steps of
Lanczos on the preconditioned matrix
$\mathbf{G}^{-1}\mathbf{A}\mathbf{G}^{-T}$. The cost of running $k$ steps of
this algorithm to generate $\mathbf{V}_{k}$ and $\mathbf{T}_{k}$ is
$k(T_{\text{apply}}+2T_{\text{prec}})+\mathcal{O}(nk)$ flops; the additional
cost of computing $\mathbf{W}_{k}$ and $\mathbf{d}^{k}_{\text{Lanczos}}$ is
$\mathcal{O}(nk^{2})$ flops. Therefore, the total cost of computing
$\mathbf{d}^{k}_{\text{Lanczos}}$ is

$T_{\text{Lanczos}}=k(T_{\text{apply}}+2T_{\text{prec}})+\mathcal{O}(nk^{2})\text{
flops}.$

Each application of $\mathbf{A}^{-1}-\mathbf{W}_{k}\mathbf{W}_{k}^{T}$ costs
$T_{\text{solve}}+\mathcal{O}(nk^{2})$, so the overall cost of computing
$\mathbf{d}^{N}_{\text{MC}}$ is

$T_{\text{MC}}=N(T_{\text{solve}}+\mathcal{O}(nk))\text{ flops}.$

The total cost of the LanczosMC approach is

$\displaystyle T_{\text{MC}}=$
$\displaystyle\>T_{\text{Lanczos}}+T_{\text{MC}}$ $\displaystyle=$
$\displaystyle\>k(T_{\text{apply}}+2T_{\text{prec}})+NT_{\text{solve}}+\mathcal{O}(nNk+nk^{2})\text{
flops}.$

In addition to the computational cost discussed here, the method requires
storing $k$ vectors of length $n$. We briefly comment on the choice of
parameters $k$ and $N$. The convergence of Monte Carlo estimators is slow,
like $N^{-1/2}$, so we take $N$ to be large say $100-1000$. However, this
computational cost can be trivially parallelized. The choice of $k$ depends on
the storage budget and accuracy. In numerical experiments, we found that
$k\sim 100$ was a reasonable choice.

###  A.4 Computational Examples

To demonstrate the efficiency of A.3 we estimate the diagonal of the inverse
of 5 different test matrices. The test matrices were taken from [10] and were
chosen to represent a wide variety of sizes and condition numbers. Each matrix
is sparse and SPD. We chose the matrices to be of relatively small size, so
that we can study the accuracy of the estimators; however, our algorithms are
applicable to larger problem sizes. In Table 4 we list the name, size, and
condition number for each test matrix.

##### Experiment 1: Effect of varying $N$

We apply Algorithm 1 to the test matrices in Table 4. We chose $k=100$ and
used the incomplete Cholesky preconditioner. We used CG to solve systems with
$\mathbf{A}$. We compute the relative error

$\text{rel
err}=\frac{\|\mathbf{d}_{\text{C}}-\mathbf{d}_{\text{True}}\|_{2}}{\|\mathbf{d}_{\text{True}}\|_{2}}$

where $\mathbf{d}_{\text{C}}$ is the approximate diagonal and
$\mathbf{d}_{\text{True}}$ is the true diagonal. The results are reported in
Figure 14. The errors are averaged over $100$ different runs for a fixed
sample size. In addition to plotting the average we also plot the $2.5$th and
$97.5$th percentiles for the Monte Carlo and LanczosMC errors as well. From
the plots it is clear that the LanczosMC estimator is accurate compared to MC
estimator. It is also seen that the accuracy of the estimators improves as the
condition number of the matrix increases. We see improvements in relative
error by $1-2$ orders of magnitude (except the matrix mesh3em5, which is well-
conditioned), suggesting that the Lanczos estimator is better for more ill-
conditioned systems.

Matrix | mesh3em5 | nos3 | Trefethen500 | 1138bus | mhdb416
---|---|---|---|---|---
Size | $289\times 289$ | $960\times 960$ | $500\times 500$ | $1138\times 1138$ | $416\times 416$
Condition Number | $5$ | $4.63\times 10^{3}$ | $7.35\times 10^{4}$ | $1.23\times 10^{7}$ | $5.052\times 10^{9}$
Table 4: Summary of the test matrices used in the numerical experiments.

Figure 14: Relative error plots for the Monte Carlo, LanczosMC, and Lanczos
estimators. The error statistics were calculated for 100 realizations for
every fixed value of $N$.

In Figure 14 the shaded area between the $2.5$th and $97.5$th quantiles for
the MC and Lanczos estimators corresponds to the 95% confidence interval. For
each test matrix we notice that the confidence interval for the LanczosMC
estimator is smaller than the MC estimator for every test matrix. Since the
confidence interval is a measure of variance in an estimator, this implies
that the LanczosMC estimator has a smaller variance.

##### Experiment 2: Effect of varying $k$

For this next experiment we analyze the impact of changing the number of
Lanczos iterations $k$. We choose $\mathbf{A}$ to be the mhdb416 matrix and
plot the mean for the error while using the LanczosMC estimator. We also plot
the mean of the error while using the Lanczos estimator. We choose the same
preconditioner, range of Monte Carlo samples, and number of trials for each
fixed number of Monte Carlo samples. The value of $k$ ranges from 20 to 200
and the results are reported in Figure 15.

Figure 15: Relative error plots for different number of Lanczos iterations
$k$: (left) mean of relative errors corresponding to Lanczos estimator, and
(right) Mean of relative errors corresponding to LanczosMC estimator.