# Conditional score-based diffusion models for solving inverse problems in
mechanics

Agnimitra Dasgupta  Harisankar Ramaswamy  Javier Murgoitio Esandi  Ken Foo
Runze Li  Qifa Zhou  Brendan F Kennedy  Assad A Oberai<EMAIL_ADDRESS>
###### Abstract

We propose a framework to perform Bayesian inference using conditional score-
based diffusion models to solve a class of inverse problems in mechanics
involving the inference of a specimen’s spatially varying material properties
from noisy measurements of its mechanical response to loading. Conditional
score-based diffusion models are generative models that learn to approximate
the score function of a conditional distribution using samples from the joint
distribution. More specifically, the score functions corresponding to multiple
realizations of the measurement are approximated using a single neural
network, the so-called score network, which is subsequently used to sample the
posterior distribution using an appropriate Markov chain Monte Carlo scheme
based on Langevin dynamics. Training the score network only requires
simulating the forward model. Hence, the proposed approach can accommodate
black-box forward models and complex measurement noise. Moreover, once the
score network has been trained, it can be re-used to solve the inverse problem
for different realizations of the measurements. We demonstrate the efficacy of
the proposed approach on a suite of high-dimensional inverse problems in
mechanics that involve inferring heterogeneous material properties from noisy
measurements. Some examples we consider involve synthetic data, while others
include data collected from actual elastography experiments. Further, our
applications demonstrate that the proposed approach can handle different
measurement modalities, complex patterns in the inferred quantities, non-
Gaussian and non-additive noise models, and nonlinear black-box forward
models. The results show that the proposed framework can solve large-scale
physics-based inverse problems efficiently.

###### keywords:

Conditional generative models, inverse problems, Bayesian inference,
diffusion-based modeling, uncertainty quantification, elastography

[label1]organization=Department of Aerospace & Mechanical Engineering,
University of Southern California, city=Los Angeles, postcode=90089,
state=California, country=USA

[label2]organization=BRITElab, Harry Perkins Institute of Medical Research,
QEII Medical Centre, city=Nedlands, postcode=6009, state=Western Australia,
country=Australia

[label3]organization=Centre for Medical Research, The University of Western
Australia, city=Perth, postcode=6009, state=Western Australia,
country=Australia

[label4]organization=Department of Electrical, Electronic & Computer
Engineering, School of Engineering, The University of Western Australia,
city=Perth, postcode=6009, state=Western Australia, country=Australia

[label7]organization=Institute of Physics, Faculty of Physics, Astronomy and
Informatics, Nicolaus Copernicus University in Toruń, city=Grudziadzka 5,
postcode=87-100, state=Toruń, country=Poland

[label5]organization=Department of Mechanical Engineering, Massachusetts
Institute of Technology, city=Cambridge, postcode=02139, state=Massachusetts,
country=USA

[label6]organization=Alfred E. Mann Department of Biomedical Engineering,
University of Southern California, city=Los Angeles, postcode=90089,
state=California, country=USA

##  1 Introduction

Generative artificial intelligence (AI) is at the forefront of the latest wave
of advances in deep learning technology. Generative AI is the backbone for
much of the technology that has pervaded our society, such as large language
models and AI-driven image and video generators. Generative models seek to
understand how data is generated and replicate the data generation process to
synthesize new data as desired. The essence of generative models is as
follows: given data sampled from an underlying probability distribution,
inaccessible at the outset, generative models learn to sample new realizations
from the same probability distribution [1]. Mathematically, a generative model
is a means to approximate the underlying probability distribution, say
$\mathrm{p}_{{\mathbf{X}}}$, using realizations of ${\mathbf{X}}$ that are
available at hand. Several types of generative models have been developed,
including generative adversarial networks (GANs), variational autoencoders
(VAEs), autoregressive models (ARMs), normalizing flows (NFs), energy-based
models (EBMs), and diffusion models. These models vary in the way they
approximate the underlying probability distribution and whether they construct
an explicit or implicit approximation [1]. Among them, diffusion models are
emerging as the leading choice and have achieved state-of-the-art performance
in image synthesis [2]. Indeed, diffusion models are used by the latest
generation of AI-based text, image, and video generators [3, 4, 5].

Sometimes, the generation process must be conditioned on inputs or covariates,
say ${\mathbf{Y}}$. For instance, when generating images based on text prompts
[4]. Mathematically framed, a conditional generative model must learn the
conditional distribution $\mathrm{p}_{{\mathbf{X}}|{\mathbf{Y}}}$ using
pairwise realizations from the underlying joint distribution of ${\mathbf{X}}$
and ${\mathbf{Y}}$. In particular, conditional diffusion models are
increasingly being used across several fields in science and engineering; see
recent reviews [6, 7, 8, 9, 10]. Some notable applications in computational
sciences and engineering include material microstructure generation based on
various microstructural descriptors [11] and desired nonlinear behavior [12],
metamaterial synthesis to obtain specified stress-strain response [13],
topology optimization under constraints [14], reconstruction of turbulent flow
fields [15], and synthesizing high-fidelity simulations of turbulent flow from
low-fidelity data [16].

Clearly, the conditional generative modeling resembles the process of solving
an inverse problem within the Bayesian paradigm where ${\mathbf{X}}$ and
${\mathbf{Y}}$ now represent the quantity that must be inferred and the
measurements available, respectively [17]. Thus, motivated by the recent wide-
scale adoption of diffusion models, we proposed to use conditional diffusion
models to solve a class of inverse problems in mechanics that involve
inferring the spatially varying material constitutive parameters from full-
field measurements of the mechanical response, such as displacements or
strains. Inverse problems of this type have many applications in material
characterization [18], structural health monitoring and damage detection [19,
20], and medical diagnosis [21, 22]. We are particularly motivated by the
applications of this type of inverse problem in medical diagnosis, which
arises when using elastography-based imaging techniques that assess the state
of tissue using its mechanical properties, mainly to differentiate between
healthy and diseased tissue [21, 23].

Solving inverse problems of the aforementioned type can be challenging. First,
the inferred quantity is a field that must be discretized to make the inverse
problem tractable, leading to a high-dimensional inverse problem. Second, the
forward model is complicated. Mechanistic models that simulate the
deformations of the specimen in consideration may be a high-fidelity black-box
model. Moreover, the measurement noise in these applications may be complex—
non-Gaussian and non-additive. These factors combine to make Bayesian
inference difficult since standard inference tools, such as Markov chain Monte
Carlo (MCMC) methods and variational inference, may not perform well in the
face of these challenges. However, the recent success of conditional diffusion
models in solving high-dimensional inverse problems offers promise [24, 25,
26, 27, 28].

We propose to use conditional score-based diffusion models — a particular type
of diffusion model that leverages deep neural networks to approximate the
_score function_ of the target posterior distribution [29, 30, 31]. In the
context of inverse problems, the score function is the logarithm of the
posterior density. The neural network approximating the score function is
called the score network [31]. Approximating the score function allows greater
flexibility in choosing the architecture for the score network in the absence
of additional constraints necessary to model a valid probability density. For
instance, the score network need not be invertible [31], unlike normalizing
flows [32]. Moreover, a single neural network can approximate the posterior
distribution for different realizations of the covariate ${\mathbf{Y}}$, which
only requires adding a channel to the input of the neural network [33, 34,
35]. This also allows for re-using the trained conditional score-based
diffusion for different realizations of ${\mathbf{Y}}$ during testing.
Significantly, we show that the score network can be trained using a
regression-type objective, see Eq. 15, which is much more stable than
adversarial losses used to train GANs. Moreover, training only requires
samples from the joint distribution between ${\mathbf{X}}$ and ${\mathbf{Y}}$.
The training dataset can be easily created by simulating the forward model for
different realizations of ${\mathbf{X}}$ sampled from a suitable prior or
collected a priori through experiments. Thus, the forward model can be
arbitrarily complex, black-box, and incompatible with automatic
differentiation. This lends great utility to the proposed approach since most
forward models for the inverse problems of practical interest are likely to
have one or more of the aforementioned properties. Moreover, the training
dataset is entirely synthetic and does not require access to experimental
data, which can be scarce in some applications. After training the score
network, we can sample the posterior distribution for any realization of
${\mathbf{Y}}$, perhaps the result of an experiment, using Langevin dynamics
MCMC [31, 36].

We note that recent work by Jacobsen et al. [37] uses diffusion models for
full-field inversion from sparse measurements. Jacobsen et al. [37] primarily
focuses on enforcing physical consistency to the realizations generated using
the conditional diffusion model, which is not the focus of the current work.
The realizations obtained using conditional score-based diffusion models can
be enhanced by enforcing physical consistency using ideas developed in [37].
Similarly, Bastek et al. [38] propose to enforce physical consistency during
training of the score network, which can also be used to augment the proposed
approach. However, enforcing physical consistency by discretizing the
underlying forward model may not be straightforward in many applications
involving nonlinear phenomena and complex geometries. Therefore, we do not
pursue informing the score network or the sampling process with physics
knowledge in this work. Indeed, the proposed approach is agnostic to the
physics of the forward process, i.e., we assume that the forward model is
black-box.

We also note that deep generative models have been used to solve inverse
problems [17, 39]. One line of work uses deep generative models as data
informative priors for solving inverse problems [40, 41, 42]. Similar to this
work, conditional generative models constructed using GANs [43, 44, 45] and
NFs [46] have been used to sample the posterior in inverse problems. NFs have
also been used for variational Bayesian inference [47, 48]. More recently,
modular Bayesian frameworks have also been developed that use coupled
generative models to perform Bayesian inference [49, 50]. Compared to
diffusion models, other generative models suffer from various limitations,
such as mode collapse in GANs, poor synthesis quality in VAEs, and the need
for specialized architecture with large memory footprints in NFs [33]. Thus,
there is now a growing body of work that utilizes diffusion models to solve
inverse problems [25, 51, 24, 26, 27, 52]. These works use the gradient of the
likelihood function to influence the sample generation process. However,
computing the gradient of the likelihood function usually requires taking
gradients of the black-box nonlinear forward model, which can be
computationally expensive. Thus, efforts have been made to approximate the
gradient of the likelihood function. In contrast, we use the score network to
directly approximate the score function, similar to [33]. Whereas Saharia et
al. [33] considers the problem of image super-resolution, we are motivated by
inverse problems arising in solid mechanics.

The remainder of this paper is organized as follows. In Section 2, we set up
the inverse problem we intend to solve and provide the necessary background on
score-based diffusion models. We introduce our proposed approach to solve
inverse problems using conditional score-based diffusion models in Section 3
where Section 3.1 develops the training objective, Section 3.2 outlines the
sampling scheme, Section 3.3 provides guidelines to choose various hyper-
parameters associated with the model and sampling, Section 3.4 provides design
details pertaining to the neural network used to approximate the score
function, and Section 3.5 provides an overview of the various steps in the
proposed approach. In Section 4, we use the proposed approach to solve various
inverse problems that concern the inference of a spatially varying
constitutive parameter from noisy measurements of sample deformation. The
applications include high-dimensional inverse problems with synthetic (Section
4.1) and experimental data (Section 4.2). Section 5 ends with conclusions and
future directions.

##  2 Background

###  2.1 Problem setup

A typical inverse problem in mechanics involves inferring the spatially
varying material constitutive parameter(s) of a specimen from measurements,
likely noisy and possibly sparse, of its mechanical response under some
controlled loading. The material constitutive parameter(s) that must be
inferred is a continuous field quantity and discretizing it helps make the
inverse problem tractable. We denote using ${\mathbf{X}}$ the
${n_{\mathcal{X}}}$-dimensional vector of material constitutive parameter(s)
that has to be inferred. The ${i}^{\text{th}}$ component of ${\mathbf{X}}$
denotes the material constitutive parameter at some unique physical coordinate
on the specimen. Similarly, we denote using ${\mathbf{Y}}$ the
${n_{\mathcal{Y}}}$-dimensional prediction vector of the specimen’s response.
Associated with the inverse problem is a mechanistic model $\mathit{F}$, which
we can evaluate to obtain the prediction for a given input ${\mathbf{X}}$.
$\mathit{F}$ is likely to be a computational model which we can query to
obtain the response of a specimen for a given input spatial distribution of
the material constitutive parameter ${\mathbf{x}}$, and initial and boundary
conditions. This computational model can be a high-fidelity black-box model,
making it incompatible with automatic differentiation and, therefore,
impossible to interface with deep learning libraries. It can also be
probabilistic, thereby accounting for the effect of the stochastic measurement
moise and model error. Finally, the goal of the inverse problem is to infer
${\mathbf{X}}$ from a measurement $\hat{{\mathbf{y}}}$, a realization of
${\mathbf{Y}}$, corrupted with noise.

###  2.2 Bayesian inference

Bayesian inference is a popular framework for solving inverse problems in a
probabilistic setting. The Bayesian treatment of an inverse problem begins
with the designation of ${\mathbf{X}}$ and ${\mathbf{Y}}$ as random vectors.
Herein, $\mathrm{p}_{{\mathbf{X}}}$ denotes the density of the probability
distribution of ${\mathbf{X}}$; the latter is popularly know as the _prior
distribution_. Now, let $\mathrm{p}_{{\mathbf{Y}}|{\mathbf{X}}}$ denote the
density corresponding to the probability distribution of ${\mathbf{Y}}$
conditioned on ${\mathbf{X}}$, also known as the _likelihood_ distribution,
induced by the forward model

${\mathbf{y}}=\mathcal{F}({\mathbf{x}};\bm{\eta}),$ (1)

where ${\mathbf{x}}$ and ${\mathbf{y}}$ denote realizations of ${\mathbf{X}}$
and ${\mathbf{Y}}$, respectively, and $\bm{\eta}$ is a random vector
representing measurement noise and modeling errors. Drawing a realization of
${\mathbf{Y}}$ given a realization of ${\mathbf{X}}$ using Eq. 1 involves
evaluating the computational model $\mathit{F}$ and sampling $\bm{\eta}$ from
an appropriate distribution. Herein, we assume that we have access to the
likelihood, which implies that we can sample $\bm{\eta}$. Now, given the
measurement $\hat{{\mathbf{y}}}$, which is a realization of the random vector
${\mathbf{Y}}$, the goal of Bayesian inference is to obtain the conditional
probability distribution, also known as the _posterior_ distribution, with
density
$\mathrm{p}_{{\mathbf{X}}|{\mathbf{Y}}}\\!\left(\cdot|{\mathbf{Y}}=\hat{{\mathbf{y}}}\right)$
as follows:

$\mathrm{p}_{{\mathbf{X}}|{\mathbf{Y}}}\\!\left({\mathbf{X}}={\mathbf{x}}|{\mathbf{Y}}=\hat{{\mathbf{y}}}\right)\propto\mathrm{p}_{{\mathbf{X}}|{\mathbf{Y}}}\\!\left({\mathbf{Y}}=\hat{{\mathbf{y}}}|{\mathbf{X}}={\mathbf{x}}\right)\mathrm{p}_{{\mathbf{X}}}\\!\left({{\mathbf{X}}={\mathbf{x}}}\right).$
(2)

In the remainder of this manuscript, we drop the random variables in the
argument of the respective densities to simplify notation. For instance,
$\mathrm{p}_{{\mathbf{X}}}\\!\left({{\mathbf{x}}}\right)$ will mean
$\mathrm{p}_{{\mathbf{X}}}\\!\left({{\mathbf{X}}={\mathbf{x}}}\right)$.

As simple as the Bayes’ rule Eq. 2 is, its application can be very
challenging. We cannot evaluate the posterior analytically unless the prior
and the likelihood distributions form a conjugate pair. A conjugate prior and
likelihood pair is unlikely to arise in applications relevant to the current
work since the noise may be non-additive and the computational model
$\mathit{F}$ may be non-linear. Hence, we must approximate the posterior
distribution. One way is to draw samples from the posterior distribution,
which we can use to compute posterior statistics and posterior predictive
quantities. However, sampling high-dimensional posteriors is far from
straightforward.

###  2.3 Langevin dynamics and the score function

Markov chain Monte Carlo (MCMC) algorithms have been the staple of Bayesian
machinery. MCMC methods are useful for sampling unnormalized probability
models such as the posterior distribution in inverse problems. Early MCMC
methods based on the Metropolis-Hastings (MH) algorithm [53, 54] use a
proposal density to advance a Markov chain with an invariant distribution
similar to the target posterior. However, in high-dimensional settings, simple
MH-type MCMC methods are known to be inefficient as randomly chosen new states
fail to explore high-likelihood regions quickly.

Langevin dynamics is an improvement over MH-MCMC, and uses a proposal density
obtained from discretizing Langevin diffusion equation with a drift term
informed by the gradient of the target posterior distribution. The resulting
stochastic differential equation is defined as:

$\mathrm{d}{\mathbf{X}}_{t}=\frac{\nabla_{{\mathbf{x}}}\log\mathrm{p}_{{\mathbf{X}}|{\mathbf{Y}}}\\!\left({{\mathbf{x}}|\hat{{\mathbf{y}}}}\right)}{2}\mathrm{d}t+\mathrm{d}{\mathbf{B}}_{t},$
(3)

where ${\mathbf{B}}$ represent as ${n_{\mathcal{X}}}$-dimensional Brownian
motion. The first term on the right hand side of Eq. 3 is popularly known as
the _score_ function and it has two noteworthy properties. First, the score
function is vector valued, i.e.,
$\nabla_{{\mathbf{x}}}\log\mathrm{p}_{{\mathbf{X}}|{\mathbf{Y}}}\\!\left({{\mathbf{x}}|\hat{{\mathbf{y}}}}\right):\mathbb{R}^{{n_{\mathcal{X}}}}\to\mathbb{R}^{{n_{\mathcal{X}}}}$.
Second, the score function does not depend on the normalization constant of
the posterior distribution, i.e.,

$\mathrm{p}_{{\mathbf{X}}|{\mathbf{Y}}}\\!\left({{\mathbf{x}}|\hat{{\mathbf{y}}}}\right)\propto\mathrm{q}_{{\mathbf{X}}|{\mathbf{Y}}}\\!\left({\mathbf{x}}|\hat{{\mathbf{y}}}\right)\implies\nabla_{\tilde{{\mathbf{x}}}}\log\mathrm{p}_{{\mathbf{X}}|{\mathbf{Y}}}\\!\left({{\mathbf{x}}|\hat{{\mathbf{y}}}}\right)=\nabla_{{\mathbf{x}}}\log\mathrm{q}_{{\mathbf{X}}|{\mathbf{Y}}}\\!\left({\mathbf{x}}|\hat{{\mathbf{y}}}\right),$
(4)

where the density
$\mathrm{q}_{{\mathbf{X}}|{\mathbf{Y}}}\\!\left({\mathbf{x}}|\hat{{\mathbf{y}}}\right)$
is unnormalized, i.e.,
$\int\mathrm{q}_{{\mathbf{X}}|{\mathbf{Y}}}\\!\left({\mathbf{x}}|\hat{{\mathbf{y}}}\right)\;\mathrm{d}{\mathbf{x}}\neq
1$. Further, a first order Euler-Maruyama discretization of Eq. 3 yields the
following iterative proposal mechanism:

${\mathbf{x}}_{t+1}={\mathbf{x}}_{t}+\frac{\epsilon}{2}\nabla_{{\mathbf{x}}}\log\mathrm{p}_{{\mathbf{X}}|{\mathbf{Y}}}\\!\left({{\mathbf{x}}_{t}|\hat{{\mathbf{y}}}}\right)+\sqrt{\epsilon}{\mathbf{z}}_{t},$
(5)

where $\epsilon$ is the integration step size and ${\mathbf{z}}_{t}$ are
independent and identical realizations of standard normal variable
${\mathbf{Z}}$. Eq. 5 converges to the invariant distribution
$\mathrm{p}_{{\mathbf{X}}|{\mathbf{Y}}}$ as time $t\to\infty$. As evident from
Eqs. 3 and 5, the score function
$\nabla_{{\mathbf{x}}}\log\mathrm{p}_{{\mathbf{X}}|{\mathbf{Y}}}$ is crucial
for sampling the posterior distribution using Langevin dynamics.

##  3 Bayesian inference using conditional score-based diffusion models

Bayesian inference methods try to approximate the posterior distribution of
${\mathbf{X}}$ for a given measurement $\hat{{\mathbf{y}}}$ using samples
drawn from it. Eqs. 3 and 5 are useful for devising an MCMC method for
sampling the posterior distribution. However, the ‘ _true_ ’ score function
remains inaccessible in practical applications. Using conditional score-based
diffusion generative models, we approximate the score function using a neural
network $\mathbf{s}({\mathbf{x}},{\mathbf{y}};\bm{\theta})$, popularly known
as the _score network_ , where $\bm{\theta}$ represents the parameters of the
neural network. In this section, we develop the training objective for the
score network, introduce the MCMC-based sampler called Annealed Langevin
dynamics [31] that is useful for sampling target posteriors using a trained
score network, provide guidelines regarding hyper-parameter selection and
design choices for the score network, and finally provide an overview of the
entire workflow.

###  3.1 Training the score network using denoising score matching

Consider a realization ${\mathbf{y}}$ of ${\mathbf{Y}}$ associated with which
is the target posterior distribution
$\mathrm{p}_{{\mathbf{X}}|{\mathbf{Y}}}\\!\left({{\mathbf{x}}|{\mathbf{y}}}\right)$
with the score function
$\nabla_{{\mathbf{x}}}\log\mathrm{p}_{{\mathbf{X}}|{\mathbf{Y}}}\\!\left({{\mathbf{x}}|{\mathbf{y}}}\right)$.
We want the score network $\mathbf{s}({\mathbf{x}},{\mathbf{y}};\bm{\theta})$
to be a good approximation to the score function over the entire support of
${\mathbf{X}}$. So, a potential candidate objective function to train the
score network is the score matching objective,

$\ell_{\mathrm{SM}}({\mathbf{y}};\bm{\theta})=\frac{1}{2}\mathbb{E}_{{\mathbf{X}}\sim\mathrm{p}_{{\mathbf{X}}|{\mathbf{Y}}}}\Big{[}\lVert\mathbf{s}({\mathbf{x}},{\mathbf{y}};\bm{\theta})-\nabla_{{\mathbf{x}}}\log\mathrm{p}_{{\mathbf{X}}|{\mathbf{Y}}}\\!\left({{\mathbf{x}}|{\mathbf{y}}}\right)\rVert_{2}^{2}\Big{]}.$
(6)

One can show that, under some mild regularity conditions, the score matching
objective in Eq. 6 is equivalent, up to a constant, to the following

$\ell_{\mathrm{SM}}({\mathbf{y}};\bm{\theta})=\mathbb{E}_{{\mathbf{X}}\sim\mathrm{p}_{{\mathbf{X}}|{\mathbf{Y}}}}\Big{[}\frac{1}{2}\lVert\mathbf{s}({\mathbf{x}},{\mathbf{y}};\bm{\theta})\rVert_{2}^{2}+\mathrm{Tr}(\nabla_{{\mathbf{x}}}\mathbf{s}({\mathbf{x}},{\mathbf{y}};\bm{\theta}))\Big{]},$
(7)

where $\mathrm{Tr}(\cdot)$ is the trace operator [29, 31]. However, Eq. 7 is
not scalable to high-dimensional problems since it involves computing the
gradient of the output of the score network with respect to its input, the
computational cost of which will scale linearly with the dimension
${n_{\mathcal{X}}}$ of the inverse problem. To circumvent this issue,
Hyvärinen and Dayan [29] proposed a _denoising score matching_ objective that
altogether bypasses computing
$\nabla_{{\mathbf{x}}}\mathbf{s}({\mathbf{x}},{\mathbf{y}};\bm{\theta})$.

Denoising score matching avoids the computation of
$\nabla_{{\mathbf{x}}}\mathbf{s}({\mathbf{x}},{\mathbf{y}};\bm{\theta})$ by
first perturbing the random variable ${\mathbf{X}}$, conditioned on
${\mathbf{Y}}$, to obtain $\tilde{{\mathbf{X}}}$ using a perturbation kernel
$\mathrm{p}_{\tilde{{\mathbf{X}}}|{\mathbf{X}}}$. Note that this makes
$\tilde{{\mathbf{X}}}$ and ${\mathbf{Y}}$ conditionally independent given
${\mathbf{X}}$. Now, the denoising score matching objective, for a given
realization ${\mathbf{y}}$ of ${\mathbf{Y}}$, is defined as

$\ell_{\mathrm{DSM}}({\mathbf{y}};\bm{\theta})=\frac{1}{2}\mathbb{E}_{\tilde{{\mathbf{X}}}\sim\mathrm{p}_{\tilde{{\mathbf{X}}}|{\mathbf{Y}}}}\Big{[}\lVert\mathbf{s}(\tilde{{\mathbf{x}}},{\mathbf{y}};\bm{\theta})-\nabla_{\tilde{{\mathbf{x}}}}\log\mathrm{p}_{\tilde{{\mathbf{X}}}|{\mathbf{Y}}}\\!(\tilde{{\mathbf{x}}}|{\mathbf{y}})\rVert_{2}^{2}\Big{]}.$
(8)

Using the fact that

$\mathrm{p}_{\tilde{{\mathbf{X}}}|{\mathbf{Y}}}(\tilde{{\mathbf{x}}}|{\mathbf{y}})=\int\mathrm{p}_{\tilde{{\mathbf{X}}}|{\mathbf{X}}}\\!\left({\tilde{{\mathbf{x}}}|{\mathbf{x}}}\right)\;\mathrm{p}_{{\mathbf{X}}|{\mathbf{Y}}}({\mathbf{x}}|{\mathbf{y}})\;\mathrm{d}{\mathbf{x}},$
(9)

one can show that the denoising score matching objective $\ell_{\mathrm{DSM}}$
is equivalent to [30]:

$\ell_{\mathrm{DSM}}({\mathbf{y}};\bm{\theta})\equiv\frac{1}{2}\mathbb{E}_{{\mathbf{X}}\sim\mathrm{p}_{{\mathbf{X}}|{\mathbf{Y}}}}\Bigg{[}\mathbb{E}_{\tilde{{\mathbf{X}}}\sim\mathrm{p}_{\tilde{{\mathbf{X}}}|{\mathbf{X}}}}\Big{[}\lVert\mathbf{s}(\tilde{{\mathbf{x}}},{\mathbf{y}};\bm{\theta})-\nabla_{\tilde{{\mathbf{x}}}}\log\mathrm{p}_{\tilde{{\mathbf{X}}}|{\mathbf{X}}}\\!\left({\tilde{{\mathbf{x}}}|{\mathbf{x}}}\right)\rVert_{2}^{2}\Big{]}\Bigg{]},$
(10)

up to an additive constant that does not depend on $\bm{\theta}$; see Appendix
A for the derivation. The optimally trained score network with parameters
$\bm{\theta}^{\ast}$, trained using Eq. 10, will satisfy
$\mathbf{s}(\tilde{{\mathbf{x}}},{\mathbf{y}};\bm{\theta}^{\ast})\approx\nabla_{\tilde{{\mathbf{x}}}}\log\mathrm{p}_{\tilde{{\mathbf{X}}}|{\mathbf{Y}}}(\tilde{{\mathbf{x}}}|{\mathbf{y}})$
and
$\mathrm{p}_{\tilde{{\mathbf{X}}}|{\mathbf{Y}}}(\tilde{{\mathbf{x}}}|{\mathbf{y}})\approx\mathrm{p}_{{\mathbf{X}}|{\mathbf{Y}}}({\mathbf{x}}|{\mathbf{y}})$,
when the perturbation kernel bandwidth goes to zero. Thus, the training
objective Eq. 10 offers a path to approximating the score function of the
target posterior distribution for a specified realization of ${\mathbf{y}}$
using a score network, provided we choose a very narrow perturbation kernel.

We note that we wish to use the score network to approximate the score of the
posterior distribution over the entire support of ${\mathbf{Y}}$. Hence, we
marginalize ${\mathbf{Y}}$ out of
$\ell_{\mathrm{DSM}}({\mathbf{y}};\bm{\theta})$ to obtain the training
objective

$\begin{split}\mathcal{L}_{\mathrm{DSM}}(\bm{\theta})&=\mathbb{E}_{{\mathbf{Y}}\sim\mathrm{p}_{{\mathbf{Y}}}}\big{[}\ell_{\mathrm{DSM}}({\mathbf{y}};\bm{\theta})\big{]}\\\
&\equiv\frac{1}{2}\mathbb{E}_{{\mathbf{Y}}\sim\mathrm{p}_{{\mathbf{Y}}}}\Bigg{[}\mathbb{E}_{{\mathbf{X}}\sim\mathrm{p}_{{\mathbf{X}}|{\mathbf{Y}}}}\bigg{[}\mathbb{E}_{\tilde{{\mathbf{X}}}\sim\mathrm{p}_{\tilde{{\mathbf{X}}}|{\mathbf{X}}}}\Big{[}\lVert\mathbf{s}(\tilde{{\mathbf{x}}},{\mathbf{y}};\bm{\theta})-\nabla_{\tilde{{\mathbf{x}}}}\log\mathrm{p}_{\tilde{{\mathbf{X}}}|{\mathbf{X}}}\\!\left({\tilde{{\mathbf{x}}}|{\mathbf{x}}}\right)\rVert_{2}^{2}\Big{]}\bigg{]}\Bigg{]}\\\
&=\frac{1}{2}\mathbb{E}_{{\mathbf{X}},{\mathbf{Y}}\sim\mathrm{p}_{{\mathbf{X}}{\mathbf{Y}}}}\bigg{[}\mathbb{E}_{\tilde{{\mathbf{X}}}\sim\mathrm{p}_{\tilde{{\mathbf{X}}}|{\mathbf{X}}}}\Big{[}\lVert\mathbf{s}(\tilde{{\mathbf{x}}},{\mathbf{y}};\bm{\theta})-\nabla_{\tilde{{\mathbf{x}}}}\log\mathrm{p}_{\tilde{{\mathbf{X}}}|{\mathbf{X}}}\\!\left({\tilde{{\mathbf{x}}}|{\mathbf{x}}}\right)\rVert_{2}^{2}\Big{]}\bigg{]},\end{split}$
(11)

which is remarkable for one very particular reason. Note, estimating the
objective $\mathcal{L}_{\mathrm{DSM}}(\bm{\theta})$ requires computing two
expectations. The outer expectation is with respect to the joint distribution
of ${\mathbf{X}}$ and ${\mathbf{Y}}$. We can easily sample the joint by first
sampling ${\mathbf{X}}$ from the prior $\mathrm{p}_{{\mathbf{X}}}$
subsequently ${\mathbf{Y}}$ from the likelihood
$\mathrm{p}_{{\mathbf{Y}}|{\mathbf{X}}}$, which simply requires the evaluation
of the forward model Eq. 1. The inner expectation is with respect to the
conditional distribution $\mathrm{p}_{\tilde{{\mathbf{X}}}|{\mathbf{X}}}$,
which we can choose to be a tractable distribution that can be easily sampled.
Thus, the Monte Carlo approximation of
$\mathcal{L}_{\mathrm{DSM}}(\bm{\theta})$, which we denote as
$\mathcal{L}_{\mathrm{eDSM}}(\bm{\theta})$, is given by

$\mathcal{L}_{\mathrm{eDSM}}(\bm{\theta})=\frac{1}{2N_{\mathrm{train}}}\sum_{i=1}^{N_{\mathrm{train}}}\mathbb{E}_{\tilde{{\mathbf{X}}}\sim\mathrm{p}_{\tilde{{\mathbf{X}}}|{\mathbf{X}}}}\Big{[}\lVert\mathbf{s}(\tilde{{\mathbf{x}}},{\mathbf{y}}^{(i)};\bm{\theta})-\nabla_{\tilde{{\mathbf{x}}}}\log\mathrm{p}_{\tilde{{\mathbf{X}}}|{\mathbf{X}}}\\!\left({\tilde{{\mathbf{x}}}|{\mathbf{x}}^{(i)}}\right)\rVert_{2}^{2}\Big{]},$
(12)

where ${\mathbf{x}}^{(i)}$ and ${\mathbf{y}}^{(i)}$ are the ${i}^{\text{th}}$
realizations of ${\mathbf{x}}$ and ${\mathbf{y}}$ sampled from the joint
$\mathrm{p}_{{\mathbf{X}}{\mathbf{Y}}}$. Moreover,
${\mathbf{y}}^{(i)}=\mathcal{F}({\mathbf{x}}^{(i)};\bm{\eta}^{(i)})$, where
$\bm{\eta}^{(i)}$ is the ${i}^{\text{th}}$ realization of the noise vector.

There are three practical choices made to aid in the approximation of the
score function by the score network. First, the perturbation kernel is assumed
to be Gaussian with covariance equal to
$\sigma^{2}\mathbb{I}_{{n_{\mathcal{X}}}}$, i.e.,

$\mathrm{p}_{\tilde{{\mathbf{X}}}|{\mathbf{X}}}\\!\left({\tilde{{\mathbf{x}}}|{\mathbf{x}}}\right)=\mathcal{N}(\tilde{{\mathbf{x}}};{\mathbf{x}},\sigma^{2}\mathbb{I}_{{n_{\mathcal{X}}}})\propto\exp\left(-\frac{\lVert\tilde{{\mathbf{x}}}-{\mathbf{x}}\rVert_{2}^{2}}{2\sigma^{2}}\right)\implies\nabla_{{\mathbf{x}}}\log\mathrm{p}_{\tilde{{\mathbf{X}}}|{\mathbf{X}}}\\!\left({\tilde{{\mathbf{x}}}|{\mathbf{x}}}\right)=\frac{{\mathbf{x}}-\tilde{{\mathbf{x}}}}{\sigma^{2}},$
(13)

where $\sigma$ is referred to as noise level (different from measurement
noise) in the literature [31, 36] because Eq. 13 represents adding independent
and identically distributed Gaussian noise with variance $\sigma^{2}$ to a
realization ${\mathbf{x}}$. Thus, the bandwidth of the perturbation kernel
$\mathrm{p}_{\tilde{{\mathbf{X}}}|{\mathbf{X}}}$ scales linearly with
$\sigma$. We make the dependence on $\sigma$ explicit by introducing the
subscript $\sigma$ to $\tilde{{\mathbf{X}}}$. For instance, herein, we rewrite
$\mathrm{p}_{\tilde{{\mathbf{X}}}|{\mathbf{X}}}$ as
$\mathrm{p}_{\tilde{{\mathbf{X}}}_{\sigma}|{\mathbf{X}}}$. Second, Song and
Ermon [31] suggest that the score function be simultaneously estimated across
different noise levels $\\{\sigma_{i}\\}_{i=1}^{L}$, where
$\sigma_{1}>\sigma_{2}>\ldots>\sigma_{L}>0$. This is akin to ‘annealing’, and,
in this case, helps sampling from
$\mathrm{p}_{\tilde{{\mathbf{X}}}|{\mathbf{Y}}}(\tilde{{\mathbf{x}}}|{\mathbf{y}})$
by ensuring that less significant modes are smoothed out initially. In fact,
sampling will commence from a near uni-modal multivariate Gaussian
distribution for a very large value of $\sigma_{1}$. Song and Ermon [31] also
demonstrate that annealing helps in the exploration of low density regions of
the joint distribution $\mathrm{p}_{{\mathbf{X}}{\mathbf{Y}}}$ when sampling.
However, working with $L$ noise levels instead of one also entails a re-
parameterization of the score network as
$\mathbf{s}(\tilde{{\mathbf{x}}},{\mathbf{y}},\sigma;\bm{\theta})$ as $\sigma$
becomes an additional input to it. However, the input $\sigma$ is used
differently from the other inputs; we discuss how the score network uses
$\sigma$ in Section 3.4. Third, the inner expectation in Eq. 12, with respect
to
$\mathbb{E}_{\tilde{{\mathbf{X}}}\sim\mathrm{p}_{\tilde{{\mathbf{X}}}|{\mathbf{X}}}}$,
is approximated using a single realization randomly sampled on the fly during
training.

On assimilating all these practices into Eq. 12, the new and final training
objective becomes

$\mathcal{L}_{\mathrm{eDSM}}(\bm{\theta})=\frac{1}{2LN_{\mathrm{train}}}\sum_{j=1}^{L}\sum_{i=1}^{N_{\mathrm{train}}}\lambda(\sigma_{j})\Bigg{\lVert}\mathbf{s}(\tilde{{\mathbf{x}}}_{\sigma_{j}}^{(i)},{\mathbf{y}}^{(i)},\sigma_{j};\bm{\theta})+\frac{\tilde{{\mathbf{x}}}_{\sigma_{j}}^{(i)}-{\mathbf{x}}^{(i)}}{\sigma_{j}^{2}}\Bigg{\rVert}_{2}^{2},$
(14)

where $\tilde{{\mathbf{x}}}_{\sigma_{j}}^{(i)}$ is the ${i}^{\text{th}}$
realization of $\tilde{{\mathbf{X}}}_{\sigma_{j}}$ sampled from
$\mathrm{p}_{\tilde{{\mathbf{X}}}_{\sigma_{j}}|{\mathbf{X}}}$ and
$\lambda(\sigma_{j})$ is a weighting factor chosen equal to $\sigma^{2}_{j}$
[31, 36], which in turn makes

$\mathcal{L}_{\mathrm{eDSM}}(\bm{\theta})=\frac{1}{2LN_{\mathrm{train}}}\sum_{j=1}^{L}\sum_{i=1}^{N_{\mathrm{train}}}\Bigg{\lVert}\sigma_{j}\mathbf{s}(\tilde{{\mathbf{x}}}_{\sigma_{j}}^{(i)},{\mathbf{y}}^{(i)},\sigma_{j};\bm{\theta})+\frac{\tilde{{\mathbf{x}}}_{\sigma_{j}}^{(i)}-{\mathbf{x}}^{(i)}}{\sigma_{j}}\Bigg{\rVert}_{2}^{2}.$
(15)

A score network of sufficient capacity trained using
$\mathcal{L}_{\mathrm{eDSM}}(\bm{\theta})$, as defined in Eq. 15, will
approximate the score function across all noise levels. Let
$\bm{\theta}^{\ast}$ denote the optimal parameters of the score network such
that

$\bm{\theta}^{\ast}=\arg\min_{\bm{\theta}}\mathcal{L}_{\mathrm{eDSM}}(\bm{\theta}).$
(16)

Using the trained score network
$\mathbf{s}(\cdot,\cdot,\cdot;\bm{\theta}^{\ast})$, Annealed Langevin dynamics
[31] can be used to sample from the target posterior distribution.

###  3.2 Sampling using Annealed Langevin Dynamics

From Eq. 9, note, $\mathrm{p}_{\tilde{{\mathbf{X}}}|{\mathbf{Y}}}$ is the
outcome of a convolution between the target posterior
$\mathrm{p}_{{\mathbf{X}}|{\mathbf{Y}}}$ and the zero-mean Gaussian kernel
with covariance matrix $\sigma^{2}\mathbb{I}_{{n_{\mathcal{X}}}}$. This is
akin to a ’ _smoothing_ ’ operation of the target posterior. In this work, we
use Annealed Langevin Dynamics, proposed by Song and Ermon [31], to sample
from the target posterior distributions using the learned score network
starting from the smoothest target and gradually reducing the amount of
smoothing; see Algorithm 1.

Input: Trained score network
$\mathbf{s}({\mathbf{x}},{\mathbf{y}},\sigma;\bm{\theta})$, measurement
$\hat{{\mathbf{y}}}$, sampling steps $T$, step size $\epsilon^{2}$, and
perturbation kernel standard deviations
$\left\\{\sigma_{i}\right\\}_{i=1}^{L}$

1 Initialize ${\mathbf{x}}_{0}$ such that its ${i}^{\text{th}}$ component
$\\{{\mathbf{x}}_{0}\\}_{i}\sim\mathcal{U}(0,1)$

2 for $i=1$ to $L$ do

3 $\alpha_{i}=\epsilon\sigma_{i}^{2}/\sigma^{2}_{L}$

4 for $j=1$ to $T$ do

5 Sample ${\mathbf{z}}\sim\mathcal{N}(\bm{0},\mathbb{I}_{{n_{\mathcal{X}}}})$

6
${\mathbf{x}}_{j}={\mathbf{x}}_{j-1}+\frac{\alpha_{i}}{2}\mathbf{s}({\mathbf{x}}_{j-1},\hat{{\mathbf{y}}},\sigma_{i};\bm{\theta})+\sqrt{\alpha_{i}}{\mathbf{z}}$

7

8 end for

9 Set ${\mathbf{x}}_{0}={\mathbf{x}}_{T}$

10

11 end for

12Denoise
${\mathbf{x}}_{\mathrm{final}}={\mathbf{x}}_{0}+\sigma_{L}^{2}s({\mathbf{x}}_{0},\hat{{\mathbf{y}}},\sigma_{L};\bm{\theta})$

Output: Realization ${\mathbf{x}}_{\mathrm{final}}$

Algorithm 1 Annealed Langevin dynamics [36]

###  3.3 Hyper-parameters

There are a few hyper-parameters associated with the conditional score-based
diffusion models. Some of these hyper-parameters are related to the training
of the score network while others are associated with the sampling algorithm
used to generate new realizations. Here we discuss how to choose the various
important hyper-parameters.

The hyper-parameters associated with score estimation include the noise levels
$\sigma_{1}$ through $\sigma_{L}$, and the number $L$ of noise levels.
Following Song and Ermon [36], we choose $\sigma_{1}$ to be larger than the
Euclidean distance between any two points of ${\mathbf{x}}$ within the
training data. This helps smooth out the effect of
$\mathrm{p}_{{\mathbf{X}}|{\mathbf{Y}}}$ from
$\mathrm{p}_{\tilde{{\mathbf{X}}}|{\mathbf{Y}}}$ and sampling can commence
from a near unimodal Gaussian distribution with covariance close to
$\sigma_{1}^{2}\mathbb{I}_{{n_{\mathcal{X}}}}$. This helps avoid particles
from getting stuck near modes of the target distribution and the sampling
space can be properly explored with a limited number of particles. We choose
$\sigma_{L}$ to be a very small value — typically $0.001$ or $0.01$ — so that
the approximation
$\mathrm{p}_{\tilde{{\mathbf{X}}}|{\mathbf{Y}}}(\tilde{{\mathbf{x}}}|{\mathbf{y}})\approx\mathrm{p}_{{\mathbf{X}}|{\mathbf{Y}}}({\mathbf{x}}|{\mathbf{y}})$
is better.

Following Song and Ermon [36], we choose $\sigma_{j}$ to be a geometric
progression between $\sigma_{1}$ and $\sigma_{L}$, i.e., we fix
$\gamma=\sigma_{j-1}/\sigma_{j}$. Moreover, we choose the number $L$ of noise
levels such that particles from
$\mathrm{p}_{\tilde{{\mathbf{X}}}_{\sigma_{j-1}}|{\mathbf{Y}}}$ are good
initializers for the Langevin dynamics MCMC when exploring
$\mathrm{p}_{\tilde{{\mathbf{X}}}_{\sigma_{j}}|{\mathbf{Y}}}$, which will be
true when their is significant overlap between the two distribution. Song and
Ermon [36] suggest choosing $\gamma$ such that:

$\Phi(\sqrt{2{n_{\mathcal{X}}}}(\gamma-1)+3\gamma)-\Phi(\sqrt{2{n_{\mathcal{X}}}}(\gamma-1)-3\gamma)\approx
0.5,$ (17)

where $\Phi$ is the cumulative distribution function of a standard normal
variable, to ensure good overlap between successive smoothened versions of the
posterior. In practice, we first choose $L$, compute $\gamma$ and then verify
that Eq. 17 is satisfied.

Apart from the aforementioned hyper-parameters, there are the usual learning
related hyper-parameters related to the training of the score network, which
includes the learning rate, batch size and number of training iterations.
These hyper-parameters are chosen to avoid over-fitting by monitoring the loss
on a test set, which can be carved out from the initial training set if
necessary.

The hyper-parameters for the Annealed Langevin dynamics include the number of
Langevin steps $T$ and the step size $\epsilon$. Increasing $L$ will increase
the total number of Langevin steps, and the total time it takes to generate
realizations from the target posterior distribution. So, we choose $T$ as
large as possible as the compute budget allows. Finally, we choose $\epsilon$
such that particles from
$\mathrm{p}_{\tilde{{\mathbf{X}}}_{\sigma_{j-1}}|{\mathbf{X}}}$ can reach
their steady state distribution
$\mathrm{p}_{\tilde{{\mathbf{X}}}_{\sigma_{j}}|{\mathbf{X}}}$ within $T$
steps. Song and Ermon [36] suggest choosing a value of $\epsilon$ for which
[36]

$\left(1-\frac{\epsilon}{\sigma^{2}_{L}}\right)^{2T}\left(\gamma^{2}-\frac{2\epsilon}{\sigma^{2}_{L}-\sigma^{2}_{L}\left(1-\frac{\epsilon}{\sigma_{L}^{2}}\right)}\right)+\frac{2\epsilon}{\sigma_{L}^{2}-\sigma^{2}_{L}\left(1-\frac{\epsilon}{\sigma_{L}^{2}}\right)}\approx
1.$ (18)

In practice, we perform a grid search over probable values of $\epsilon$ and
choose the value for which the left hand side of Eq. 18 is closest to 1.

###  3.4 Design choices for the score network

The design choices for the score networks we use in this work follow
suggestions made by Song and Ermon [36]. Primarily, the score network treats
$\sigma_{j}$ unlike an usual input such as ${\mathbf{x}}$ or ${\mathbf{y}}$.
Rather, $\sigma_{j}$ is used to scale the outputs from the network i.e.,
$\sigma_{j}$ divides the outcome from the last layer of the score network.
Moreover, we treat ${\mathbf{x}}$ and ${\mathbf{y}}$ as images as they are
discretized forms of continuous field quantities: ${\mathbf{x}}$ is the
spatially varying material constitutive parameter over the specimen
discretized over an appropriate Cartesian grid, and ${\mathbf{y}}$ is the
measured mechanical response of the specimen over the same Cartesian grid.
Therefore, we use convolution neural network type architectures to construct
the score functions with ${\mathbf{x}}$ and ${\mathbf{y}}$ as separate
channels. More specifically, we use RefineNets [55] to construct the score
networks in this work. If there are multiple measurement types available or
the measurement type itself has more than one component, we simply add
additional channels. For instance, we treat vertical and horizontal
displacement measurements as separate channels. Similarly, in the case of a
time-dependent problem, if vertical displacement measurements are made at
different time instants, the instantaneous measurements are concatenated as
additional channels.

###  3.5 Workflow

In this section, we describe the workflow for solving an inverse problem using
conditional score-based diffusion models. The basic workflow can be divided
into three steps; see Fig. 1.

Figure 1: Schematic diagram of the proposed framework for solving inverse
problems in mechanics using conditional diffusion models. Step 1 involves
creating the dataset by sampling realizations of $({\mathbf{X}},{\mathbf{Y}})$
from the joint distribution $\mathrm{p}_{{\mathbf{X}}{\mathbf{Y}}}$, which
requires multiple evaluations of the forward model $\mathcal{F}$. Step 2
involves the training of the score network, which progressively learns to
predict the noise in realizations of $\tilde{{\mathbf{X}}}$ conditioned on
${\mathbf{Y}}$ for all levels of noise $\sigma$. Finally, Annealed Langevin
dynamics is used to sample from the target posterior distribution in Step 3

####  3.5.1 Step 1: Create the training dataset

The first step involves the creation of the training dataset. For this step,
we first develop a parameteric prior model for ${\mathbf{X}}$, resembling the
first step in a more traditional application of Bayesian inference.
Alternatively, multiple instances of ${\mathbf{X}}$ would also suffice. Next,
for each realization ${\mathbf{x}}$ of ${\mathbf{X}}$ we solve the forward
problem Eq. 1 to obtain a corresponding realization of ${\mathbf{y}}$ of
${\mathbf{Y}}$. Note that, to solve the forward model Eq. 1 we require access
to the distribution of the noise $\bm{\eta}$ in the measurement model. Any
other modeling errors must also be accounted for by Eq. 1. Also note that
simulating Eq. 1 will require simulating the underlying physics-based
mechanics solver, which can be black-box and incompatible with auto-
differentiation. We denote using $\mathcal{S}$ the training dataset consisting
$N_{\mathrm{train}}$ realizations of ${\mathbf{x}}$ and ${\mathbf{y}}$, i.e.,
$\mathcal{S}=\left\\{{\mathbf{x}}^{(i)},{\mathbf{y}}^{(i)}\right\\}_{i=1}^{N_{\mathrm{train}}}$.

######  Remark 1.

We assume that we have access to a sufficient number of realizations of
${\mathbf{X}}$ to build the training dataset $\mathcal{S}$. This may be data
that was previously acquired, or the result of sampling from a parametric
prior developed to capture the uncertainty in ${\mathbf{X}}$. We also assume
that we can sample from the distribution for $\bm{\eta}$ to account for the
measurement noise or modeling errors in the prediction ${\mathbf{Y}}$.

####  3.5.2 Step 2: Train the score network using the training dataset

After creating the training dataset $\mathcal{S}$, we train the score network.
First, we determine the following hyper-parameters— $\sigma_{1}$, $\sigma_{L}$
and $L$. Then we train the score network using a suitable optimization
algorithm and the training objective Eq. 15. As is typical in machine
learning, at each iteration we do not estimate
$\mathcal{L}_{\mathrm{eDSM}}(\bm{\theta})$ using the entire training dataset,
rather we use a small batch of training data. Thus,
$\mathcal{L}_{\mathrm{eDSM}}(\bm{\theta})$ is estimated using
$N_{\mathrm{batch}}$ training data points instead of $N_{\mathrm{train}}$ in
Eq. 15, where $N_{\mathrm{batch}}$ is the batch size. Additionally, the
exponential moving averages of the weights are maintained. Let
$\bm{\theta}^{(i)}$ denote the parameters after the ${i}^{\text{th}}$
iteration and $\bm{\theta}^{\prime}$ denote an independent copy of those
parameters. After every iteration, we update $\bm{\theta}^{\prime}$ as:

$\bm{\theta}^{\prime}=m\bm{\theta}^{\prime}+(1-m)\bm{\theta}^{(i)},$ (19)

where $m=0.999$ is a momentum parameter. This helps suppress any sudden
fluctuations in the training loss, as is typically observed when training
score networks [36]. At the end of training, we use the averaged parameters as
the optimal parameters for the score network.

####  3.5.3 Step 3: Generate samples from the posterior using Annealed
Langevin dynamics

In this step, we sample from the target posterior via the Annealed Langevin
dynamics algorithm shown in Algorithm 1. First, we specify the number of
Langevin steps $T$ and step size $\epsilon$, and use the trained score network
$\mathbf{s}(\cdot,\cdot,\cdot;\bm{\theta}^{\prime})$ to generate new
realizations of ${\mathbf{X}}$ for a given measurement $\hat{{\mathbf{y}}}$.
Algorithm 1 is parallelizable and multiple realizations can be generated
together. We use the generated realizations to estimate the posterior mean and
posterior standard deviation of ${\mathbf{X}}$, with the former serving as the
reconstruction and the latter providing a measure of uncertainty.

######  Remark 2.

Using a trained score network, we can sample the posterior distributions for
different realizations of ${\mathbf{Y}}$ without re-training the score
network. This makes the proposed approach particularly suitable for high
throughput applications requiring repeated solving of the same inverse
problem.

##  4 Results

In this section, we use conditional score-based diffusion models to solve
various large-scale inverse problems in mechanics where the goal is to infer
the Young’s modulus (or a related constitutive parameter) field of a material
specimen subject to external loading. We consider the following problems:

  1. 1.

Section 4.1.1 — Quasi-static elastography problem of inferring Young’s modulus
field of specimens containing a stiff inclusion within a soft homogeneous
background media from displacement measurements.

  2. 2.

Section 4.1.2 — Time-harmonic elastography problem of inferring the shear
modulus field of regions surrounding the human optic nerve head from complex
measurements of a wave field as the specimen is subject to harmonic
excitation.

  3. 3.

Section 4.2.1 — Quasi-static elastography problem similar to Section 4.1.1 but
using experimental data from tissue-mimicking hydrogels [56].

  4. 4.

Section 4.2.2 — Time-dependent elastography problem of inferring Young’s
modulus of a dual layer tissue-mimicking phantom excited by acoustic radiation
force using snapshots of the displacement field taken at regular time
intervals from a laboratory experiment [57].

  5. 5.

Section 4.2.3 — Optical coherence elastography of tumor spheroids (clusters of
cells within a hydrogel) wherein the Young’s modulus field of the specimen is
inferred from the phase difference between optical coherence microscopy scans
as the specimen undergoes compression. This example also involves experimental
data and its details can be found in Foo et al. [58].

Table 1: Summary of inverse problems considered in this work | Inverse problem |  Forward physics | Dimensionality |  Training dataset size $N_{\mathrm{train}}$
---|---|---|---|---
| ${n_{\mathcal{X}}}$ | ${n_{\mathcal{Y}}}$
Synthetic data |  Circular inclusion (Section 4.1.1) |  Linear, quasi-static |  56$\times$56 |  56$\times$56 | 10,000
Optic nerve head (Section 4.1.2) |  Linear, time-harmonic |  64$\times$64 |  2$\times$64$\times$64 | 12,000
Experimental data |  Circular inclusion (Section 4.2.1) |  Linear, quasi-static |  56$\times$56 |  56$\times$56 | 10,000
Dual-layer phantom (Section 4.2.2) |  Linear, time-dependent |  64$\times$64 |  10$\times$64$\times$64 |  3,000
Tumor spheroids (Section 4.2.3) |  Linear, quasi-static |  256$\times$256 |  256$\times$256 | 24,000

Table 1 provides an overview of the various inverse problems we tackle
including the type of measurements, the nature of the forward models, the
dimensionality ${n_{\mathcal{X}}}$ and ${n_{\mathcal{Y}}}$ of the inferred
quantity and measurements, respectively, and the corresponding size of the
training dataset used to train the score network. More specifically, the
numerical examples have the following variety:

  1. 1.

_Synthetic and real measurements_ — In Section 4.1, we consider two inverse
problems with synthetic measurements. These problems provide the opportunity
to validate the proposed method and, for the quasi-static elastography inverse
problem in Section 4.1.1 with a relatively simple latent prior distribution,
compare it against Monte Carlo inference. In Section 4.2, we apply conditional
score-based diffusion models to solve inverse problems with measurements
coming from physical experiments. However, note that even in such cases, we
train the score network using synthetic data generated from the joint
distribution that we model.

  2. 2.

_Forward physics_ — In general, the governing mechanics model are the
following equations of equilibrium and linear elastic constitutive law

$\displaystyle\nabla\cdot\bm{\sigma}+\rho{\mathbf{b}}=\frac{\partial^{2}{\mathbf{u}}}{\partial
t^{2}},$ (20)
$\displaystyle\bm{\sigma}=2\mu\nabla^{\mathrm{s}}{\mathbf{u}}+\lambda(\nabla\cdot{\mathbf{u}})\bm{I},$
(21)

for an isotropic solid, along with appropriate boundary conditions that we
specify later. In Eqs. 20 and 21, $\bm{\sigma}$ is the Cauchy stress tensor,
$\rho$ denotes the material density, ${\mathbf{b}}$ denotes the body forces,
${\mathbf{u}}$ denotes the displacement field,
$\nabla^{\mathrm{s}}{\mathbf{u}}$ denotes the strain tensor, and $\lambda$ and
$\mu$ are the first and second Lamé parameters, respectively. For the inverse
problems in Sections 4.1.1, 4.2.1 and 4.2.3, which involve quasi-static
conditions, Eq. 20 simplifies to

$\nabla\cdot\bm{\sigma}=0,$ (22)

in the absence of body forces. However, the inverse problem in Section 4.2.2
requires the equilibrium equation Eq. 20 to account for inertia and time-
dependent body forces. For the inverse problem in Section 4.1.2, the forward
physics relating the complex-valued displacement field with the material shear
modulus field, $\mu$, is modeled using the Helmholtz equation:

$-\omega^{2}\left(u_{\mathrm{R}}+iu_{\mathrm{I}}\right)-\nabla\cdot\left[\frac{\mu}{\rho}\left(1+i\alpha\omega\right)\nabla\cdot\left(u_{\mathrm{R}}+iu_{\mathrm{I}}\right)\right]=0,$
(23)

where $\omega$ is the frequency of excitation, $u_{\mathrm{R}}$ and
$u_{\mathrm{I}}$ denote the real and imaginary components of the displacement,
$\mu$ denotes the shear modulus field, $\alpha$ is the wave dissipation
coefficient, and $\rho$ is the material density.

  3. 3.

_Measurement noise and modeling errors_ — In Sections 4.1.1, 4.1.2 and 4.2.1,
we model the measurement noise as homoscedastic zero-mean Gaussian random
variables. The measurement noise is non-Gaussian and non-additive in the
inverse problems that we consider in Sections 4.2.2 and 4.2.3. This shows that
the proposed method can handle complex measurement noise models. Additionally,
the forward model in Section 4.2.2 contains random variables that account for
the modeling error, which further highlights the flexibility of the proposed
method in accommodating complex forward models.

  4. 4.

_Black-box forward models_ — In addition to the complex noise, we use finite
element models to simulate the response of a specimen. These models are
implemented using suitable finite element packages that are used as a black-
box.

  5. 5.

_Types of measurements_ — For the inverse problems in Sections 4.1.1, 4.2.1
and 4.2.3, the underlying physics is quasi-static and we consider only the
vertical deformations. Therefore, the measurements comprise an instantaneous
displacement field after deformation (as in Sections 4.1.1 and 4.2.1) or a
quantity derived from the relative deformations of the specimen prior to and
after undergoing deformation (as in Section 4.2.3). The measurements are akin
to an image with a single channel in these cases. For the inverse problem in
Section 4.1.2, the measurement is a discretized complex-valued displacement
field, and its real and imaginary components are treated as two separate
channels of an image. Similarly, for the inverse problem in Section 4.2.2,
each instantaneous displacement field is treated as a separate channel,
leading to a total of 10 channels.

  6. 6.

_High dimensionality_ — All the inverse problems involve the inference of a
high-dimensional quantity as a result of discretizing the Young’s modulus
field over a cartesian grid. In this regard, the inverse problem in Section
4.2.3 where ${n_{\mathcal{X}}}=$256$\times$256 is the most challenging.

##### Details of implementation

We implement the proposed method on PyTorch [59]. Before training the score
network, we min-max normalize the training dataset to the range $[0,1]$. We
perform this normalization independently for the realizations of
${\mathbf{X}}$ and ${\mathbf{Y}}$ in the training dataset. For the synthetic
problems in Section 4.1, we present all results in normalized units. For the
inverse problems where we have experimental measurements, we present all
results in physical units after inverting the normalization. Additionally, for
the tumor spheroid application in Section 4.2.3, we take the logarithm of the
Young’s modulus field relative to the Young’s modulus of the surrounding media
(hydrogel) before normalizing the realizations of ${\mathbf{X}}$ in the
training dataset. This helps improve the scaling of the training data because
the Young’s modulus may vary by many orders of magnitude between the cells
within a tumor spheroid and the surrounding media. We list the training and
sampling hyper-parameters for the various inverse problems in Table 2.

Table 2: Training and sampling hyper-parameters for the conditional diffusion models Hyper-parameters | Inverse problem
---|---
Synthetic data | Experimental data
|  Circular inclusion (§ 4.1.1) |  Optic nerve head (§ 4.1.2) |  Circular inclusion (§ 4.2.1) |  Dual-layer phantom (§ 4.2.2) |  Tumor spheroids (§ 4.2.3)
First noise level $\sigma_{1}$ | 20 | 50 | 50 | 50 | 120
Final noise level $\sigma_{L}$ | 0.01 | 0.01 | 0.01 | 0.01 | 0.01
No. of noise levels $L$ | 128 | 256 | 256 | 256 | 680
No. of Langevin steps $T$ | 5 | 5 | 5 | 5 | 5
Step size $\epsilon$ |  5.7$\times$10$-$6 |  5.7$\times$10$-$6 |  5.7$\times$10$-$6 |  5.7$\times$10$-$6 |  5.7$\times$10$-$6
Learning rate | 0.0001 | 0.0001 | 0.0001 | 0.0001 | 0.0001
Batch size | 128 | 64 | 128 | 128 | 32
Total training epochs | 100,000 | 300,000 | 400,000 | 300,000 | 300,000

###  4.1 Experiments on synthetic data

In this section, we apply the conditional score-based diffusion models (cSDMs)
to two inverse problems with synthetic data. These examples help validate the
proposed approach by comparing the solutions with Monte Carlo simulation-based
inference.

####  4.1.1 Quasi-static elastography

The first synthetic example concerns quasi-static elastography, a medical
imaging technique in which tissues undergoing deformations under an external
load are imaged using ultrasound [23, 60]. The ultrasound images are used to
determine the displacements of the heterogeneous specimen, and the goal is to
infer the spatially varying shear modulus. This study features a tissue
phantom embedded with a stiffer circular inclusion within a softer surrounding
matrix, such as the realizations shown in the first row of Fig. 2. In
elastography applications, ultrasound can accurately measure the displacement
component along the transducer axis, which in this case is the vertical
direction. Thus, the inverse problem comprises inferring the spatial
distribution of the shear modulus from noisy measurements of vertical
displacements of the phantom. The mechanics is governed by the linear elastic
model Eq. 22 and 21. Moreover, we assume that the specimen is incompressible
and in plane stress condition.

Figure 2: Five typical realizations of ${\mathbf{X}}$ and ${\mathbf{Y}}$
sampled from the joint distribution for the synthetic quasi-static
elastography example and belong to the training dataset. The first row shows
the shear modulus field, while the second row shows the corresponding noisy
vertical displacement field used as measurements. In this figure, the standard
deviation of the Gaussian noise $\sigma_{\eta}=0.025$. Moreover, realizations
of ${\mathbf{X}}$ are ${\mathbf{Y}}$ are min-max normalized to $[0,1]$ scale

In this example, the dimensions of the specimen is 1$\times$1 cm2. To simulate
compression, we fix the bottom edge and impose a downwards vertical
displacement of 0.01 cm on the top edge. The bottom edge is traction-free in
the horizontal direction, whereas the left and right edges are kept traction-
free in both directions. Additionally, we pin the bottom left corner to
prevent rigid body motion. A realization of $\bm{X}$ corresponds to the
spatial distribution map of the shear modulus of the specimen discretized over
a 56$\times$56 Cartesian grid. We formulate a parametric prior that represents
a uniformly stiff circular inclusion set against a homogeneous background: the
coordinates of the inclusion’s center are uniformly distributed over the
domain such that the inclusion fits within the boundaries. We intentionally
work with a simple prior, consisting of only two random variables, to
facilitate Monte Carlo simulations (MCS). See Table 3 for details regarding
the prior distribution for this example. We fix the inclusion’s radius to be
0.12 cm and choose the shear modulus of the inclusion and the background to be
1.5 and 0.1 kPa, respectively.

Table 3: Random variables comprising the parametric prior for the synthetic quasi-static elastography problem. Note, $\mathcal{U}(a,b)$ denotes an uniform random variable between $a$ and $b$ Random variable | Distribution
---|---
Distance of the inclusion’s center from the left edge (mm) | $\mathcal{U}(0.2,0.8)$
Distance of the inclusion’s center from the bottom edge (mm) | $\mathcal{U}(0.2,0.8)$

The training dataset comprises 10,000 pairwise realizations of ${\mathbf{X}}$
and ${\mathbf{Y}}$. For a given realization ${\mathbf{x}}$ of ${\mathbf{X}}$,
we generate the corresponding measurement ${\mathbf{y}}$, which is a
realization of ${\mathbf{Y}}$, by utilizing ${\mathbf{x}}$ in a finite element
model consisting of linear triangular elements. We use FEniCS [61] to conduct
the finite element analysis. We discretize the measurement over the same
56$\times$56 Cartesian grid as the shear modulus field. Therefore,
${n_{\mathcal{X}}}={n_{\mathcal{Y}}}=$ 56$\times$56 in this example.
Subsequently, we add measurement noise — homoscedastic Gaussian noise with
standard deviation $\sigma_{\eta}\times u_{\mathrm{max}}$ truncated to within
the $[-3\sigma_{\eta}\times u_{\mathrm{max}},3\sigma_{\eta}\times
u_{\mathrm{max}}]$ range, where $u_{\mathrm{max}}$ is the maximum vertical
displacement across all specimens in the training dataset. Fig. 2 illustrates
five randomly selected realizations from the training dataset. Note that in
the training data, both ${\mathbf{X}}$ and ${\mathbf{Y}}$, are min-max
normalized to the range $[0,1]$. We work with three values of
$\sigma_{\eta}\times u_{\mathrm{max}}$ corresponding to varying intensity of
measurement noise — 2.5%, 5%, and 10% of $u_{\mathrm{max}}$ — to create three
training data sets, which we then use to train the score network. Note, we
train a separate score network for each level of measurement noise. Table 2
lists the training and sampling hyper-parameters we use for this problem.
After training the score network, we employ Annealed Langevin dynamics, as
described in Section 3.2, to sample the posterior distribution.

We also estimate the posterior statistics using Monte Carlo simulation (MCS),
which we treat as the true statistics, to validate the results obtained using
the cSDMs. We use a sample size equal to 500,000 for the Monte Carlo
simulations, which is significantly more than the size of the training
dataset. Note that the measurement noise is additive in this case and the
likelihood function can be used to obtain the importance weights, which we
subsequently self-normalize. We compare the posterior statistics estimated
using cSDM and MCS on two test samples that were not part of the training
dataset; Fig. 3 shows the results. In Fig. 3(a) and (b), column 1 shows the
test phantom while column 2 shows the corresponding measurement for three
different noise levels. Note that the ground truth is not readily discernible
through visual inspection of the measurements in column 2 of Fig. 3.

(a) Test sample 1(b) Test sample 2 Figure 3: Posterior statistics – pixel-wise
posterior mean and standard deviation – estimated using the conditional
diffusion models (cSDMs) and Monte Carlo simulation (MCS) for the synthetic
quasi-static elastography example on two test samples. In subfigures (a) and
(b), the three rows show the results corresponding to three different levels
of measurement noise $\sigma_{\eta}$ but for the same ground truth. All values
are normalized to $[0,1]$ scale Table 4: Root mean squared error (RMSE)
between the pixel-wise posterior statistics estimated using cSDM and MCS, for
the test phantoms in Fig. 3 RMSE in posterior statistics
---
Measurement noise ($\sigma_{\eta}$) | Mean | Std. Dev.
Test sample 1 | Test sample 2 | Test sample 1 | Test sample 2
10.0% | 0.036 | 0.027 | 0.038 | 0.029
5.0% | 0.039 | 0.025 | 0.032 | 0.029
2.5% | 0.031 | 0.022 | 0.028 | 0.027

Fig. 3 also shows the pixel-wise posterior mean and standard deviation
estimated using the cSDM and Monte Carlo simulations. We generate 10000, 4000,
and 1000 posterior realizations to compute these statistics using the trained
score network for $\sigma_{\eta}=$10%, 5%, and 2.5%, respectively. These
sample sizes were roughly close to the effective sample size of the Monte
Carlo samples for the three cases. Fig. 3 shows that the cSDM mean
successfully locates the inclusion at all measurement noise levels, similar to
MCS. Furthermore, the posterior statistics estimated using cSDM are
qualitatively similar to those obtained using MCS, with the root mean squared
errors (RMSE) listed in Table 4. Both these methods predict larger uncertainty
around the edges of the inclusion, and the uncertainty predictably grows with
the level of measurement noise. For the lowest level of measurement noise, the
reconstructed inclusion in the posterior mean is sharp as there is very little
uncertainty, both in spread and magnitude, surrounding the inclusion. The cSDM
yields results similar to Monte Carlo inference for other test phantoms as
well. Table 5 tabulates the RMSE between the posterior statistics estimated
using cSDM and MCS across ten test phantoms. The consistently low average RMSE
values suggest that the proposed method is accurate and robust to the amount
of measurement noise.

Table 5: Root mean squared error (RMSE) between the pixel-wise posterior
statistics estimated using cSDM and MCS, averaged over ten test phantoms for
the synthetic quasi-static elastography example Average RMSE in posterior
statistics
---
Measurement noise ($\sigma_{\eta}$) | Mean | Std. Dev.
10.0% | 0.051 | 0.051
5.0% | 0.031 | 0.032
2.5% | 0.025 | 0.029

####  4.1.2 Time-harmonic elastography

In this example, we consider the inverse Helmholtz problem, which arises in
elastography applications wherein the heterogeneous mechanical properties of a
specimen must be recovered from measurements corresponding to time-harmonic
excitation. More specifically, in this example, the shear modulus of the
tissue surrounding a typical human optic nerve head (ONH) is to be inferred
from measurements of propagating shear waves. Fig. 4 shows a typical
realization from the prior and labels the various parts surrounding the ONH.
We adapt this example from a previous study by Ray et al. [45].

Figure 4: A realization of the optic nerve head (ONH) sampled from the
parametric prior distribution showing the (1) peripapillary sclera, (2) lamina
cribrosa, (3) pia matter, (4) optic nerve, and (5) retina. This figure also
appears in [45]

The forward model comprises of the Helmholtz equation, Eq. 23, where we choose
$\alpha=$ 5$\times$10$-$5 and $\rho=$ 1000 kg/m3. The physical domain of
interest is a 1.75$\times$1.75 mm2 square, but we consider a larger
circumscribing domain to allow for wave dissipation and to avoid reflection.
We pad the left edge by 2.6 mm, and the top and bottom edges by 1.75 mm. The
right edge represents a line of symmetry and is not padded. We set
$u_{\mathrm{R}}=$ 0.02 mm and $u_{\mathrm{I}}=$ 0 on the right edge, and
$u_{\mathrm{R}}=u_{\mathrm{I}}=$ 0 on all other edges of the larger domain. We
leverage a parametric prior $\mathrm{p}_{{\mathbf{X}}}$, which we sample from
to generate realizations of the spatial distribution of the shear modulus
around the ONH. The parametric prior controls the geometry of the ONH and the
spatial distribution of the shear modulus. The components of the ONH that we
model include the peripapillary sclera, lamina cribrosa, optic nerve, and the
pia matter. We take guidance from existence literature [62, 63, 64, 65, 66,
67] to model the geometry. Table 6 shows the ranges of the dimensions for
various components of the ONH geometry resulting from the prior. Similarly,
Table 7 shows the random variables that control the magnitude of the shear
modulus field. We choose the mean and variance for these random variables
based on previous studies [68, 69, 70].

Table 6: Ranges of the dimensions for components of the optic nerve head geometry resulting from the prior distribution. These ranges result from combining multiple random variables in some cases Random variable | Range
---|---
Lamina cribrosa radius (width in cross-section) (mm) | $(0.55,1.35)$
Lamina cribrosa thickness (mm) | $(0.16,0.44)$
Peripapillary sclera thickness (mm) | $(0.45,1.15)$
Peripapillary sclera and retinal radius (mm) | $(8.0,16.0)$
Retinal thickness (mm) | $(0.20,0.40)$
Optic nerve radius (at bottom of cross-section) (mm) | $(1.10,2.20)$
Pia matter thickness (mm) | $(0.06,0.10)$
Rotation of the geometry (rad) | $(-\pi/12,\pi/12)$

In this study, a realization of ${\mathbf{X}}$ is a shear modulus field that
we discretize over a 64$\times$64 Cartesian grid. We utilize realizations of
${\mathbf{X}}$ in the finite element solver FEniCS [61] to obtain noise-free
predictions of the real and imaginary components of the displacement field. To
this, we add homoscedastic Gaussian noise with a pre-specified standard
deviation to obtain a corresponding realization of ${\mathbf{Y}}$. Therefore,
the realizations of ${\mathbf{Y}}$ comprise the noisy real and imaginary
components of the displacement field, which we denote using
${\mathbf{u}}_{\mathrm{R}}$ and ${\mathbf{u}}_{\mathrm{I}}$, respectively.
Therefore, in this case, ${n_{\mathcal{X}}}=$ 64$\times$64 and
${n_{\mathcal{Y}}}=$ 2$\times$64$\times$64\. We build a training dataset
comprising $N_{\mathrm{train}}=$ 12,000 pairs of realizations of
${\mathbf{X}}$ and ${\mathbf{Y}}$. We choose the measurement noise standard
deviation to be equal to 4% of the maximum value of
${\mathbf{u}}_{\mathrm{R}}$ and ${\mathbf{u}}_{\mathrm{I}}$ across the
training dataset. As before, we min-max normalize the training data to $[0,1]$
for better conditioning. Fig. 5 shows five randomly sampled training data
points, where we observe significant variability in the geometry and the
material property within the prior distribution. Subsequently, we train the
score network using the training dataset. See Table 2 for the training hyper-
parameters for this problem.

Table 7: Random variables comprising the parametric prior for the optic nerve head (ONH). Note, $\mathcal{N}(\varpi,\varsigma^{2})$ denotes a normal random variable with mean $\varpi$ and variance $\varsigma^{2}$ Random variable | Distribution
---|---
Optic nerve shear modulus (kPa) | $\mathcal{N}(9.8,3.34^{2})^{*}$
Sclera shear modulus (kPa) | $\mathcal{N}(125,5^{2})^{*}$
Pia matter shear modulus (kPa) | $\mathcal{N}(125,50^{2})^{*}$
Retina shear modulus (kPa) | $\mathcal{N}(9.8,3.34^{2})^{*}$
Lamina cribrosa shear modulus (kPa) | $\mathcal{N}(73.1,46.9^{2})^{*}$
Background shear modulus (kPa) | $0.1$
∗The distributions are truncated to have support between $(0,2\varsigma]$
Figure 5: Five typical realizations of ${\mathbf{X}}$ and ${\mathbf{Y}}$
sampled from the joint distribution that form the training dataset for the
synthetic time-harmonic elastography example. Note that ${\mathbf{Y}}$
comprises of ${\mathbf{u}}_{\mathrm{R}}$ and ${\mathbf{u}}_{\mathrm{I}}$ in
this case. The first row shows the spatial distribution of the shear modulus
field around the optic nerve head, the second and third rows show the
corresponding noisy measurements ${\mathbf{u}}_{\mathrm{R}}$ and
${\mathbf{u}}_{\mathrm{I}}$, respectively. All values are normalized to
$[0,1]$

We sample the posterior distribution for two test cases that were not part of
the training dataset using the trained score network and Annealed Langevin
dynamics. For these cases, the first column in Fig. 6 shows the ground truth
shear modulus distribution. The second and third columns in Fig. 6 show the
corresponding noisy measurements. Note that the shear modulus distribution and
the corresponding measurements are all plotted in the normalized $[0,1]$
scale. We sample the target posterior distribution for both cases using
Annealed Langevin dynamics and generate 3,000 realizations for each case;
Table 2 lists the sampling hyper-parameters for this problem. We use the
generated samples to estimate the posterior pixel-wise mean and standard
deviation, which we plot in the fourth and fifth columns of Fig. 6,
respectively. The last column in Fig. 6 shows the absolute pixel-wise error
between the posterior mean estimated using the conditional diffusion model and
the corresponding ground truth. The RMSE between the posterior mean and the
ground truth are 0.021 and 0.032 for the two test cases, which shows that the
posterior mean provides a good reconstruction for the ground truth. We observe
relatively larger uncertainty in regions of sharp change in the elastic
modulus field, which also coincides with areas of larger reconstruction
errors.

Figure 6: Posterior statistics estimated using the conditional diffusion model
for the inverse Helmholtz problem with synthetic data. Note, the ground truth
distribution of the shear modulus, corresponding measurements and posterior
statistics are all plotted in the normalized $[0,1]$ scale

###  4.2 Applications with experimental data

In this section, we apply the cSDMs to three inverse problems with
experimental data. We train the score network on synthetic data but perform
inference using experimental data. These examples help demonstrate the utility
of the proposed approach to problems of practical relevance.

####  4.2.1 Quasi-static elastography of a circular inclusion within
homogeneous media

Similar to the example in Section 4.1.1, this example also concerns an
application of the proposed methodology to quasi-static elastography. We adapt
this example from [44]. The inverse problem at hand is that of inferring the
spatial distribution of the shear modulus of a specimen from noisy full-field
measurements of the vertical displacements.

For an elastic medium, the relation between the displacements and shear
modulus is given by equilibrium and constitutive equations Eqs. 22 and 21,
respectively. We also make plane stress and incompressibility assumptions. In
this example, the specimen is 34.608$\times$26.297 mm2, with traction free
boundaries along the left and right edges. Additionally, the top and bottom
surfaces are traction free along the horizontal direction. To simulate
compression of the specimen, the top and bottom edges are subject to 0.084 mm
and 0.392 mm vertical displacements, respectively.

In this example, a realization of ${\mathbf{X}}$ corresponds to a realization
of the spatial distribution of the shear modulus of the specimen discretized
over a 56$\times$56 Cartesian grid. The parametric prior for the spatial
distribution of the shear modulus models a stiff circular inclusion in a
uniform background of 4.7 kPa. Thus, the parametric prior consists of four
random variables, which control the coordinates of the center of the
inclusion, the radius of the inclusion, and the ratio of the shear modulus of
the inclusion with respect to the shear modulus of the background. Table 8
contains details of the parametric prior for this application [44]. We obtain
a realization of ${\mathbf{Y}}$ by propagating the corresponding realization
of ${\mathbf{X}}$ through a finite element model with linear triangular
elements, and adding homoscedastic Gaussian noise with standard deviation
equal to 0.001 mm to the predicted vertical displacement. In this manner, we
build a training dataset consisting of $N_{\mathrm{train}}=8000$ pairwise
realizations of ${\mathbf{X}}$ and ${\mathbf{Y}}$. Fig. 7 shows five data
points randomly sampled from the training dataset. Using this dataset we train
the score network; see Table 2 for more details regarding the various training
and sampling hyper-parameters.

Table 8: Random variables comprising the parametric prior distribution for ${\mathbf{X}}$ in the quasi-static elastography application Random variable | Distribution
---|---
Distance of the inclusion’s center from the left edge (mm) | $\mathcal{U}(7.1,19.2)$
Distance of the inclusion’s center from the bottom edge (mm) | $\mathcal{U}(7.1,27.6)$
Radius of the inclusion (mm) | $\mathcal{U}(3.5,7.0)$
Ratio between the inclusion’s and background’s shear modulus (mm) | $\mathcal{U}(1,8)$
Figure 7: Five typical realizations of ${\mathbf{X}}$ and ${\mathbf{Y}}$
sampled from the joint distribution that forms the training dataset for the
quasi-static elastography application. The first row shows the spatial
distribution of the shear modulus field, and the second row shows the
corresponding noisy measurements of the vertical component of the noisy
displacement field ${\mathbf{u}}$

Next, we use the trained score network and Annealed Langevin dynamics to
sample the posterior distribution and infer the spatial distribution of the
shear modulus for two test cases that were not part of the training dataset.
Fig. 8 shows the results for these tests. In Fig. 8, the first and second
columns show the true spatial distribution for the two test cases and the
corresponding full-field measurements, respectively. We use Annealed Langevin
dynamics sampling to generate 800 posterior realizations for each test case
and then utilize them to estimate the pixel-wise posterior mean and standard
deviation. We show the posterior statistics in the third and fourth columns of
Fig. 8. Finally, the last column in Fig. 8 shows the absolute pixel-wise error
between the posterior mean estimated using the conditional diffusion model and
the corresponding ground truth. The RMSE between the posterior mean and the
ground truth are 0.46 and 0.63 for the two test cases. We observe relatively
larger uncertainty around the edges of the stiff inclusion and comparatively
less uncertainty about the stiffness of the inclusion.

Figure 8: Posterior statistics estimated using the conditional diffusion model
on test samples for the quasi-static elastography application

Finally, we use the trained score network to infer the spatially varying shear
modulus of a tissue-mimicking phantom specimen [56]. The phantom consists of a
stiff inclusion embedded within a softer substrate and was manufactured using
a mixture of gelatin, agar, and oil. The experimental measurements were
obtained from ultrasound scans as the phantom was gently compressed. The left
column in Fig. 9 shows the measured vertical displacement. The second and
third columns show the pixel-wise posterior mean and standard deviation
estimate using 800 posterior samples obtained from Annealed Langevin dynamics
sampling performed with the trained score network. From the posterior mean,
the average stiffness and diameter of the inclusion are 12.94 kPa and 10.8 mm,
respectively. These are close to the corresponding experimental measurements
of 10.7 kPa and 10.0 mm, respectively. These predictions are consistent with a
previous study [44].

Figure 9: Posterior statistics estimated using the conditional diffusion model
on experimental data for the quasi-static elastography application

####  4.2.2 Time-dependent elastography of a dual-layer tissue-mimicking
phantom

This problem consists of inferring the spatial distribution of Young’s modulus
in a dual-layer tissue-mimicking phantom from sequential measurements of the
displacement field made in an ultrasound elastography experiment [57]. The
phantom is a dual-layer cylinder of 10 mm diameter. The phantom was fabricated
using 5% gelatin concentration for the top layer and 15% concentration for the
bottom layer. The top and bottom layers were 1 mm and 3 mm thick,
respectively. In the experiment, acoustic radiation force (ARF) is applied for
200 $\mu$s at sufficient depth to excite a shear wave in both layers. The ARF
pushes the phantom through the central axis of the cylinder. The vertical
component of the displacement field corresponding to the propagating shear
wave is measured with an ultrasound transducer for the next 10 ms. The
measurements are recorded at a frequency of 10 kHz and cover a two-dimensional
cross-sectional plane through the center of the specimen. Based on the
geometry of the specimen, the boundary conditions, and the spatial
distribution of material properties and the ARF force, we assume that the
problem is asymmetric about a vertical axis running through the center of the
specimen. We wish to recover the Young’s modulus distribution within the
phantom using these measurements.

Table 9: Random variables comprising the parametric prior for the time-dependent elastography application Parameter | Distribution
---|---
Young modulus of top phantom (kPa) | $\mathcal{U}(1,25)$
Young modulus of bottom phantom (kPa) | $\mathcal{U}(1,25)$
Height of the interface from bottom surface(mm) | $\mathcal{U}(2.0,3.6)$

In this case, realizations of ${\mathbf{X}}$ correspond to the spatial
distribution of the Young’s modulus discretized on a 64$\times$64 Cartesian
grid over a 3.9$\times$3.9 mm2 region. Thus, ${n_{\mathcal{X}}}=64\times 64$
in this application. We assume that each layer in the phantom is homogeneous.
Hence, the parametric prior model consists of three random variables: the
Young’s Modulus of the top and bottom layers and the location of the
interface. See Table 9 for details of the parametric prior. The forward model
consists of a finite element [71] solution to the equations of elastodynamics
with a linear elastic constitutive law, Eqs. 20 and 21, respectively. In this
example, ${\mathbf{b}}$ in Eq. 20 is the body force due to the ARF, and
$\bm{\sigma}$ is the stress tensor. Additionally, we assume axisymmetry
conditions for the specimen. Hence, we consider and plot only half of the
specimen subsequently, i.e., the left edge of the specimen in Fig. 10
corresponds to the axis of rotation. We model the ARF using a body force given
by

${\mathbf{b}}(s_{1},s_{2})=-F\exp\left\\{-\frac{s_{1}^{2}}{2l_{1}^{2}}-\frac{(s_{2}-\bar{s}_{2})^{2}}{2l_{2}^{2}}\right\\}\bm{e}_{2},$
(24)

where $s_{1}$ and $s_{2}$ are the physical Cartesian coordinates, $\bm{e}_{2}$
is a unit vector in the vertical direction, $F=-$5$\times$105 is a
multiplicative factor used to achieve the desired magnitude of displacement,
$\bar{s}_{2}=$ 2.6 mm is center of the ARF along the vertical direction, and
$l_{1}=$ 0.2 mm and $l_{2}=$ 0.6 mm denote the attenuation lengths of the ARF
along the horizontal and vertical directions, respectively. The ARF is applied
over 200 $\mu$s to replicate the experimental setup. Once the forcing is
removed, the propagation of the shear wave is modeled for 0.9 ms and ten
snapshots of the vertical displacement field are recorded on a 64$\times$64
grid at a frequency of 10 kHz. This leads to ${n_{\mathcal{Y}}}=10\times
64\times 64$.

We model errors in the forward model by introducing uncertainty in the modulus
of the two layers and the location of the interface. This is accomplished by
assigning the computed displacement field to a slightly perturbed version of
the specimen for which it was calculated. The perturbed specimen is obtained
by (a) scaling the modulus of each layer by a factor sampled from the uniform
distribution in the interval (0.8, 1.2), and (b) translating the interface
vertically by an amount sampled from a uniform distribution in the interval
($-$0.2, 0.2) mm. These model errors are included in the forward model to
account for the uncertainty in the manufacturing process. Finally, to obtain
the noisy measurements, a Gaussian blurring kernel of width equal to 6 pixels
is applied, and homoscedastic Gaussian noise with standard deviation equal to
5$\times$10$-$5 mm is added. In this manner, a total of 3,000 training data
points are generated. Fig. 10 shows two typical realizations of ${\mathbf{X}}$
and ${\mathbf{Y}}$ sampled from the joint distribution that are part of the
training dataset. We use this dataset to train the score network and
subsequently sample the posterior distribution using Annealed Langevin
dynamics.

Fig. 11 shows the results for two test samples where we use synthetic
measurements to infer the spatial distribution of the Young’s moduli. The
first column shows the ground truth for the Young’s modulus distribution of
the phantom. The second to fourth columns show three posterior realizations.
The fifth and sixth columns show the pixel-wise posterior mean and standard
deviation estimated using the conditional score-based diffusion model with 500
posterior realizations. The posterior mean provides a good estimate for the
ground truth both in terms of the Young’s modulus for each layer and the
location of the interface. Table 10 quantitatively compares the estimated
posterior mean values of the Young’s modulus and the location of the interface
to their corresponding true values.

Figure 10: Two samples from the training dataset are shown. The Young’s modulus distribution is shown on the left and the displacements during the $0.9$ ms after the forcing is removed are shown on the grid to the right Figure 11: Posterior statistics estimated using the conditional diffusion model on test data for the time-dependent elastography application Table 10: Posterior mean estimates of the Young’s modulus of the top and bottom layers of the dual-layer phantom, and the location of the interface obtained using cSDM compared to their true value for two test samples Test sample |  Young’s modulus of top layer (kPa) |  Young’s modulus of bottom layer (kPa) |  Location of interface (mm)
---|---|---|---
True value |  Posterior mean |  True value |  Posterior mean |  True value |  Posterior mean
1 | 14.74 | 14.20 | 6.66 | 6.01 | 2.48 | 2.57
2 | 21.81 | 23.42 | 9.70 | 8.48 | 2.11 | 2.05
Figure 12: Ultrasound measurements of the displacement in $\mu$m of a dual
layer phantom during an interval of $0.9$ ms. Data collected by Lu et al. [57]
Figure 13: Posterior realizations (first three columns) and statistics
estimated using the conditional diffusion model on experimental data for the
ultrasound elastography application

Fig. 13 shows the results of the inference for the experimental data collected
by Lu et al. [57] and shown in Fig. 12. This figure shows three typical
samples from the posterior distribution and the pixel-wise mean and standard
deviation computed using 500 realizations from the posterior distribution. The
posterior mean estimate for the Young’s modulus of the top and bottom layers
is 6.98 kPa and 19.73 kPa, respectively. Uniaxial quasi-static mechanical
tests performed on homogenous phantoms with the same gelatin concentrations as
the two layers yielded values of 2.1 $\pm$ 0.2 kPa and 20.3 $\pm$ 0.4 kPa.
Thus, the cSDM estimates the modulus of the lower layer with good accuracy.
However, there is a significant error in the estimate for the modulus of the
top layer, which may be because this layer is relatively thin and, therefore,
its contribution to the measured displacement is unlikely to be significant.
The phantom was fabricated such that the interface would be at a height of 3.0
mm. The cSDM estimate for the mean interface height is 3.21 mm with an average
standard deviation of 0.10 mm; the true location of the interface is close to
two standard deviations of the estimated posterior mean.

####  4.2.3 Optical coherence elastography of tumor spheroids

Optical coherence tomography (OCT) is a technique that leverages optics,
analogous to acoustics in ultrasound imaging, to produce high-resolution
tomograms of tissue. A functional extension of OCT that enables imaging based
on mechanical contrast is optical coherence elastography (OCE) [72], which
involves capturing the phase data alongside intensity measurements of the
interference of light backscattered from the specimen with light reflected
from a reference mirror at varied axial depths. One variant of OCE is phase-
sensitive compression OCE [73], wherein the phase shifts observed between
B-scans (cross-sectional images) of the specimen as it undergoes compressive
loading quantify the specimen’s deformation along the optical axis. The
mechanical properties of the specimen can be inferred from these phase shifts
because they encode within them information about the displacement component
along the optical axis. In this study, we consider a high resolution variant
of phase-sensitive compression OCE, termed mechano-microscopy [74, 75], which
utilizes a high-resolution variant of OCT known as optical coherence
microscopy (OCM). A conditional score-based diffusion model is used to solve
the inverse problem of determining elasticity profiles, specifically the
spatial distribution of the Young’s modulus ${\mathbf{X}}$ of tumor spheroids
from noisy phase difference measurements ${\mathbf{Y}}$. The data for this
experiment were obtained from measurements and datasets originally compiled by
Foo et al. [58].

We provide a brief summary of the experimental setup here (Fig. 14) and the
reader is referred to [58] for more details. In the experiment, a custom-built
spectral-domain OCM system in dual-arm configuration images tumor spheroids
embedded in hydrogel specimens from the top. Following the objective lens, the
sample arm contains, from top to bottom: a rigid glass coverslip acting as an
imaging window, a 300 $\mu$m thick sample containing tumor spheroids, a 1 mm
thick compliant silicone layer, and a piezoelectric actuator to apply
compressive displacement. The purpose of the compliant layer is to provide
surface stress measurements for mechano-microscopy [74]; we do not utilize
this data in this study. Two volumetric images at different compression levels
were acquired at each location to produce a volumetric phase difference image.
The plane corresponding to the central cross-section of the tumor spheroid was
identified in the volumetric image and was cropped to a 256$\times$256 $\mu$m2
window around the spheroid. This two-dimensional phase difference image was
used as measurement in solving the inverse problem.

Figure 14: Schematic showing the experimental setup in the tumor spheroid
application

The specimens were prepared by encapsulating MCF7 (non-metastatic breast
cancer) cells within a hydrogel matrix composed of lyophilized gelatin
methacryloyl (GelMA). The cells were cultured until they formed spheroids.
Each tumor spheroid contains multiple cells, and each cell domain can be
broadly classified into two components: the cytoplasm and the nucleus (see
Fig. 15). We take inspiration from this morphology to develop the parametric
prior for this study. We refer interested readers to [76] for typical three-
dimensional immuno-fluorescence images of tumor spheroids taken at different
time points.

The synthetic samples, which we generate from the parametric prior
distribution, represent a spheroid within a hydrogel on a 256$\times$256
$\mu$m2 domain (first column in Fig. 16). We model the spheroid as an ellipse,
and straight lines subdivide the spheroid into cells. Each cell contains a
circle representing a nucleus, surrounded by the cytoplasm of the cell. We
treat the following parameters as truncated normal random variables: spheroid
width, spheroid aspect ratio, spheroid position, number of nuclei per unit
area, nuclei radii, nuclei elasticity, cytoplasm elasticity, and hydrogel
elasticity. Table 11 provides details on these parameter’s distributions,
which we derive from preliminary analyses of spheroid OCM images. We determine
the nuclei locations (defined by their center coordinates) using a custom
algorithm that iteratively accumulates nuclei coordinates by creating a list
of valid coordinates for the current nucleus to be positioned, randomly
selecting a coordinate from this list for the current nucleus, removing the
coordinates occupied by the current nucleus from the list of valid coordinates
for the next nucleus, and repeating until the required number of nuclei for
the spheroid is reached or there are no more valid coordinates for the next
nucleus. We constrain the sizes and positions of the nuclei to ensure a
minimum spacing of 2 $\mu$m between nuclei and between nuclei and the ellipse
boundary. We subsequently determine the cell boundaries using a power diagram,
also known as a Laguerre–Voronoi diagram [77].

Table 11: Random variables comprising the parametric prior for the tumor
spheroid application

Parameter | Variable | Minimum | Maximum
---|---|---|---
Spheroid width ($\mu$m) | $\mathcal{N}(105,8.33)$ | 80 | 130
Spheroid height to width ratio | $\mathcal{N}(0.90,0.0167)$ | 0.85 | 0.95
Spheroid (center) horizontal position ($\mu$m) | $\mathcal{N}(0,13.33)$ | -40 | 40
Spheroid (center) vertical position ($\mu$m) | $\mathcal{N}(0,13.33)$ | -40 | 40
Number of nuclei per square area ($\text{mm}^{-2}$) | $\mathcal{N}(4750,250)$ | 4000 | 5500
Radius of nuclei ($\mu$m) | $\mathcal{N}(7.5,0.5)$ | 6 | 9
Young’s modulus of nuclei (kPa) | $\mathcal{N}(32.5,5.833)$ | 15 | 50
Young’s modulus of cytoplasm (kPa) | $\mathcal{N}(10,1.67)$ | 5 | 15
Young’s modulus of hydrogel (kPa) | $\mathcal{N}(7.5,0.833)$ | 5 | 10
Displacement applied to top edge ($\mu$m) | $\mathcal{N}(-2.2,0.6)$ | $-0.4$ | $-4$

*In this example, the normal distributions were truncated at three standard deviations from the mean, with truncation limits explicitly defined by the specified minimum and maximum values.

Fig. 15 shows the steps for generating realizations of ${\mathbf{Y}}$ from
realizations of ${\mathbf{X}}$. Given a realization ${\mathbf{x}}$ of the
Young’s modulus field, we simulate uniaxial compression on the specimen via
finite element analysis (FEA) in ABAQUS [71] to obtain the vertical
displacement field (third column in Fig. 15). To mitigate edge effects, we
initially expand each sample to an 800$\times$800 $\mu\text{m}^{2}$ domain by
padding with hydrogel on all four sides. We displace the top edge downward
while the bottom-left corner remains fixed, and the bottom edge is constrained
vertically. We sample the displacement on the top edge from a truncated normal
distribution; see Table 11. The magnitude of the displacement of the top edge
ensures that the corresponding strain values are within the order of a few
millistrains, consistent with observations from mechano-microscopy
experiments. For this analysis, given the low levels of strain involved, we
assume that the specimen is linear and elastic (Eq. 21), under plane strain
conditions (Eq. 22), and nearly incompressible with Poisson’s ratio uniformly
set to 0.49. We develop a CAD model for the specimen (second column in Fig.
15), which we subsequently use to develop the FEA mesh geometry consisting of
quadrilateral elements for the main region and triangular elements for the
padded regions. Once the FEA is complete, we obtain 2D arrays of the Young’s
modulus $E$ and vertical displacement $u_{2}$ by interpolating the main
256$\times$256 $\mu\text{m}^{2}$ domain window to a 256$\times$256 Cartesian
grid. We convert the displacement fields to phase difference fields
$\Delta\phi$ using the following equation:

$\Delta\phi=\frac{4\pi n_{r}u_{2}}{\lambda_{0}}\;,$ (25)

where $n_{r}=1.4$ is the refractive index of the hydrogel, and
$\lambda_{0}=800$ nm is the central wavelength of the OCM light source.

Figure 15: Workflow for generating the training data in the tumor spheroid
application

Subsequently, we add non-Gaussian measurement noise $\phi_{\eta}$ to the phase
difference field to obtain a noisy measurement image (fourth column in Fig.
15). The noise in the phase difference varies with depth and is modeled to
follow a probability law characterized by the following density function [58]:

$\text{p}_{\Phi_{\eta}}(\phi_{\eta})=\frac{1}{2\pi}\exp{\left(-\frac{k^{2}}{2}\right)}+\sqrt{\frac{1}{8\pi}}k\exp{\left(-\frac{k^{2}}{2}\sin^{2}{\phi_{\eta}}\right)}\cos{\phi_{\eta}}\left[1+\operatorname{erf}\left(\frac{k\cos{\phi_{\eta}}}{\sqrt{2}}\right)\right],$
(26)

where $k$ is related to the variance of this distribution. The variance varies
with depth along the optical axis in a B-scan due to optical attenuation and
the confocal function of the objective lens, and was determined as a function
of the depth by curve fitting based on experimental calibration conducted with
hydrogel specimens [58]. After the addition of noise, the phase difference is
wrapped to ensure it remains within the range $(-\pi,\pi]$. We obtain the
final wrapped phase difference image $\Delta\phi_{\rm w}$ using the following
formula

$\Delta\phi_{\rm
w}=\arg\left(\exp\left[i\left(\Delta\phi+\phi_{\eta}\right)\right]\right).$
(27)

This wrapped and noisy phase difference field is the realization
${\mathbf{y}}$ of the measurements corresponding to the realization
${\mathbf{x}}$ of the Young’s modulus field.

In this manner, we initially synthesize a synthetic dataset comprising 5000
pairs of $\bm{X}$ and $\bm{Y}$ realizations. We augment this dataset sixfold
by shifting the 256$\times$256 $\mu\text{m}^{2}$ image window 20 $\mu$m in
each horizontal and vertical direction (left, right, up, and down), as well as
applying a horizontal flip to each image. We divide the resulting dataset into
a training set of 24000 images and a test set of 6000 images. As a
preprocessing step, we rescale the elasticity images to enhance contrast and
facilitate better training convergence. The rescaling involves the following
transformation

$E^{\prime}=\log_{10}\left(\frac{E}{E_{\text{hydrogel}}}\right),$ (28)

where $E^{\prime}$ is the rescaled Young’s modulus. This transformation means
that it will only be possible to recover the Young’s modulus field up to the
normalization constant $E_{\text{hydrogel}}$, which is the Young’s modulus of
the hydrogel. Fig. 16 illustrates five randomly selected samples from the
dataset. We utilize the training data to train the score network. After
training the score network, we employ Annealed Langevin dynamics to sample
elasticity images from the posterior. See Table 2 for additional information
regarding the training and sampling hyperparameters.

Figure 16: Five typical realizations of $\bm{X}$ and $\bm{Y}$ sampled from the
joint distribution that form the training dataset for the tumor spheroid
experiment. In the first row, we have the log-normalized Young’s module
fields, and in the second row, the corresponding noisy measurements are
displayed. All values are normalized to $[0,1]$ scale

Fig. 17 shows realizations from the posterior, and its statistics, estimated
using the conditional score-based diffusion model on both synthetic and
experimental data. The first two columns correspond to samples from the test
dataset previously unseen by the model, while columns three and four depict
predictions on experimental measurements. The first row shows the measurements
used for sampling using the trained score network. Rows two to four show three
realizations sampled by the model for the corresponding measurements. Rows
five and six show the pixel-wise posterior mean and standard deviation,
respectively, estimated using 200 realizations. For the synthetic samples, we
compare the estimated posterior means for the Young’s modulus distribution
with the corresponding ground truth. In the case of experimental data,
comparisons were made against OCM images that were co-registered during the
acquisition of the phase difference images. Note that all images are presented
at their true scales. Recall that the recovery of the Young’s modulus is
possible up to a multiplicative constant: the hydrogel modulus
$E_{\text{hydrogel}}$ in this case. We assume $E_{\text{hydrogel}}=$ 10 kPa
for the test samples. For the experimental samples, $E_{\text{hydrogel}}$ was
estimated from mechano-microscopy measurements [58].

Figure 17: Posterior statistics estimated for selected samples using the cSDM
for the tumor spheroid problem. The field of view in each image is
256$\times$256 $\mu\text{m}^{2}$

From the results for the synthetic data, we note that each realization from
the posterior distribution correctly captures the volume and the location of
the spheroid. Further, there are significant differences between each
realization in the modulus distribution predicted within each spheroid. We can
attribute these differences to the ill-posed nature of the inverse problem
[78, 79], which stipulates that there are multiple modulus distributions that
are consistent with a given noisy measured displacement field. The differences
in the modulus field within the spheroid lead to a fuzzy mean modulus field.
Remarkably, when we compare this field to the ground truth, we observe that
the location of several nuclei, which correspond to regions with elevated
stiffness, is captured accurately. The uncertainty in the modulus field within
the spheroid directly translates to large values of standard deviation within
the spheroid. We also observe that regions of the largest standard deviation
(uncertainty) are associated with the boundaries of the nuclei.

From the results for the experimental data, we observe that the measured phase
field is more heterogeneous and noisy than its synthetic counterpart. Thus, we
anticipate more variability and uncertainty in the predictions, which is
precisely what we observe when examining the three realizations from the
posterior distribution. This translates to a fuzzy image for the mean
stiffness field and large standard deviation values (when compared with
synthetic data) within the spheroid. The underlying ground truth image is
unavailable in the experimental case. Therefore, we cannot quantitatively
evaluate the accuracy of the modulus field. However, from the OCM images (last
row of Fig. 17), we can infer the location and the size of the spheroid, and
when we compare this with the mean modulus image we conclude that the location
and the size of the spheroid have been captured accurately.

We conclude this section by noting that although the initial tumor spheroid
results are encouraging, they point to several extensions that are currently
underway. These include improving the forward model for measurement noise so
that the synthetic measurements resemble their experimental counterparts more
closely and performing inference on tissue-mimicking phantoms of tumor
spheroids for a quantitative comparison between the inferred results and the
ground truth.

##  5 Conclusions, outlook, and extensions

Large-scale Bayesian inference using black-box models is challenging. In this
work, we explore the use of conditional score-based diffusion models for
large-scale Bayesian inference. Fundamental to the score-based diffusion
models is the score network that can be trained using realizations from the
joint distribution between the inferred quantity and the measurements. Thus,
we only need to simulate from black-box models to create the training dataset.
The trained score network can be used to sample the target posterior
distribution using Annealed Langevin dynamics. The score network learns the
score function of the target posterior ‘averaged’ over the joint distribution.
Therefore, there is no need to re-train the score function when a new
measurement is made. This makes conditional score-based diffusion models
particularly suitable for high-throughput applications.

We apply conditional score-based diffusion models to solve several inverse
problems in mechanics. We validate the proposed approach on two synthetic
examples; for these problems the quality of the inference is similar to Monte
Carlo simulation. Further, we demonstrate the utility of the proposed approach
on three applications where experimental measurements are available. Overall,
we show that conditional score-based diffusion models can tackle high-
dimensional inverse problems involving complex prior distributions, different
types of forward physics simulated using black-box models, complex measurement
noise models and modeling errors, and various measurement modalities.

Although the results in this manuscript are promising, we identify several
interesting future research directions that will help with the adoption of
this technology as an inference toolbox. First, the primary computational
bottleneck in applying conditional score-based diffusion models, and any
conditional generative algorithm for that matter, will be the cost of curating
the training dataset, which requires multiple simulations of the black-box
forward model. Thus, there is a need for the development of data-efficient
conditional diffusion models. Moreover, training using multi-fidelity datasets
promises to be another interesting research direction. Further, we consider
full-field measurements in this work. Extensions to sparse measurements using
masking-based training strategies [80] is another extension that will increase
the utility of the proposed method. Finally, applying conditional score-based
diffusion models to physics-based inverse problems that arise in different
science and engineering domains will be interesting, as it will reveal
insights into the behavior of conditional generative algorithms of this type.

##  6 Acknowledgments

The authors acknowledge support from ARO, United States, USA grant
W911NF2010050, and from the National Institutes of Health grant
1R01EY032229-01. BFK acknowledges funding from the Australian Research
Council, the NAWA Chair programme of the Polish National Agency for Academic
Exchange and from the National Science Centre, Poland. The authors acknowledge
the following colleagues for their contribution to the results presented in
Section 4.2.3: Alireza Mowla, Matt Hepburn, Jiayue Li, Danielle Vahala,
Sebastian Amos, Liisa Hirvonen, Samuel Maher, and Yu Suk Choi. The authors
also acknowledge the Center for Advanced Research Computing (CARC,
carc.usc.edu) at the University of Southern California for providing computing
resources that have contributed to the research results reported within this
publication.

## References

  * Murphy [2023] K. P. Murphy, Probabilistic Machine Learning: Advanced Topics, MIT Press, 2023. URL: http://probml.github.io/book2.
  * Dhariwal and Nichol [2021] P. Dhariwal, A. Nichol,  Diffusion models beat GANs on image synthesis,  Advances in Neural Information Processing Systems 34 (2021) 8780–8794.
  * Achiam et al. [2023] J. Achiam, S. Adler, S. Agarwal, L. Ahmad, I. Akkaya, F. L. Aleman, D. Almeida, J. Altenschmidt, S. Altman, S. Anadkat, et al.,  GPT-4 Technical Report,  arXiv preprint arXiv:2303.08774 (2023).
  * Ramesh et al. [2022] A. Ramesh, P. Dhariwal, A. Nichol, C. Chu, M. Chen,  Hierarchical text-conditional image generation with CLIP latents,  arXiv preprint arXiv:2204.06125 1 (2022) 3.
  * Esser et al. [2024] P. Esser, S. Kulal, A. Blattmann, R. Entezari, J. Müller, H. Saini, Y. Levi, D. Lorenz, A. Sauer, F. Boesel, et al.,  Scaling rectified flow transformers for high-resolution image synthesis,  arXiv preprint arXiv:2403.03206 (2024).
  * Yang et al. [2023] L. Yang, Z. Zhang, Y. Song, S. Hong, R. Xu, Y. Zhao, W. Zhang, B. Cui, M.-H. Yang,  Diffusion models: A comprehensive survey of methods and applications,  ACM Computing Surveys 56 (2023) 1–39.
  * Li et al. [2023] Y. Li, K. Zhou, W. X. Zhao, J.-R. Wen,  Diffusion models for non-autoregressive text generation: A survey,  arXiv preprint arXiv:2303.06574 (2023).
  * Croitoru et al. [2023] F.-A. Croitoru, V. Hondru, R. T. Ionescu, M. Shah,  Diffusion models in vision: A survey,  IEEE Transactions on Pattern Analysis and Machine Intelligence (2023).
  * Kazerouni et al. [2022] A. Kazerouni, E. K. Aghdam, M. Heidari, R. Azad, M. Fayyaz, I. Hacihaliloglu, D. Merhof,  Diffusion models for medical image analysis: A comprehensive survey,  arXiv preprint arXiv:2211.07804 (2022).
  * Guo et al. [2024] Z. Guo, J. Liu, Y. Wang, M. Chen, D. Wang, D. Xu, J. Cheng,  Diffusion models in bioinformatics and computational biology,  Nature Reviews Bioengineering 2 (2024) 136–154.
  * Düreth et al. [2023] C. Düreth, P. Seibert, D. Rücker, S. Handford, M. Kästner, M. Gude,  Conditional diffusion-based microstructure reconstruction,  Materials Today Communications 35 (2023) 105608.
  * Vlassis and Sun [2023] N. N. Vlassis, W. Sun,  Denoising diffusion algorithm for inverse design of microstructures with fine-tuned nonlinear material properties,  Computer Methods in Applied Mechanics and Engineering 413 (2023) 116126\.
  * Bastek and Kochmann [2023] J.-H. Bastek, D. M. Kochmann,  Inverse design of nonlinear mechanical metamaterials via video denoising diffusion models,  Nature Machine Intelligence 5 (2023) 1466–1475.
  * Mazé and Ahmed [2023] F. Mazé, F. Ahmed,  Diffusion models beat GANs on topology optimization,  in: Proceedings of the AAAI conference on artificial intelligence, volume 37, 2023, pp. 9108–9116.
  * Li et al. [2023] T. Li, A. S. Lanotte, M. Buzzicotti, F. Bonaccorso, L. Biferale,  Multi-scale reconstruction of turbulent rotating flows with generative diffusion models,  Atmosphere 15 (2023) 60.
  * Shu et al. [2023] D. Shu, Z. Li, A. B. Farimani,  A physics-informed diffusion model for high-fidelity flow field reconstruction,  Journal of Computational Physics 478 (2023) 111972.
  * Dimakis et al. [2022] A. G. Dimakis, A. Bora, D. Van Veen, A. Jalal, S. Vishwanath, E. Price,  Deep generative models and inverse problems,  Mathematical Aspects of Deep Learning (2022) 400.
  * Avril et al. [2008] S. Avril, M. Bonnet, A.-S. Bretelle, M. Grédiac, F. Hild, P. Ienny, F. Latourte, D. Lemosse, S. Pagano, E. Pagnacco, et al.,  Overview of identification methods of mechanical parameters based on full-field measurements,  Experimental Mechanics 48 (2008) 381–402.
  * Yang et al. [2017] Y. Yang, P. Sun, S. Nagarajaiah, S. M. Bachilo, R. B. Weisman,  Full-field, high-spatial-resolution detection of local structural damage from low-resolution random strain field measurements,  Journal of Sound and Vibration 399 (2017) 75–85.
  * Dong et al. [2020] C.-Z. Dong, O. Celik, F. N. Catbas, E. J. O’Brien, S. Taylor,  Structural displacement monitoring using deep learning-based full field optical flow methods,  Structure and Infrastructure Engineering 16 (2020) 51–71.
  * Sarvazyan [1998] A. Sarvazyan,  Mechanical imaging: A new technology for medical diagnostics,  International journal of medical informatics 49 (1998) 195–216.
  * Patel et al. [2019] D. Patel, R. Tibrewala, A. Vega, L. Dong, N. Hugenberg, A. A. Oberai,  Circumventing the solution of inverse problems in mechanics through deep learning: Application to elasticity imaging,  Computer Methods in Applied Mechanics and Engineering 353 (2019) 448–466.
  * Barbone and Oberai [2007] P. E. Barbone, A. A. Oberai,  Elastic modulus imaging: some exact solutions of the compressible elastography inverse problem,  Physics in Medicine & Biology 52 (2007) 1577.
  * Song et al. [2021] Y. Song, L. Shen, L. Xing, S. Ermon,  Solving inverse problems in medical imaging with score-based generative models,  arXiv preprint arXiv:2111.08005 (2021).
  * Song et al. [2022] J. Song, A. Vahdat, M. Mardani, J. Kautz,  Pseudoinverse-guided diffusion models for inverse problems,  in: International Conference on Learning Representations, 2022.
  * Chung et al. [2022a] H. Chung, B. Sim, D. Ryu, J. C. Ye,  Improving diffusion models for inverse problems using manifold constraints,  Advances in Neural Information Processing Systems 35 (2022a) 25683–25696.
  * Chung et al. [2022b] H. Chung, J. Kim, M. T. Mccann, M. L. Klasky, J. C. Ye,  Diffusion posterior sampling for general noisy inverse problems,  arXiv preprint arXiv:2209.14687 (2022b).
  * Chung et al. [2022c] H. Chung, B. Sim, J. C. Ye,  Come-closer-diffuse-faster: Accelerating conditional diffusion models for inverse problems through stochastic contraction,  in: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, 2022c, pp. 12413–12422.
  * Hyvärinen and Dayan [2005] A. Hyvärinen, P. Dayan,  Estimation of non-normalized statistical models by score matching,  Journal of Machine Learning Research 6 (2005).
  * Vincent [2011] P. Vincent,  A connection between score matching and denoising autoencoders,  Neural computation 23 (2011) 1661–1674.
  * Song and Ermon [2019] Y. Song, S. Ermon,  Generative modeling by estimating gradients of the data distribution,  Advances in Neural Information Processing Systems 32 (2019).
  * Kobyzev et al. [2020] I. Kobyzev, S. J. Prince, M. A. Brubaker,  Normalizing flows: An introduction and review of current methods,  IEEE Transactions on Pattern Analysis and Machine Intelligence 43 (2020) 3964–3979.
  * Saharia et al. [2022] C. Saharia, J. Ho, W. Chan, T. Salimans, D. J. Fleet, M. Norouzi,  Image super-resolution via iterative refinement,  IEEE Transactions on Pattern Analysis and Machine Intelligence 45 (2022) 4713–4726.
  * Song et al. [2020] Y. Song, J. Sohl-Dickstein, D. P. Kingma, A. Kumar, S. Ermon, B. Poole,  Score-based generative modeling through stochastic differential equations,  arXiv preprint arXiv:2011.13456 (2020).
  * Batzolis et al. [2021] G. Batzolis, J. Stanczuk, C.-B. Schönlieb, C. Etmann,  Conditional image generation with score-based diffusion models,  arXiv preprint arXiv:2111.13606 (2021).
  * Song and Ermon [2020] Y. Song, S. Ermon,  Improved techniques for training score-based generative models,  Advances in Neural Information Processing Systems 33 (2020) 12438–12448.
  * Jacobsen et al. [2023] C. Jacobsen, Y. Zhuang, K. Duraisamy,  CoCoGen: Physically-consistent and conditioned score-based generative models for forward and inverse problems,  arXiv preprint arXiv:2312.10527 (2023).
  * Bastek et al. [2024] J.-H. Bastek, W. C. Sun, D. M. Kochmann,  Physics-informed diffusion models,  arXiv preprint arXiv:2403.14404 (2024).
  * Ongie et al. [2020] G. Ongie, A. Jalal, C. A. Metzler, R. G. Baraniuk, A. G. Dimakis, R. Willett,  Deep learning techniques for inverse problems in imaging,  IEEE Journal on Selected Areas in Information Theory 1 (2020) 39–56.
  * Patel et al. [2022] D. V. Patel, D. Ray, A. A. Oberai,  Solution of physics-based bayesian inverse problems with deep generative priors,  Computer Methods in Applied Mechanics and Engineering 400 (2022) 115428\.
  * Bohra et al. [2022] P. Bohra, T.-a. Pham, J. Dong, M. Unser,  Bayesian inversion for nonlinear imaging models using deep generative priors,  arXiv preprint arXiv:2203.10078 (2022).
  * Whang et al. [2021] J. Whang, E. Lindgren, A. Dimakis,  Composing normalizing flows for inverse problems,  in: International Conference on Machine Learning, PMLR, 2021, pp. 11158–11169.
  * Adler and Öktem [2018] J. Adler, O. Öktem,  Deep Bayesian inversion,  arXiv preprint arXiv:1811.05910 (2018).
  * Ray et al. [2022] D. Ray, H. Ramaswamy, D. V. Patel, A. A. Oberai,  The efficacy and generalizability of conditional GANs for posterior inference in physics-based inverse problems,  Numerical Algebra, Control and Optimization (2022).
  * Ray et al. [2023] D. Ray, J. Murgoitio-Esandi, A. Dasgupta, A. A. Oberai,  Solution of physics-based inverse problems using conditional generative adversarial networks with full gradient penalty,  Computer Methods in Applied Mechanics and Engineering 417 (2023) 116338\.
  * Padmanabha and Zabaras [2021] G. A. Padmanabha, N. Zabaras,  Solving inverse problems using conditional invertible neural networks,  Journal of Computational Physics 433 (2021) 110194.
  * Sun and Bouman [2021] H. Sun, K. L. Bouman,  Deep probabilistic imaging: Uncertainty quantification and multi-modal solution characterization for computational imaging,  in: Proceedings of the AAAI Conference on Artificial Intelligence, volume 35, 2021, pp. 2628–2637.
  * Dasgupta and Di [2021] A. Dasgupta, Z. W. Di,  Uncertainty quantification for ptychography using normalizing flows,  in: Fourth Workshop on Machine Learning and the Physical Sciences, 2021.
  * Feng et al. [2023] Y. Feng, K. Tang, X. Wan, Q. Liao,  Dimension-reduced KRnet maps for high-dimensional inverse problems,  arXiv preprint arXiv:2303.00573 (2023).
  * Dasgupta et al. [2024] A. Dasgupta, D. V. Patel, D. Ray, E. A. Johnson, A. A. Oberai,  A dimension-reduced variational approach for solving physics-based inverse problems using generative adversarial network priors and normalizing flows,  Computer Methods in Applied Mechanics and Engineering 420 (2024) 116682\.
  * Graikos et al. [2022] A. Graikos, N. Malkin, N. Jojic, D. Samaras,  Diffusion models as plug-and-play priors,  Advances in Neural Information Processing Systems 35 (2022) 14715–14728.
  * Mardani et al. [2023] M. Mardani, J. Song, J. Kautz, A. Vahdat,  A variational perspective on solving inverse problems with diffusion models,  arXiv preprint arXiv:2305.04391 (2023).
  * Metropolis et al. [1953] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, E. Teller,  Equation of state calculations by fast computing machines,  The Journal of Chemical Physics 21 (1953) 1087–1092.
  * Hastings [1970] W. K. Hastings,  Monte Carlo sampling methods using Markov chains and their applications (1970).
  * Lin et al. [2017] G. Lin, A. Milan, C. Shen, I. Reid,  Refinenet: Multi-path refinement networks for high-resolution semantic segmentation,  in: Proceedings of the IEEE conference on computer vision and pattern recognition, 2017, pp. 1925–1934.
  * Pavan et al. [2012] T. Z. Pavan, E. L. Madsen, G. R. Frank, J. Jiang, A. A. Carneiro, T. J. Hall,  A nonlinear elasticity phantom containing spherical inclusions,  Physics in medicine & biology 57 (2012) 4787.
  * Lu et al. [2021] G. Lu, R. Li, X. Qian, R. Chen, L. Jiang, Z. Chen, K. K. Shung, M. S. Humayun, Q. Zhou,  Layer-specific ultrasound elastography using a multi-layered shear wave dispersion model for assessing the viscoelastic properties,  Physics in Medicine & Biology 66 (2021) 035003.
  * Foo et al. [2024] K. Y. Foo, B. Shaddy, J. Murgoitio-Esandi, M. S. Hepburn, J. Li, A. Mowla, D. Vahala, S. E. Amos, Y. S. Choi, A. A. Oberai, et al.,  Tumor spheroid elasticity estimation using mechano-microscopy combined with a conditional generative adversarial network,  bioRxiv (2024) 2024–04.
  * Paszke et al. [2019] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, A. Desmaison, A. Kopf, E. Yang, Z. DeVito, M. Raison, A. Tejani, S. Chilamkurthy, B. Steiner, L. Fang, J. Bai, S. Chintala,  PyTorch: An imperative style, high-performance deep learning library,  in: Advances in Neural Information Processing Systems, 2019, pp. 8024–8035.
  * Barbone and Oberai [2010] P. E. Barbone, A. A. Oberai,  A review of the mathematical and computational foundations of biomechanical imaging,  Computational Modeling in Biomechanics (2010) 375–408.
  * Logg et al. [2012] A. Logg, K.-A. Mardal, G. N. Wells, et al., Automated Solution of Differential Equations by the Finite Element Method, Springer, 2012. doi:10.1007/978-3-642-23099-8.
  * Jonas et al. [1991] J. Jonas, C. Y. Mardin, U. Schlötzer-Schrehardt, G. Naumann,  Morphometry of the human lamina cribrosa surface.,  Investigative ophthalmology & visual science 32 (1991) 401–405.
  * Park et al. [2015] S. C. Park, J. Brumm, R. L. Furlanetto, C. Netto, Y. Liu, C. Tello, J. M. Liebmann, R. Ritch,  Lamina cribrosa depth in different stages of glaucoma,  Investigative ophthalmology & visual science 56 (2015) 2059–2064.
  * Vurgese et al. [2012] S. Vurgese, S. Panda-Jonas, J. B. Jonas,  Scleral thickness in human eyes,  PloS one 7 (2012) e29692.
  * Alamouti and Funk [2003] B. Alamouti, J. Funk,  Retinal thickness decreases with age: an oct study,  British journal of ophthalmology 87 (2003) 899–901.
  * Yic et al. [2023] C. D. Yic, J. Pontet, M. Mercado, M. Muñoz, A. Biestro,  Ultrasonographic measurement of the optic nerve sheath diameter to detect intracranial hypertension: an observational study,  The Ultrasound Journal 15 (2023) 4.
  * Hua et al. [2017] Y. Hua, J. Tong, D. Ghate, S. Kedar, L. Gu,  Intracranial pressure influences the behavior of the optic nerve head,  Journal of biomechanical engineering 139 (2017) 031003.
  * Şahan et al. [2019] M. H. Şahan, A. Doğan, M. İnal, M. Alpua, N. Asal,  Evaluation of the optic nerve by strain and shear wave elastography in patients with migraine,  Journal of Ultrasound in Medicine 38 (2019) 1153–1161.
  * Qian et al. [2021] X. Qian, R. Li, G. Lu, L. Jiang, H. Kang, K. K. Shung, M. S. Humayun, Q. Zhou,  Ultrasonic elastography to assess biomechanical properties of the optic nerve head and peripapillary sclera of the eye,  Ultrasonics 110 (2021) 106263.
  * Zhang et al. [2020] L. Zhang, M. R. Beotra, M. Baskaran, T. A. Tun, X. Wang, S. A. Perera, N. G. Strouthidis, T. Aung, C. Boote, M. J. Girard,  In vivo measurements of prelamina and lamina cribrosa biomechanical properties in humans,  Investigative Ophthalmology & Visual Science 61 (2020) 27–27.
  * Dassault Systèmes [2022] Dassault Systèmes, Simulia 2022 User Assistance, 2022.
  * Kennedy et al. [2014a] B. F. Kennedy, K. M. Kennedy, D. D. Sampson,  A review of optical coherence elastography: Fundamentals, techniques and prospects,  IEEE Journal of Selected Topics in Quantum Electronics 20 (2014a) 272–288. doi:10.1109/JSTQE.2013.2291445.
  * Kennedy et al. [2014b] B. F. Kennedy, R. A. McLaughlin, K. M. Kennedy, L. Chin, A. Curatolo, A. Tien, B. Latham, C. M. Saunders, D. D. Sampson,  Optical coherence micro-elastography: mechanical-contrast imaging of tissue microstructure,  Biomedical optics express 5 (2014b) 2113–2124.
  * Mowla et al. [2022] A. Mowla, J. Li, M. S. Hepburn, S. Maher, L. Chin, G. C. Yeoh, Y. S. Choi, B. F. Kennedy,  Subcellular mechano-microscopy: High resolution three-dimensional elasticity mapping using optical coherence microscopy,  Optics Letters 47 (2022) 3303–3306. doi:10.1364/OL.451681.
  * Mowla et al. [2024] A. Mowla, M. S. Hepburn, J. Li, D. Vahala, S. E. Amos, L. M. Hirvonen, R. W. Sanderson, P. Wijesinghe, S. Maher, Y. S. Choi, B. F. Kennedy,  Multimodal mechano-microscopy reveals mechanical phenotypes of breast cancer spheroids in three dimensions,  bioRxiv (2024) 2024.04.05.588260. doi:10.1101/2024.04.05.588260.
  * Olofsson et al. [2018] K. Olofsson, V. Carannante, M. Ohlin, T. Frisk, K. Kushiro, M. Takai, A. Lundqvist, B. Önfelt, M. Wiklund,  Acoustic formation of multicellular tumor spheroids enabling on-chip functional and structural imaging,  Lab on a Chip 18 (2018) 2466–2476.
  * Imai et al. [1985] H. Imai, M. Iri, K. Murota,  Voronoi diagram in the laguerre geometry and its applications,  SIAM Journal on Computing 14 (1985) 93–105. doi:10.1137/0214006.
  * Barbone and Gokhale [2004] P. E. Barbone, N. H. Gokhale,  Elastic modulus imaging: on the uniqueness and nonuniqueness of the elastography inverse problem in two dimensions,  Inverse problems 20 (2004) 283.
  * Ferreira et al. [2012] E. R. Ferreira, A. A. Oberai, P. E. Barbone,  Uniqueness of the elastography inverse problem for incompressible nonlinear planar hyperelasticity,  Inverse problems 28 (2012) 065008.
  * He et al. [2022] K. He, X. Chen, S. Xie, Y. Li, P. Dollár, R. Girshick,  Masked autoencoders are scalable vision learners,  in: Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, 2022, pp. 16000–16009.

##  Appendix Appendix A Derivation of Eq. 10 from Eq. 8

Using Eq. 9,

$\mathrm{p}_{\tilde{{\mathbf{X}}}|{\mathbf{Y}}}(\tilde{{\mathbf{x}}}|{\mathbf{y}})=\int\mathrm{p}_{\tilde{{\mathbf{X}}}|{\mathbf{X}}}\\!\left({\tilde{{\mathbf{x}}}|{\mathbf{x}}}\right)\;\mathrm{p}_{{\mathbf{X}}|{\mathbf{Y}}}({\mathbf{x}}|{\mathbf{y}})\;\mathrm{d}{\mathbf{x}},$
(29)

we wish to show the denoising score matching objective Eq. 8

$\ell_{\mathrm{DSM}}({\mathbf{y}};\bm{\theta})=\frac{1}{2}\mathbb{E}_{\tilde{{\mathbf{X}}}\sim\mathrm{p}_{\tilde{{\mathbf{X}}}|{\mathbf{Y}}}}\Big{[}\lVert\mathbf{s}(\tilde{{\mathbf{x}}},{\mathbf{y}};\bm{\theta})-\nabla_{\tilde{{\mathbf{x}}}}\log\mathrm{p}_{\tilde{{\mathbf{X}}}|{\mathbf{Y}}}\\!(\tilde{{\mathbf{x}}}|{\mathbf{y}})\rVert_{2}^{2}\Big{]},$
(30)

is equivalent to Eq. 10, i.e.,

$\ell_{\mathrm{DSM}}({\mathbf{y}};\bm{\theta})=\frac{1}{2}\mathbb{E}_{{\mathbf{X}}\sim\mathrm{p}_{{\mathbf{X}}|{\mathbf{Y}}}}\Bigg{[}\mathbb{E}_{\tilde{{\mathbf{X}}}\sim\mathrm{p}_{\tilde{{\mathbf{X}}}|{\mathbf{X}}}}\Big{[}\lVert\mathbf{s}(\tilde{{\mathbf{x}}},{\mathbf{y}};\bm{\theta})-\nabla_{\tilde{{\mathbf{x}}}}\log\mathrm{p}_{\tilde{{\mathbf{X}}}|{\mathbf{X}}}\\!\left({\tilde{{\mathbf{x}}}|{\mathbf{x}}}\right)\rVert_{2}^{2}\Big{]}\Bigg{]}+\mathcal{C},$
(31)

where $\mathcal{C}$ is an additive constant that does not depend on
$\bm{\theta}$. To show this, we start with Eq. 30 and expand the integrand

$\begin{split}\ell_{\mathrm{DSM}}({\mathbf{y}};\bm{\theta})&=\frac{1}{2}\mathbb{E}_{\tilde{{\mathbf{X}}}\sim\mathrm{p}_{\tilde{{\mathbf{X}}}|{\mathbf{Y}}}}\Big{[}\lVert\mathbf{s}(\tilde{{\mathbf{x}}},{\mathbf{y}};\bm{\theta})-\nabla_{\tilde{{\mathbf{x}}}}\log\mathrm{p}_{\tilde{{\mathbf{X}}}|{\mathbf{Y}}}\\!(\tilde{{\mathbf{x}}}|{\mathbf{y}})\rVert_{2}^{2}\Big{]}\\\
&=\frac{1}{2}\mathbb{E}_{\tilde{{\mathbf{X}}}\sim\mathrm{p}_{\tilde{{\mathbf{X}}}|{\mathbf{Y}}}}\Big{[}\Big{<}\mathbf{s}(\tilde{{\mathbf{x}}},{\mathbf{y}};\bm{\theta}),\mathbf{s}(\tilde{{\mathbf{x}}},{\mathbf{y}};\bm{\theta})\Big{>}-2\Big{<}\mathbf{s}(\tilde{{\mathbf{x}}},{\mathbf{y}};\bm{\theta}),\nabla_{\tilde{{\mathbf{x}}}}\log\mathrm{p}_{\tilde{{\mathbf{X}}}|{\mathbf{Y}}}\\!(\tilde{{\mathbf{x}}}|{\mathbf{y}})\Big{>}\\\
&\hphantom{=\frac{1}{2}\mathbb{E}_{\tilde{{\mathbf{X}}}\sim\mathrm{p}_{\tilde{{\mathbf{X}}}|{\mathbf{Y}}}}\Big{[}\Big{<}\mathbf{s}(\tilde{{\mathbf{x}}},{\mathbf{y}};\bm{\theta}),\mathbf{s}(\tilde{{\mathbf{x}}},{\mathbf{y}};\bm{\theta})\Big{>}}+\Big{<}\nabla_{\tilde{{\mathbf{x}}}}\log\mathrm{p}_{\tilde{{\mathbf{X}}}|{\mathbf{Y}}}\\!(\tilde{{\mathbf{x}}}|{\mathbf{y}}),\nabla_{\tilde{{\mathbf{x}}}}\log\mathrm{p}_{\tilde{{\mathbf{X}}}|{\mathbf{Y}}}\\!(\tilde{{\mathbf{x}}}|{\mathbf{y}})\Big{>}\Big{]},\end{split}$
(32)

where $<\cdot,\cdot>$ denotes the Euclidean inner product in
$\mathbb{R}^{n_{\mathcal{X}}}$, i.e.,
$<{\mathbf{a}},{\mathbf{a}}>\,={\mathbf{a}}^{\mathrm{T}}{\mathbf{a}}$ . Note
that the last term Eq. 32 does not depend on $\bm{\theta}$, so, herein, we
will denote the last term using $\mathcal{C}_{1}^{\prime}$. Now, we expand the
expectations in Eq. 32 and utilize Eq. 29 to rewrite Eq. 32 as follows:

$\begin{split}\ell_{\mathrm{DSM}}({\mathbf{y}};\bm{\theta})&=\frac{1}{2}\Bigg{\\{}\left[\int\Big{<}\mathbf{s}(\tilde{{\mathbf{x}}},{\mathbf{y}};\bm{\theta}),\mathbf{s}(\tilde{{\mathbf{x}}},{\mathbf{y}};\bm{\theta})\Big{>}\;\mathrm{p}_{\tilde{{\mathbf{X}}}|{\mathbf{Y}}}\\!\left({\tilde{{\mathbf{x}}}|{\mathbf{y}}}\right)\mathrm{d}\tilde{{\mathbf{x}}}\right]\\\
&\hphantom{=\mathbf{s}(\tilde{{\mathbf{x}}},{\mathbf{y}};\bm{\theta})}-2\left[\int\Big{<}\mathbf{s}(\tilde{{\mathbf{x}}},{\mathbf{y}};\bm{\theta}),\nabla_{\tilde{{\mathbf{x}}}}\log\mathrm{p}_{\tilde{{\mathbf{X}}}|{\mathbf{Y}}}\\!(\tilde{{\mathbf{x}}}|{\mathbf{y}})\Big{>}\;\mathrm{p}_{\tilde{{\mathbf{X}}}|{\mathbf{Y}}}\\!\left({\tilde{{\mathbf{x}}}|{\mathbf{y}}}\right)\mathrm{d}\tilde{{\mathbf{x}}}\right]\Bigg{\\}}+\mathcal{C}_{1},\end{split}$
(33)

where
$\mathcal{C}_{1}=\frac{1}{2}\mathbb{E}_{\tilde{{\mathbf{X}}}\sim\mathrm{p}_{\tilde{{\mathbf{X}}}|{\mathbf{Y}}}}\\{\mathcal{C}_{1}^{\prime}\\}$.
Utilizing

$\nabla_{\tilde{{\mathbf{x}}}}\log\mathrm{p}_{\tilde{{\mathbf{X}}}|{\mathbf{Y}}}\\!(\tilde{{\mathbf{x}}}|{\mathbf{y}})=\frac{1}{\mathrm{p}_{\tilde{{\mathbf{X}}}|{\mathbf{Y}}}\\!(\tilde{{\mathbf{x}}}|{\mathbf{y}})}\nabla_{\tilde{{\mathbf{x}}}}\mathrm{p}_{\tilde{{\mathbf{X}}}|{\mathbf{Y}}}\\!(\tilde{{\mathbf{x}}}|{\mathbf{y}}),$
(34)
$\nabla_{\tilde{{\mathbf{x}}}}\log\mathrm{p}_{\tilde{{\mathbf{X}}}|{\mathbf{X}}}\\!(\tilde{{\mathbf{x}}}|{\mathbf{x}})=\frac{1}{\mathrm{p}_{\tilde{{\mathbf{X}}}|{\mathbf{X}}}\\!(\tilde{{\mathbf{x}}}|{\mathbf{x}})}\nabla_{\tilde{{\mathbf{x}}}}\mathrm{p}_{\tilde{{\mathbf{X}}}|{\mathbf{X}}}\\!(\tilde{{\mathbf{x}}}|{\mathbf{x}}),$
(35)

and Eq. 29, we can simplify the second term in Eq. 33 as follows:

$\begin{split}&\int\Big{<}\mathbf{s}(\tilde{{\mathbf{x}}},{\mathbf{y}};\bm{\theta}),\nabla_{\tilde{{\mathbf{x}}}}\log\mathrm{p}_{\tilde{{\mathbf{X}}}|{\mathbf{Y}}}\\!(\tilde{{\mathbf{x}}}|{\mathbf{y}})\Big{>}\mathrm{p}_{\tilde{{\mathbf{X}}}|{\mathbf{Y}}}\\!\left({\tilde{{\mathbf{x}}}|{\mathbf{y}}}\right)\mathrm{d}\tilde{{\mathbf{x}}}=\int\Big{<}\mathbf{s}(\tilde{{\mathbf{x}}},{\mathbf{y}};\bm{\theta}),\nabla_{\tilde{{\mathbf{x}}}}\mathrm{p}_{\tilde{{\mathbf{X}}}|{\mathbf{Y}}}\\!(\tilde{{\mathbf{x}}}|{\mathbf{y}})\Big{>}\;\mathrm{d}\tilde{{\mathbf{x}}}\\\
&\hphantom{\int\Big{<}\mathbf{s}(\tilde{{\mathbf{x}}},{\mathbf{y}};\bm{\theta}),\nabla_{\tilde{{\mathbf{x}}}}\log\mathrm{p}_{\tilde{{\mathbf{X}}}|{\mathbf{Y}}}\\!(\tilde{{\mathbf{x}}}|{\mathbf{y}})\Big{>}}=\int\Big{<}\mathbf{s}(\tilde{{\mathbf{x}}},{\mathbf{y}};\bm{\theta}),\nabla_{\tilde{{\mathbf{x}}}}\underbrace{\int\mathrm{p}_{\tilde{{\mathbf{X}}}|{\mathbf{X}}}\\!\left({\tilde{{\mathbf{x}}}|{\mathbf{x}}}\right)\;\mathrm{p}_{{\mathbf{X}}|{\mathbf{Y}}}({\mathbf{x}}|{\mathbf{y}})\;\mathrm{d}{\mathbf{x}}}_{\mathrm{p}_{\tilde{{\mathbf{X}}}|{\mathbf{Y}}}\\!(\tilde{{\mathbf{x}}}|{\mathbf{y}})}\Big{>}\;\mathrm{d}\tilde{{\mathbf{x}}}\\\