# Bayesian Model Calibration for Block Copolymer Self-Assembly: Likelihood-
Free Inference and Expected Information Gain Computation via Measure Transport

Ricardo Baptista11footnotemark: 1<EMAIL_ADDRESS>Lianghao Cao22footnotemark: 2
<EMAIL_ADDRESS>Joshua Chen33footnotemark: 3<EMAIL_ADDRESS>Omar Ghattas<EMAIL_ADDRESS>Fengyi Li44footnotemark: 4<EMAIL_ADDRESS>Youssef M. Marzouk<EMAIL_ADDRESS>J. Tinsley Oden<EMAIL_ADDRESS>Massachusetts Institute of Technology, Cambridge, MA 02139, USA  The
University of Texas at Austin, Austin, TX 78712, USA

###### Abstract

We consider the Bayesian calibration of models describing the phenomenon of
block copolymer (BCP) self-assembly using image data produced by microscopy or
X-ray scattering techniques. To account for the long-range disorder in BCP
equilibrium structures we introduce auxiliary variables to represent this
aleatory uncertainty. These variables, however, result in an integrated
likelihood for high-dimensional image data that is generally intractable to
evaluate. We tackle this challenging Bayesian inference problem using a
likelihood-free approach based on measure transport together with the
construction of summary statistics for the image data. We also show that
expected information gains (EIG) from the observed data about the model
parameters can be computed with no significant additional cost. Lastly, we
present a numerical case study based on the Ohta–Kawasaki model for diblock
copolymer thin film self-assembly and top-down microscopy characterization.
For calibration, we introduce several domain-specific energy- and Fourier-
based summary statistics and quantify their informativeness using EIG. We
demonstrate the power of the proposed approach to study the effect of various
data corruptions and experimental designs on the calibration results.

###### keywords:

block copolymers , material self-assembly , uncertainty quantification ,
likelihood-free inference , summary statistics , measure transport , expected
information gain , the Ohta–Kawasaki model

fn1fn1footnotetext: Authors contributed equally to this manuscript.

###### Contents

  1. 1 Introduction
    1. 1.1 The self-assembly of block copolymers: Application and modeling
    2. 1.2 Calibration of block copolymer self-assembly models: Uncertainties and the Bayesian approach
    3. 1.3 Likelihood-free inference and expected information gain computation via measure transport
    4. 1.4 Related work
    5. 1.5 Paper layout
    6. 1.6 Contributions
  2. 2 A probabilistic data model of block copolymer self-assembly
    1. 2.1 Statistical models of experimental characterization
    2. 2.2 The aleatoric uncertainty and statistical models of block copolymer self-assembly
    3. 2.3 A generic data model of block copolymer self-assembly
  3. 3 Bayesian inference and quantifying uncertainty
    1. 3.1 The posterior and joint posterior
    2. 3.2 Entropy and information gain
    3. 3.3 Expected information gain
  4. 4 Addressing the challenges of auxiliary variables and high-dimensional image data via likelihood-free inference
    1. 4.1 Evaluating integrated likelihoods
    2. 4.2 Inference methods for intractable likelihoods
    3. 4.3 Likelihood-free inference: Related work and practical desiderata
    4. 4.4 Informative summary statistics of high-dimensional image data
  5. 5 Likelihood-free inference via triangular transport maps
    1. 5.1 Measure transport
    2. 5.2 The Knothe–Rosenblatt rearrangement
    3. 5.3 Estimating the Knothe–Rosenblatt rearrangement
    4. 5.4 Conditional density estimation via triangular transport maps
  6. 6 Estimating expected information gain via triangular transport maps
    1. 6.1 Related work
    2. 6.2 Likelihood-free estimators for expectation information gain
  7. 7 Case study: Model definitions
    1. 7.1 The Ohta–Kawasaki model
    2. 7.2 The forward operator of the Ohta–Kawasaki model
    3. 7.3 A model of top-down microscopy characterization
    4. 7.4 Summary statistics for top-down microscope images of diblock copolymer thin films
      1. 7.4.1 Energy-based summary statistics
      2. 7.4.2 Fourier-based summary statistics
    5. 7.5 Prior distribution
  8. 8 Case study: Numerical results
    1. 8.1 Posterior and posterior predictives
    2. 8.2 Expected information gain: Energy-based summary statistics
    3. 8.3 Expected information gain: Fourier-based summary statistics
    4. 8.4 Expected information gain: Levels of magnification
    5. 8.5 Expected information gain: Known volume ratio of monomer blocks
  9. 9 Conclusion
  10. Acknowledgement
  11. A The auxiliary variable: a transformed Gaussian random field
  12. B The discretization of the state space, the restriction operators, and image interpolation
  13. C Parameterization and computation of transport maps

##  1 Introduction

###  1.1 The self-assembly of block copolymers: Application and modeling

Block copolymer (BCP) melts are large collections of one or more types of
BCPs, or polymers consisting of several blocks of distinct monomers joined
together by covalent bonds. Originally mixed in high temperature equilibrium,
the distinct monomer blocks in a BCP melt spontaneously segregate upon thermal
or solvent annealing as a result of their thermodynamic incompatibility.
Periodically ordered nanoscale structures emerge as the melt achieves a new
equilibrium configuration. This physical phenomenon is known as the self-
assembly of BCPs [1, 2]. Examples of the equilibrium structures for the self-
assembly of diblock copolymers (Di-BCPs), i.e., BCPs with two distinct blocks,
are lamellae, spheres, cylinders, etc.; see Fig. 1 (_left_) along with an
experimental phase diagram. The massive attention given to BCP self-assembly
stems from the flexibility in their achievable equilibrium structures when
subject to external guidance. Specifically, carefully engineered techniques,
such as substrate-based techniques [3, 4, 5] or techniques using thermal and
electric fields [6, 7, 8], can influence the self-assembly process to form
equilibrium structures with desired features. The application of these
techniques is known as _directed self-assembly_ (DSA) of BCPs. The DSA of BCPs
has quickly attracted attention for its wide-ranging commercial applications,
the most notable being _patterning_ for lithography at the sub-10 nm scale [9,
10, 11, 12].

|
---|---
Figure 1: (_left_) The experimental phase diagram for the self-assembly of
Polyisoprene-Polystyrene diblock copolymers. The parameter $f_{PI}$ denotes
the volume ratio of the Polyisoprene block, $\chi$ denotes the Flory–Huggins
interaction parameter, and $N$ denotes the degree of polymerization. Reprinted
(adapted) with permission from [13]. Copyright ©1995 American Chemical
Society. (_right_) Top-down scanning electron microscope images of lamella
phase PS-b-PMMA copolymer films with (A) a chemically neutral surface; (B) a
chemically nanopatterned surface with period smaller than the natural period
of the copolymer, thus leading to defects; (C) a chemically nanopatterned
surface with period close to the natural period of the copolymer. Reprinted
(adapted) with permission from [4]. Copyright ©2003 Springer Nature.

Computer simulation of BCP self-assembly based on mathematical models has
played an important role in the study of self-assembly. These simulations have
proven to be valuable in determining the effects of material properties, film
thickness, and polymer-substrate interactions on the self-assembly [14]. The
popular models and simulation techniques for self-assembly are predominantly
derived from statistical polymer theory, such as Monte Carlo simulations [15,
16], theoretically-informed coarse-grained (TICG) simulation [17], self-
consistent field theory (SCFT) [18, 19], and density functional theory (DFT)
[20, 21]—listed in decreasing order by granularity55footnotemark: 5. While
model-based studies of the self-assembly mainly serve as aides to interpreting
empirical findings, there is an increased interest in using models to make
potentially consequential predictions. Most noticeably, recent efforts have
been made on the _inverse design_ 666We define the granularity of a model as
the finest length scale at which it can faithfully resolve polymer
conformation, monomer interaction, and thermal fluctuation.of DSA, where model
simulations and iterative optimization schemes suggest designs of external
guidance and BCP material properties that are best suited for achieving a
target equilibrium pattern [22, 23, 24, 25, 26]. Model-based inverse design is
often considered the critical pathway towards transitioning DSA technology
from the laboratory into commercial manufacturing of integrated circuit
devices [27].

###  1.2 Calibration of block copolymer self-assembly models: Uncertainties
and the Bayesian approach

With the increasing usage and growing impact of model-based predictions for
BCP self-assembly, it is increasingly important to ensure the reliability of
models. Model parameter estimation, or model calibration, is necessary in
order to produce reliable model-based predictions [28] and is the focus of
this work. In particular, we formulate and propose a suitable methodology for
calibrating models of BCP self-assembly. Models are calibrated according to
observational data obtained by experimental characterizations [29, Vol. 2],
via microscopy or X-ray scattering techniques (see, for example, Fig. 2), of
BCP melts, and various intrinsic uncertainties are also explicitly accounted
for.

One source of intrinsic uncertainty stems from corruptions introduced in the
characterization process. This is often modeled as stochastic sensor noise. We
also address an _aleatoric uncertainty_ , related to the well-known phenomenon
that unguided or weakly guided self-assembly of BCPs results in a range of
possible equilibrium morphologies with long-range disorder. A characteristic
example for such long-range disorder is the ‘finger-print’ patterns for the
lamellae phase of Di-BCPs thin films under unguided self-assembly. For DSA,
the aleatoric uncertainty results in _defects_ in the equilibrium patterns.
Examples of a ‘finger-print’ pattern and DSA pattern defects are shown in Fig.
1 (_right_). Pattern defects have significant implications on product yield
and quality in the nanolithography application of BCP self-assembly [11, 30].

Statistical polymer theory [31] accounts for the impact of aleatoric
uncertainty by describing the variety of achievable equilibrium morphologies
for a particular polymer melt system via a _statistical ensemble_. The
statistical ensemble assigns a probability of occurrence for each achievable
equilibrium morphology [32]. In this work, we consider all BCP self-assembly
models as _statistical models_ , through which the computer-generated
solutions are samples from a model-specific statistical ensemble. We assume
the stochastic component of the models is accounted for by _auxiliary random
variables_ which are possibly infinite-dimensional and/or parameter-dependent.

To perform model calibration that accounts for these uncertainties, we
advocate for a Bayesian probabilistic approach that uses probability
distributions to describe the uncertainty in model parameters. Under such an
approach, the output of the model calibration is a posterior probability
distribution of the model parameters conditioned on the observed data. In
addition to performing parameter estimation, we are also interested in
studying the effects of instrument-specific parameters that describe the
strength of various kinds of corruptions introduced in the characterization
process. These effects can be quantified by the amount of information content,
i.e., the Shannon entropy [33], preserved by the characterization process.
Doing so informs decision-making for appropriate specifications of
characterization instruments, so that any experimental data are as effective
as possible for model calibration. For the numerical case study presented in
Sections 7 and 8, we focus on top-down microscopy characterization of Di-BCP
thin films, and study the effect of (i) various degrees of image blurring
and/or noise intensities, and (ii) various levels of magnification on the
calibration of the Ohta–Kawasaki (OK) model [34, 21, 35], a DFT model, for Di-
BCP self-assembly.

###  1.3 Likelihood-free inference and expected information gain computation
via measure transport

To solve the proposed Bayesian model calibration problem, we must first tackle
key challenges that may seem insurmountable, especially if addressed via
conventional likelihood-based methods for Bayesian inference. Specifically,
analytical evaluations of integrated likelihood functions [36] are deemed
intractable due to (i) the high-dimensional auxiliary random variables, and
(ii) the potentially complex forms of sensor noise present in the
characterization data. To address these challenges, we adopt the measure
transport approach to likelihood-free inference (LFI) that was considered in
[37, 38]. This approach addresses the challenges of intractable likelihoods,
given that it only requires samples from the joint distribution of the
parameters and model-generated data. Moreover, the generation of these samples
can be performed offline in a scalable manner.

Performing LFI with high-dimensional image data can be challenging, however,
as the regions of the parameter space with large likelihood values are
typically highly concentrated. Therefore, we advocate constructing low-
dimensional summary statistics of the images to capture information in the
data that is relevant to the model parameters of interest. This approach is
commonly used for the analysis of the image data for BCP self-assembly [39,
40, 41]. Furthermore, using summary statistics is natural for the model
calibration of strongly guided Di-BCP thin film self-assembly using X-ray
scattering image data [24, 42], as these scattering images are inherently low-
dimensional. In our numerical case study in Section 7, we propose several
energy-based and Fourier-based summary statistics of top-down microscope
images of Di-BCP thin film self-assembly, and their effectiveness for model
calibration is quantitatively studied.

The ability to produce posterior densities for Bayesian calibration of BCP
models is a prerequisite for the more computationally difficult task of
computing the expected information gain (EIG) from summary statistics of image
data under various forms of corruptions introduced by the characterization
process. To do so, we present a novel likelihood-free estimator for EIG based
on measure transport estimators for the posterior distribution. While several
EIG estimators based on variational approximations have been recently proposed
[43, 44], the estimators we use in this work possess nice theoretical
properties, such as having the potential of being consistent in the infinite
sample setting. Furthermore, we show that the estimation of the EIG can be
computed using the proposed LFI approach with no significant additional cost.

###  1.4 Related work

Parameter estimation of models for self-assembly based on experimental
characterization of the equilibrium patterns has been previously studied.
Khaira et. al [45] determined material and geometric parameters in the TICG
model for chemical substrate guided self-assembly of lamellae phase Di-BCP
thin films. The parameter estimation is based on data obtained from resonant
critical-dimension small-angle X-ray scattering (res-CDSAXS) and grazing-
incidence small-angle X-ray scattering (GISAXS) characterization data. Hannon
et. al [42] accomplished similar calibration tasks for the SCFT model with
res-CDSAXS data.

The parameter estimation in previous works generally follow a _deterministic_
approach. Since deterministic approaches do not quantify the impact of various
intrinsic uncertainties on calibration confidence, the reliability of
calibrated models may be poor. In particular, by calibrating only to
characterization data of Di-BCP equilibrium states strongly influenced by DSA
external guidance, the aleatoric uncertainty can be reasonably omitted, since
it has a negligible effect on the data. However, in the scenario where the
self-assembly is unguided or weakly guided, the aleatoric uncertainty cannot
be neglected, as model predictions and experimental characterizations can
differ in nontrivial ways due to the long-range disorder or abundance of
defected states.

We present related works in LFI and estimation of EIG in Sections 4.3 and 6.1,
respectively.

###  1.5 Paper layout

The paper is organized as follows. In Section 2, the relevant components of
probabilistic data models for BCP self-assembly are defined. These components
include the forward operator of computer simulations and the statistical model
for experimental characterization. In Section 3, we introduce the Bayesian
framework for the calibration of BCP models and provide a short review of
relevant concepts from information theory. In Section 4, we outline the
challenges in the Bayesian calibration of BCP models and discuss how these
challenges can be addressed with the LFI approach and constructing summary
statistics of image data. In Section 5, we present a measure transport
approach to LFI for estimating the posterior distribution. In Section 6, we
introduce a methodology for estimating EIG via triangular transport maps. In
Section 7, we describe the data models for the numerical case study: Bayesian
calibration of the OK model for the self-assembly of Di-BCP thin films with
top-down microscopy image data. Several energy- and Fourier-based summary
statistics of the images are proposed. The prior distribution is defined via a
careful analysis of an admissible parameter regime for the OK model. Lastly,
in Section 8 we present numerical results for the case study obtained through
the proposed LFI and EIG estimation methods.

###  1.6 Contributions

The contributions of this work are summarized, in order of appearance, as
follows.

  1. 1.

We introduce a Bayesian approach for calibration of a BCP self-assembly model
that accounts for both measurement and aleatoric uncertainties.

  2. 2.

We overview a quantitative approach for understanding the informativeness of
experimental characterizations about parameters of BCP self-assembly models.

  3. 3.

We highlight the computational challenges of integrated likelihood evaluations
for calibrating a BCP self-assembly model and demonstrate that a LFI approach
is suitable for overcoming such challenges.

  4. 4.

We overview one methodology for LFI, based on measure transport, that is well-
suited for inference problems with generative models, given that it only
requires a limited number of offline realizations of parameters and model
evaluations.

  5. 5.

We introduce the notion of extracting summary statistics from BCP
characterization image data and propose an efficient approach to measure the
utility of selected summary statistics in the LFI setting via a noval EIG
estimator that is, in principle, consistent.

  6. 6.

We present a numerical case study on the Bayesian calibration and EIG
computation of the OK model for Di-BCP thin film self-assembly with top-down
microscopy characterization.

##  2 A probabilistic data model of block copolymer self-assembly

In this section, we conceptualize and formulate generic probabilistic data
models of BCP self-assembly that consist of (i) statistical models of
experimental characterization, introduced in Section 2.1, and (ii) statistical
models of BCP self-assembly, introduced in Section 2.2. The former account for
measurement uncertainty in the characterization data, while the latter account
for the aleatoric uncertainty. A particular example of the data model that is
used for a numerical case study is provided in Section 7.

###  2.1 Statistical models of experimental characterization

We consider the process that generates an image, $\bm{y}\in\mathbb{R}^{y}$,
that characterizes a latent BCP melt $u\in V^{u}(\mathcal{D})$ in a physical
domain $\mathcal{D}\subset\mathbb{R}^{3}$. The material state $u$ may be
described by discrete spatial positions of monomers or coarse-grained monomer
segments. On the meso- or continuum scale, the material state can also be
modeled as the local fractions or normalized densities fields of distinct
monomers. We denote the set of all possible Di-BCP material states as
$V^{u}(\mathcal{D})$.

Let the characterization process refer to the map from $u\mapsto\bm{y}$. We
model this process as the composition of a state-to-observable map
$\mathcal{J}$ and a statistical sensor operator $\bm{\mathcal{S}}$,
respectively representing signal formation and image formation in the
characterization process. That is,

$\bm{y}=\bm{\mathcal{S}}(\mathcal{J}(u),\bm{n}),$ (1)

where $\bm{n}$ represents the stochastic noise introduced during the
characterization process.

|
---|---
Figure 2: (_left_) A schematic drawing of the GISAXS experimental setup for a
Di-BCP film with perpendicular lamellae. Reprinted (adapted) with permission
from [46]. Copyright ©2007 American Chemical Society. (_right_) A schematic
drawing of res-CDSAXS experimental setup and a reconstructed scattering image
of PS–b–PMMA film strongly guided by a chemically nanopatterned substrate
[40].

The state-to-observable map can be further decomposed into a restriction
operator $\mathcal{R}$ and an imaging operator $\mathcal{I}$, i.e.,
$\mathcal{J}\coloneqq\mathcal{I}\circ\mathcal{R}$. The restriction
operator777In many Bayesian inference problems associated with real-space
observations of a physical domain, this restriction operator is a point-wise
observation operator. See [47]. $\mathcal{R}\colon V^{u}(\mathcal{D})\to
V^{u}\big{(}R(\mathcal{D}))$ extracts the observed portion of the material
state in the subdomain $R(\mathcal{D})\subset\mathcal{D}$. For microscopy
characterization, $R(\mathcal{D})$ may represent regions of top or cross-
section surfaces in the field of view of the microscope. For X-ray scattering
characterization, it may represent the subdomain around the incident X-ray
beam that affects its trajectory.

The operator $\mathcal{I}\colon V^{u}\big{(}R(\mathcal{D})\big{)}\to V^{o}$
models the imaging mechanism, which depends on the experimental setup and
certain properties of the latent BCP material, such as the refractive indices
of the distinct monomers that are relevant for X-ray scattering imaging [48].
Most microscopy imaging methods are real-space methods, where the imaging
operator is often constructed using identity maps. One may also include
effects related to image resolution, contrast, brightness, depth of field,
etc., in the microscopy imaging operator [49]. For X-ray scattering
characterization, one may additionally incorporate transformations to the
reciprocal space (e.g. via Fourier transform) [50, 39, 51, 52] in the imaging
operator for reproducing scattering patterns from given material states.

The statistical sensor operator $\bm{\mathcal{S}}\colon
V^{o}\times\mathbb{R}^{n}\to\mathbb{R}^{y}$ models the discretization and
noise corruptions introduced by image-forming sensors. The sensor noise
$\bm{n}\in\mathbb{R}^{n}$ is a realization of the random variable $\bm{N}$,
often attributed to infinitesimal fluctuations in sensor voltage. A simple
additive model with independent noise is often considered:

$\bm{\mathcal{S}}\left(\mathcal{J}(u),\bm{n}\right)=\bm{\mathcal{D}}\left(\mathcal{J}(u)\right)+\bm{n},$
(2)

where the components of the stochastic noise are independent and identically
distributed (i.i.d.), and $\bm{\mathcal{D}}\colon V^{o}\to\mathbb{R}^{y}$
represents the local averaging of signals introduced by image-forming sensors.
Let us note that realistic sensor noise models also account for more complex
spatial and signal intensity dependence [53, 54, 55].

###  2.2 The aleatoric uncertainty and statistical models of block copolymer
self-assembly

While the self-assembly of block copolymers leads to the formation of ordered
periodic structures, such as those depicted in the experimental phase diagram
in Fig. 1 (_left_), they are typically not realizable in experiments without
strong external guidance, even though they are often considered to be the most
stable configurations. Instead, a rich variety of states that exhibit long-
range disorder, often referred to as defected states or metasable states are
observed; see Fig. 1 (_right_).

To account for the impact of aleatoric uncertainty on the experimentally
observed equilibrium states of BCP self-assembly, we assume that the various
states attainable under the same experimental scenario $S_{E}$ to be
realizations of a $V^{u}$-valued random variable, $U(S_{E})$. We refer to
$U(S_{E})$ as the state random variable. The state random variable gives rise
to a probability distribution of the achievable states, whose structure
depends on the material, geometrical, thermal, and various other controllable
macroscopic factors specified by $S_{E}$. When the self-assembly is strongly
guided, the distribution of the state random variable is close to a Dirac
mass, meaning that only a handful of BCP equilibrium states are attainable.
When the self-assembly is unguided or weakly guided, the distribution can be
quite complex, yet it typically concentrates on states in the same phase or
mixture of phases that have similar structures, such as similar radius and
periodicity length of spheres in the sphere phase.

A statistical model of BCP self-assembly can be generalized as a mathematical
model derived from statistical polymer theory for specifying the probability
distribution, i.e., the statistical ensemble, of $U(S_{E})$ for a range of
experimental scenarios. It typically contains a model of polymer chains,
models of the inter-polymer and polymer–environment interactions, and possibly
a model of thermal fluctuations and/or dynamics of polymers during the self-
assembly. A comprehensive review of these models can be found in [29, Vol. 1],
[31, 56, 19], and references therein.

Computer simulations of BCP self-assembly can be treated as stochastic
simulators for generating independent samples from model-specific statistical
ensembles of equilibrium states. Let us denote the forward operator of a
statistical model of BCP self-assembly as
$\mathcal{F}\colon\mathbb{R}^{x}\times V^{z}\to V^{u}$. This model assumes the
equality in distribution

$U(S_{E})\stackrel{{\scriptstyle d}}{{=}}\mathcal{F}(\bm{x},Z),$ (3)

where $\bm{x}\in\mathbb{R}^{x}$ are parameters888This work considers the model
parameters $\bm{x}$ to be finite dimensional and belonging to a Euclidean
space. See [57] for a general treatment in which $\bm{x}$ belongs to an
infinite-dimensional function space. that represents the material,
geometrical, thermal, and/or other macroscopic factors of the experimental
scenario $S_{E}$ that are considered in the model. The $V^{z}$-valued random
variable $Z$ represents the stochastic component of the forward operator.

The origin of the stochastic component $Z$ may be purely numerical. For
example, $Z$ can partly represent random initial configurations of the
polymeric system and the sampling mechanism for accepting perturbed
configurations in Monte Carlo simulations of equilibrium BCP melts [58]. It
can also be high-dimensional or even infinite-dimensional. For example, $Z$
can be specified by the stochastic forcing term that resides in infinite-
dimensional function spaces for dynamical DFT models [20]. We refer to $Z$ as
an auxiliary random variable, or nuisance variable, as it is either non-
physical or represents uncontrollable microscopic factors in experiments that
are difficult to observe. For the numerical case study in Section 7 and 8, we
use the OK model, in which the stochastic component is specified by random
initial states for free energy minimization.

###  2.3 A generic data model of block copolymer self-assembly

We call the combination of the model of BCP self-assembly and model of
experimental characterizations of equilibrium states as the _data model_. For
notational convenience, we first define the parameter-to-observable map
$\mathcal{O}\coloneqq\mathcal{J}\circ\mathcal{F}$, where
$\mathcal{O}(\bm{X},Z)$ is the observable random variable. For example, for
X-ray scattering characterization, $\mathcal{O}(\bm{x},Z)$ may be the (noise-
free) model-predicted latent scattering pattern where variability comes from
the model-predicted random BCP equilibrium structures. We define the image
data as a $\mathbb{R}^{y}$-valued random variable arising from both the
measurement uncertainty and the aleatoric uncertainty. The data model
$\bm{\mathcal{G}}:\mathbb{R}^{x}\times
V^{z}\times\mathbb{R}^{n}\to\mathbb{R}^{y}$ summarizes the relationship
between model parameters and data as

$\bm{Y}=\bm{\mathcal{G}}(\bm{x},Z,\bm{N})\coloneqq\bm{\mathcal{S}}\big{(}\mathcal{O}(\bm{x},Z),\bm{N}\big{)}.$
(4)

##  3 Bayesian inference and quantifying uncertainty

In this section, we present the Bayesian inference framework for calibrating
models for BCP self-assembly. In Section 3.1 we introduce the posterior
distributions that arise from the likelihood and the, so called, conditional
likelihood. To summarize the uncertainty in calibrated parameters, we
introduce some information-theoretic quantities in Section 3.2. Lastly, in
Section 3.3, we present the expected information gain as a natural quantity to
consider in a Bayesian inference context for measuring the information
contained in the data about the unknown model parameters.

###  3.1 The posterior and joint posterior

For a Bayesian formulation to model calibration, we introduce a random
variable for the parameters $\bm{X}\in\mathbb{R}^{x}$ and seek to characterize
its distribution given an observation of $\bm{Y}$. Initially, the parameter
$\bm{X}$ is distributed according to the _prior_ measure
$\nu_{\bm{X}}\in\mathcal{P}(\mathbb{R}^{x})$ with density $\pi_{\bm{X}}$. The
prior measure represents our _epistemic uncertainty_ in the model parameters
before collecting observations. Similarly, the (possibly infinite-dimensional)
auxiliary random variable is distributed according to the conditional prior
measure $\nu_{Z|\bm{X}}$ that may depend on the model parameter. In our
numerical case study, a prior distribution for the OK model parameters is
constructed in Section 7.5 and the distribution for the auxiliary variable is
described in A.

The data model in Section 2.3 defines the conditional distribution for the
observations. Let $\bm{Y}=\bm{\mathcal{G}}(\bm{X},Z,\bm{N})$ be a data random
variable and assume the sensor noise random variable $\bm{N}\sim\pi_{\bm{N}}$
is independent of ($\bm{X}$, $Z$). This noise yields a conditional density for
the data given the inputs $\bm{x}$ and $z$ that we denote by
$\pi_{\bm{Y}|\bm{X},Z}(\cdot|\bm{x},z)$. For example, if the sensor operator
in the data model (4) is given by the additive form in (2), then this
conditional density is given by
$\pi_{\bm{N}}(\cdot-\bm{\mathcal{D}}(\mathcal{O}(\bm{x},z)))$.

For the BCP model calibration problem, however, we are not interested in
inferring the auxiliary variable $Z$. Thus, we seek to describe the
conditional distribution of the data given only the model parameters. The
uncertainties in the sensor noise and the auxiliary variable define a marginal
conditional distribution for $\bm{Y}$ given $\bm{X}$, whose density we denote
by $\pi_{\bm{Y}|\bm{X}}(\cdot|\bm{x})$. For these conditional distributions of
the data, we define both the _conditional likelihood_ and the _likelihood_ for
an observation $\bm{y}^{*}$ as the functions which arise when fixing the data
to $\bm{y}^{*}$ in the conditional densities. That is,
$\pi_{\bm{y}^{*}|\bm{X},Z}(\bm{x},z)\coloneqq\pi_{\bm{Y}|\bm{X},Z}(\bm{y}^{*}|\bm{x},z)$
is the conditional likelihood for both the model parameters and the auxiliary
variable, while
$\pi_{\bm{y}^{*}|\bm{X}}(\bm{x})\coloneqq\pi_{\bm{Y}|\bm{X}}(\bm{y}^{*}|\bm{x})$
is the likelihood for only the model parameters. The two likelihood functions
represent, respectively, the likelihood of observing an image $\bm{y}^{*}$ for
(i) a particular model-predicted equilibrium state generated at $(\bm{x},z)$,
and (ii) a variety of random model-predicted equilibrium states generated at
$\bm{x}$. We will discuss how both are defined in depth in Section 4.

Given a prior and likelihood function for an observation $\bm{y}^{*}$, Bayes’
theorem provides a rule for producing the conditional distributions of the
parameters. The _joint posterior_ density for the parameters and the auxiliary
variable is given by999Similar to the notation of the likelihood function, we
use $\pi_{\bm{X},Z|\bm{y}^{*}}(\bm{x},z)$ to emphasize that the posterior is
conditioned on a particular observation. For the more general case of
arbitrary observations $\bm{y}\in\mathbb{R}^{y}$, we use the notation
$\pi_{\bm{X},Z|\bm{Y}}(\bm{x},z|\bm{y})$ instead. Note that in the setting
where the auxiliary variable $Z$ is infinite dimensional, Bayes’ rule in (5)
is no longer mathematically valid and it should be expressed using the
Radon–Nikodym derivative between the posterior and the prior [57].

$\displaystyle\pi_{\bm{X},Z|\bm{y}^{*}}(\bm{x},z)=\frac{\pi_{\bm{y}^{*}|\bm{X},Z}(\bm{x},z)\pi_{\bm{X},Z}(\bm{x},z)}{C_{\bm{y}^{*}}},$
(5)

where
$C_{\bm{y}^{*}}=\pi_{\bm{Y}}(\bm{y}^{*})=\int_{\mathbb{R}^{x}}\int_{V^{z}}\pi_{\bm{y}^{*}|\bm{X},Z}(\bm{x},z)\textrm{d}\nu_{Z,\bm{X}}(z|\bm{x})\pi_{\bm{X}}(\bm{x})\textrm{d}\bm{x}>0$
is the normalization constant for the posterior density and is often referred
to as the _model evidence_ or the _marginal likelihood_. To calibrate the
model parameters, our interest is in the _marginal posterior_ density of the
parameters alone given by

$\displaystyle\pi_{\bm{X}|\bm{y}^{*}}(\bm{x})=\frac{\pi_{\bm{y}^{*}|\bm{X}}(\bm{x})\pi_{\bm{X}}(\bm{x})}{C_{\bm{y}^{*}}}.$
(6)

Many computational tasks in Bayesian inference rely on evaluating these
densities to draw samples from the posterior distributions. While evaluations
of (5) depend on evaluations of the conditional likelihood
$\pi_{\bm{y}^{*}|\bm{X},Z}$, evaluations of (6) inherit difficulties from
evaluating the likelihood $\pi_{\bm{y}^{*}|\bm{X}}$ that stem from the
influence of the nuisance variable $Z$ on the observations. We discuss these
computational difficulties in Section 4. For a thorough treatment on these
challenges in Bayesian inference, in general, see [59, 60].

###  3.2 Entropy and information gain

To understand and quantify the effectiveness of Bayesian model calibration, we
turn to information theory that formally quantifies the information content in
random variables. One information-theoretic measure for the uncertainty of a
parameter $\bm{X}$ is the _Shannon entropy_
$\mathrm{H}\colon\mathcal{P}(\mathbb{R}^{x})\to\mathbb{R}$ given by

$\mathrm{H}(\bm{X})\coloneqq\mathbb{E}_{\pi_{\bm{X}}}[-\log\pi_{\bm{X}}].$ (7)

The entropy is larger for a parameter that is less certain. For example, if a
parameter is known, i.e., its distribution is a Dirac mass, then its entropy
is zero. After conditioning on the data $\bm{y}^{*}$, we can measure the
entropy of the parameters by computing the conditional entropy
$\mathrm{H}(\bm{X}|\bm{Y}=\bm{y}^{*})$. If the entropy decreases after
conditioning, i.e., $\mathrm{H}(\bm{X}|\bm{Y}=\bm{y}^{*})<\mathrm{H}(\bm{X})$,
we say that observing the data $\bm{y}^{*}$ has _decreased our uncertainty_ or
_increased our confidence_ in the values of model parameters. The expectation
of $\mathrm{H}(\bm{X}|\bm{Y}=\bm{y})$ over the distribution of the data random
variable defines the conditional entropy of the parameters given the data
$\mathrm{H}(\bm{X}|\bm{Y})$ and it always satisfies
$\mathrm{H}(\bm{X}|\bm{Y})\leq\mathrm{H}(\bm{X})$. If the conditional entropy
is _less_ than the prior’s entropy, we say that the data is, on average,
informative of the parameters.

To quantify the information gained from performing Bayesian inference, we
measure the change in the entropy of the parameters $\bm{X}$ before and after
conditioning on data. This change quantifies the value of acquiring certain
noisy, corrupted, and/or partial observations of the system to learn about the
unknown model parameters. In the context of Bayesian inference, this change
compares the prior density $\pi_{\bm{X}}$ and the posterior density
$\pi_{\bm{X}|\bm{y}^{*}}$ for the parameters. An information-theoretic
quantity for measuring information gain based on the change in entropy is the
Kullback–Leibler (KL) divergence $\mathcal{D}_{\textrm{KL}}$. The KL
divergence from the prior to the posterior is given by

$\mathcal{D}_{\textrm{KL}}(\pi_{\bm{X}|\bm{y}^{*}}||\pi_{\bm{X}})=\mathbb{E}_{\pi_{\bm{X}|\bm{y}^{*}}}\left[\log\frac{\pi_{\bm{X}|\bm{y}^{*}}}{\pi_{\bm{X}}}\right].$
(8)

The KL divergence is always non-negative and is zero if and only if the
posterior is equal to the prior almost everywhere. Furthermore, a large value
for the KL divergence implies that the posterior has lower entropy relative to
the prior and its density is typically more concentrated than the prior
density. A large departure of the posterior from the prior indicates that we
gained more information from the observations about the parameters.

###  3.3 Expected information gain

Before observing any particular data, we can measure the information gained
about the parameter in a Bayesian inference problem by considering the
expected information gain (EIG)
$I\colon\mathcal{P}(\mathbb{R}^{x})\times\mathcal{P}(\mathbb{R}^{y})\to\mathbb{R}_{+}$
between the random variables $\bm{X}$ and $\bm{Y}$. We say that
$I(\bm{X};\bm{Y})$ represents the the EIG for learning about $\bm{X}$ from
$\bm{Y}$. The EIG can be defined by the integration of the KL divergence in
(8) with respect the data $\bm{y}$. That is,

$I(\bm{X};\bm{Y})=\int_{\mathbb{R}^{y}}\mathcal{D}_{\textrm{KL}}(\pi_{\bm{X}|\bm{y}}||\pi_{\bm{X}})\pi_{\bm{Y}}(\bm{y})\textrm{d}\bm{y}.$
(9)

Within the field of Bayesian optimal experimental design (BOED), the EIG is
commonly used to assess the utility of an experimental measurement for
reducing uncertainty in the model parameters [61]. In this context, the EIG is
typically a function of the experimental design variables, and is maximized to
find the optimal design. An example of BOED is the optimal design of a
discrete observation operator, e.g., choosing sensor locations that maximize
the EIG on certain quantities of interest [62]. In the context of microscopy
characterization of BCP melts, for example, we can measure $I(\bm{X};\bm{Y})$
for different restriction operators $\mathcal{R}$, introduced in Section 2.1,
to determine the optimal level of magnification for inferring a particular
subset of the model parameters. In our numerical case study in Section 7 and
8, we will demonstrate results on the computation of EIG that enables such a
type of BOED study.

We note that the expression for the EIG $I(\bm{X};\bm{Y})$ in (9) is
equivalent to the mutual information (MI) between random variables $\bm{X}$
and $\bm{Y}$. The MI is always a non-negative random variable. While a zero
value implies that the two random variables are independent, a larger value
indicates that the data is more informative on average for inferring the
parameter. Using the form of the KL divergence in (8) and Bayes’ rule in (6),
the EIG can also be written in the two equivalent ways

$\displaystyle I(\bm{X};\bm{Y})$
$\displaystyle=\mathbb{E}_{\nu_{\bm{X},\bm{Y}}}\left[\log\frac{\pi_{\bm{X}|\bm{Y}}}{\pi_{\bm{X}}}\right]$
(10a)
$\displaystyle=\mathbb{E}_{\nu_{\bm{X},\bm{Y}}}\left[\log\frac{\pi_{\bm{Y}|\bm{X}}}{\pi_{\bm{Y}}}\right].$
(10b)

While the first expression uses the ratio of the posterior and prior
densities, the second expression uses the ratio of the conditional likelihood
and the marginal density for the observations. We note that as compared to
standard Bayesian inference problems that characterize the posterior for a
single realization of the data, computing the EIG requires evaluating the
posterior densities for many possible realizations of the data. In Section 6,
we propose an efficient estimator for the EIG in (10a) that is based on a
variational approximation for all of the posterior densities.

##  4 Addressing the challenges of auxiliary variables and high-dimensional
image data via likelihood-free inference

Performing Bayesian calibration and computing EIG requires _characterizing
posteriors_ , i.e., producing representative samples from the the posterior
distribution or providing an exact or approximate expression of the density
$\pi_{\bm{X}|\bm{y}^{*}}$. In this section, we outline distinctive challenges
for the Bayesian calibration of BCP self-assembly models. In Section 4.1 we
describe intractable integrated likelihood functions that arise when
accounting for the aleatoric uncertainty in the data model. Sections 4.2 and
4.3 discuss computational methods for dealing with integrated likelihoods,
including LFI methods. Lastly, in Section 4.4 we address inference with high-
dimensional image data by presenting an approach that constructs low-
dimensional summary statistics of the data to capture relevant information
about the model parameters of interest.

###  4.1 Evaluating integrated likelihoods

To calibrate the model parameters $\bm{x}$ given image data $\bm{y}^{*}$, it
is necessary to prescribe the likelihood function $\pi_{\bm{y}^{*}|\bm{X}}$ in
(6). This likelihood is derived from the stochastic procedure for generating
the image data using the data model for BCP self-assembly in (4). The
procedure depends on, not only the sensor noise $\bm{N}$, but also the
auxiliary variable $Z$ entering the parameter-to-observable map. In
particular, the likelihood of an observed image $\bm{y}^{*}\in\mathbb{R}^{y}$
given model parameters $\bm{x}$ accounts for _all_ realizations of the
auxiliary variable and stochastic noise that could possibly generate the
particular image. Mathematically, this is achieved by summing the
contributions of the conditional likelihoods
$\pi_{\bm{y}^{*}|\bm{X},Z}(\bm{x},z)$ weighted by the probability of
generating each auxiliary variable $z$. Hence, the likelihood function for the
parameters is defined as the following _integrated likelihood_ [36] that
marginalizes out the auxiliary variables

$\pi_{\bm{y}^{*}|\bm{X}}(\bm{x})=\int_{V^{z}}\pi_{\bm{y}^{*}|\bm{X},Z}(\bm{x},z)\textrm{d}\nu_{Z|\bm{X}}(z|\bm{x}).$
(11)

Using this likelihood leads to the posterior density for the parameters
$\bm{x}$ in (6) that is a low-dimensional marginal of the joint posterior
distribution for $(\bm{X},Z)$101010Characterizing the joint posterior involves
a tractable likelihood $\pi_{\bm{y}^{*}|\bm{X},Z}$ that does not require
marginalization, but requires performing inference on a potentially much
higher-dimensional distribution..

The integrated likelihood in (11) has a closed-form expression only in very
special situations. One setting is when both the conditional distribution for
the auxiliary variables $\nu_{Z|\bm{X}}(\cdot|\bm{x})$ and the conditional
density for the observations $\pi_{\bm{Y}|\bm{X},Z}(\cdot|\bm{x},z)$ are
Gaussian, and the parameter-to-observable map is affine. In general, however,
the integral we consider over the infinite-dimensional auxiliary variables $z$
is unavailable in closed form111111The conditional likelihood is an example of
a _semi-parametric_ model that contains both finite-dimensional and infinite-
dimensional inputs; our work seeks to build on inference for semi-parametric
models, where the infinite-dimensional auxiliary variables can include
stochastic forcing terms, stochastic initial conditions or random fields,
which are not relevant in inference. We refer the reader to [63] for more
details on these models.. In these cases, we say that the integrated
likelihood is _unavailable_ and consider it _intractable_ to evaluate.

Despite their intractability, integrated likelihoods are frequently
encountered in many fields. Other interesting inference problems that involve
integrated likelihoods include time, phase, and distance marginalization in
gravitational wave inference [64], boundary condition marginalization for
linear parabolic PDEs [65], and hyper-parameter discovery for hierarchical
models in image deblurring and X–ray computational tomography [66]. Therefore,
developing efficient computational inference methods that work with integrated
likelihoods is an important ongoing area of research.

###  4.2 Inference methods for intractable likelihoods

In settings with intractable likelihoods, it is not possible to compute
posteriors using likelihood-based sampling algorithms, such Markov chain Monte
Carlo (MCMC), or other variational methods which require evaluating the
likelihood at multiple values of the parameters. To satisfy this constraint,
three classes of inference methods have been developed that do not require
explicit likelihood evaluations. These classes include surrogate likelihood
functions, pseudo-marginal Monte Carlo methods [67, 68], which assume
tractable conditional likelihood function evaluations, and likelihood-free
inference, often referred to as simulation-based inference [69].

The _profile likelihood_ method [70] is one popular choice of surrogate that
replaces integration in (11) with a maximization over $z$. Specifically, the
profile likelihood $\widehat{\pi}_{\bm{y}^{*}|\bm{X}}$ is defined as the
conditional likelihood evaluated at the maximum likelihood estimate (MLE) for
$z$ as a function of $\bm{x}$, i.e.,
$z_{\textrm{MLE}}(\bm{x})=\operatorname*{arg\,max}_{z\in
V^{z}}\pi_{\bm{y}^{*}|\bm{X},Z}(\bm{x},z)$. Thus, if a well-defined unique
$z_{\textrm{MLE}}$ exists, the profile likelihood can be tractably evaluated
as
$\widehat{\pi}_{\bm{y}^{*}|\bm{X}}(\bm{x})\coloneqq\pi_{\bm{y}^{*}|\bm{X},Z}(\bm{x},z_{\textrm{MLE}}(\bm{x}))$.
This approach, however, replaces $z$ by a point estimate and ignores the
impact of its uncertainty on the likelihood. As a result, the profile
likelihood is inconsistent and often requires adding factors to correctly
characterize the posterior uncertainty [71].

If conditional likelihood evaluations are available, pseudo-marginal methods
can perform inference using Monte Carlo estimates of the integrated likelihood
function in (11). These approaches generate $N$ samples of the auxiliary
variable $\\{z^{i}\\}_{i=1}^{N}\sim\nu_{Z|\bm{X}}(\cdot|\bm{x})$ and estimate
the integrated likelihood as
$\widehat{\pi}_{\bm{y}^{*}|\bm{X}}(\bm{x})\approx\frac{1}{N}\sum_{i=1}^{N}\pi_{\bm{y}^{*}|\bm{X},Z}(\bm{x},z^{i})$.
If the variance
$\mathbb{V}\textrm{ar}(\pi_{\bm{y}^{*}|\bm{X},Z}(\bm{x},\cdot))$ is finite,
the standard deviation of this estimator decreases with more samples at the
rate $\mathcal{O}(1/\sqrt{N})$. Thus, we can reduce the sample variance of
$\widehat{\pi}_{\bm{y}^{*}|\bm{X}}$ by increasing the number of samples until
sufficient accuracy of the integration has been reached. While unbiased Monte
Carlo estimators of the integrated likelihood have been shown to yield
consistent posterior sampling [67], these methods may necessitate a large
number of samples to achieve a reasonable computational efficiency in practice
[72]. On the other hand, the evaluation of the conditional likelihood function
at each sample of the auxiliary variable requires evaluating the forward
operator, which can be computationally expensive. Furthermore, estimates of
the integrated likelihood are not shared for different parameter values and
need to be recomputed at each iteration of some sampling methods, such as
MCMC.

In addition to the challenge of marginalization, the conditional likelihood
$\pi_{\bm{y}^{*}|\bm{X},Z}$ can also become analytically intractable to
evaluate for many non-trivial sensor noise models,
$\bm{\mathcal{S}}(\cdot,\bm{N})$, which rules out the use of pseudo-marginal
methods and the profile likelihood. The analytical expression for the
conditional likelihood stems from a change of variables between the noise
random variable and the data random variable. For example, the sensor operator
with multiplicative noise, i.e.,
$\bm{\mathcal{S}}(\cdot,\bm{N})=\bm{\mathcal{D}}(\cdot)\bm{N}$, with noise
$\bm{N}\sim\pi_{\bm{N}}$ has the conditional likelihood

$\displaystyle\pi_{\bm{y}^{*}|\bm{X},Z}(\bm{x},z)=\pi_{\bm{N}}\Big{(}\bm{y}^{*}/\bm{\mathcal{D}}\big{(}\mathcal{O}(\bm{x},z)\big{)}\Big{)}\Big{(}\textrm{det}\big{|}\bm{\mathcal{D}}(\mathcal{O}(\bm{x},z))\big{|}\Big{)}^{-1}\;.$

In general, evaluating the conditional likelihood requires computing the
inverse and Jacobian inverse of the sensor noise model, which is often
unfeasible when it has a complicated form121212We note that,
$\pi_{\bm{Y}|\bm{X},Z}(\bm{y}|\bm{x},z)=\sum_{A_{k}}\mathds{1}_{\mathcal{S}_{\bm{x},z}(A_{k})}(\bm{y})\pi_{\bm{N}}(\mathcal{S}_{\bm{x},z}^{-1}(\bm{y}))\det\big{|}\nabla\mathcal{S}_{\bm{x},z}^{-1}(\bm{y})\big{|}$
where
$\mathcal{S}_{\bm{x},z}\coloneqq\bm{\mathcal{S}}(\mathcal{O}(\bm{x},z),\cdot):\mathbb{R}^{n}\to\mathbb{R}^{y}$
is presumed to be a generally non-monotonic differentiable transformation, and
$\pi_{\bm{Y}|\bm{X},Z}(\bm{y}|\bm{x},z)$ is a weighted sum of densities over a
partition of the input space $A_{k}\subset\mathbb{R}^{n}$ into regions for
which $\mathcal{S}_{\bm{x},z}$ is monotonic. Unless this reduces to a simpler
closed–form expression with explicit computations, the evaluation of
$\pi_{\bm{Y}|\bm{X},Z}(\bm{y}|\bm{x},z)$ can be quite complex for arbitrary
$\bm{\mathcal{S}}$, due to the difficulties of evaluating generalized inverses
$\mathcal{S}_{\bm{x},z}^{-1}$, computing Jacobian determinants (generally
$O(n^{3})$ computational complexity, where $n\gg 1$ for our sensor noise), and
determining the partitions $\\{A_{k}\\}$.. When the data arises from a “black-
box” simulator, i.e., when the data model $\bm{\mathcal{G}}$ relies on legacy
scientific code whose mathematical expressions may be unavailable, it is not
possible to derive the conditional likelihood analytically. In polymer
modeling and characterization, complex mathematical expressions and scientific
codes underlay both the forward operators $\mathcal{F}$ and realistic sensor
operators $\bm{\mathcal{S}}$. Therefore, we consider the evaluation of
conditional likelihoods to be generally unfeasible.

On the other hand, likelihood-free inference (LFI) methods enable statistical
inference in generative models. They rely _only_ on sampling from the joint
distribution of model parameters and model-generated data. In the context of
the integrated likelihood function in (11), this is accomplished by

  1. 1.

sampling the model parameters, $\bm{x}^{i}\sim\pi_{\bm{X}}$,

  2. 2.

sampling the auxiliary variable, $z^{i}\sim\nu_{Z|\bm{X}}(\cdot|\bm{x}^{i})$,

  3. 3.

evaluating the parameter–to–observable map,
$o^{i}=\mathcal{O}(\bm{x}^{i},z^{i})$,

  4. 4.

sampling the sensor noise $\bm{n}^{i}\sim\pi_{\bm{N}}$, and

  5. 5.

evaluating the statistical sensor model,
$\bm{y}^{i}=\bm{\mathcal{S}}(o^{i},\bm{n}^{i})$.

The resulting samples $(\bm{x}^{i},\bm{y}^{i})$, with $z^{i}$ and $\bm{n}^{i}$
discarded, are samples from the marginal distribution of $(\bm{X},\bm{Y})$.
This collection of marginal samples avoids the step of explicitly
marginalizing the conditional likelihood $\pi_{\bm{Y}|\bm{X},Z}$ to evaluate
the integrated likelihood in (11).

Given a set of i.i.d. joint samples $(\bm{x}^{i},\bm{y}^{i})$ that is large
enough to be representative of their joint distribution, the goal of LFI
methods is to approximate the posterior density $\pi_{\bm{X}|\bm{y}^{*}}$ or
to sample from it, without evaluating the likelihood function. Thus, LFI is
agnostic to the form of the likelihood and applicable for inference with
various sources of intractability, including the high-dimensional integration
with respect to the auxiliary variable. This makes LFI appealing for the BCP
model calibration problem.

###  4.3 Likelihood-free inference: Related work and practical desiderata

The family of LFI methods for posterior sampling includes both classical
techniques such as approximate Bayesian computation (ABC) and synthetic
likelihoods, as well as modern machine learning methods. ABC methods [73, 74,
75] define a non-parametric approximation to the likelihood function, and
thereby the posterior, using an accept/reject sampling scheme. For each pair
of joint samples $(\bm{x}^{i},\bm{y}^{i})$, this scheme accepts parameters
whose simulated data $\bm{y}^{i}$ is close to the observed data $\bm{y}^{*}$
with respect to a discrepancy function $d$ and tolerance level $\epsilon$,
i.e., ABC accepts $\bm{x}^{i}$ if $d(\bm{y}^{i},\bm{y}^{*})<\epsilon$. While
ABC can be shown to be a consistent approximation of the posterior [76] in the
limit as $\epsilon\rightarrow 0$, small values for $\epsilon$ lead to
rejection of many candidate parameters, and often an impractical number of
simulations is required to accurately estimate posterior statistics.
Furthermore, the quality of the non-parametric approximation in ABC
deteriorates as the data dimension increases [77].

An alternative to ABC that is biased, but typically more efficient, is a
parametric approximation to the likelihood function. As an example, the
Bayesian synthetic likelihood (BSL) method uses samples
$\bm{y}^{i}\sim\nu_{\bm{Y}|\bm{X}}(\cdot|\bm{x})$ to create a multivariate
Gaussian approximation to the likelihood at each parameter value $\bm{x}$
given by
$\widehat{\pi}_{\bm{Y}|\bm{X}}=\mathcal{N}(\overline{\bm{y}}(\bm{x}),C(\bm{x}))$
where $\overline{\bm{y}}(\bm{x})$ and $C(\bm{x})$ denote the estimated mean
and covariance matrix. The approximate likelihood evaluations
$\widehat{\pi}_{\bm{y}^{*}|\bm{x}}(\bm{x})=\widehat{\pi}_{\bm{Y}|\bm{x}}(\bm{y}^{*})$
can then be efficiently used within MCMC algorithms [78, 79]. While BSL can be
more efficient than ABC for high-dimensional data, this approach requires
sampling multiple data, i.e., evaluating the data model, for each parameter in
order to determine $\overline{\bm{y}}(\bm{x})$ and $C(\bm{x})$. As a result,
both ABC and synthetic likelihood approaches are less suitable for
applications that involve computationally expensive forward operators. Lastly,
the inference results of both methods are applicable only to the particular
set of observed data $\bm{y}^{*}$ and thus require repeating the computation
for any new set of data.

To amortize the cost of performing inference across many sets of observed
data, recent machine learning methods propose to learn models that are
functions of both the parameters and the data [69]. These methods use the
collection of joint samples $(\bm{x}^{i},\bm{y}^{i})$ to learn variational
approximations for either the posterior density $\pi_{\bm{X}|\bm{Y}}$ [37,
80], the likelihood function $\pi_{\bm{Y}|\bm{X}}$ [81], or the likelihood
ratio $\pi_{\bm{X}|\bm{Y}}/\pi_{\bm{X}}$ [82]. The resulting models are then
used in an inference step to define a tractable approximation to the posterior
density or to generate approximate posterior samples for any observed data
$\bm{y}^{*}$. This enables the inference results to be used for downstream
tasks that require access to the posteriors for all data realizations, such as
EIG computation and BOED. To avoid constraining the structural forms for these
models, many techniques use flexible neural network architectures that can
capture complex posterior distributions and can learn the structure of high-
dimensional and domain-specific data (e.g., using architectures designed for
images). Many of these models, however, often require a large upfront cost to
learn them reliably and adapting these methods to small-sample settings
remains challenging. To select a LFI method for the calibration of BCP self-
assembly models, we have the following desiderata:

  1. 1.

Sample efficiency: Method provides reasonable posterior approximations in
settings where it is computationally expensive to generate simulated data.

  2. 2.

Consistency: Method does not impose arbitrary approximations on the form of
the likelihood or the posterior distribution in order to capture non-
Gaussianity without biases from the inference procedure.

  3. 3.

Amortized inference: Method provides posterior densities
$\pi_{\bm{X}|\bm{y}^{*}}$ or samples for almost _all_ values of the data
without needing to repeat inference for each realization of $\bm{y}^{*}$.

To produce _consistent_ posterior approximations, without discarding or
rejecting samples, we consider ABC and BSL to be less appealing for our
application. To also benefit from amortization similarly to the machine
learning methods, we propose to use a measure transport approach for LFI in
this work, which relies on a careful parameterization for settings with
limited sample sizes. The details of the proposed approach are discussed in
Section 5.

###  4.4 Informative summary statistics of high-dimensional image data

When performing LFI with high-dimensional data, it is common to construct and
work with low-dimensional summary statistics of the data [83],
$\bm{Q}=\mathcal{Q}(\bm{Y})\in\mathbb{R}^{q}$. Working with this type of
summary statistics has two inherent advantages: (i) it can improve the quality
and tractability of LFI, and (ii) it can provide insight into what information
contained in the data is relevant for inference.

In principle, arbitrarily high-dimensional data could be used to perform
inference. This data, however, can lead to concentrated regions of the
parameter space with large likelihood values, which hampers the effectiveness
of many LFI methods. Here we provide two numerical examples to visualize this
phenomenon. For both examples, we consider the calibration of the OK model for
Di-BCP thin film self-assembly using top-down microscopy image data. We
consider the characterization model with additive Gaussian noise
$\bm{N}\sim\mathcal{N}(0,C_{\bm{N}})$ as described in (2). That is

$\bm{Y}=\bm{\mathcal{D}}\left(\mathcal{F}(\bm{X},Z)\right)+\bm{N},$ (12)

where the state-to-observable map is taken as the identity resulting in
$\mathcal{O}\equiv\mathcal{F}$. The OK model, its associated forward operator
$\mathcal{F}$, and sensor signal discretization $\bm{\mathcal{D}}$ are defined
in Section 7. Here we use the squared noise-covariance–weighted $l^{2}$ data
misfit131313The subscript $C_{\bm{N}}^{-1}$ stands for $l^{2}$ inner product
weighted by the matrix $C_{\bm{N}}^{-1}$. Here,
$C_{\bm{N}}=0.04\mathrm{Id}_{y}$ is a diagonal matrix, corresponding to a
signal-to-noise ratio of $5$.$\mathcal{M}_{f}$ defined as

$\mathcal{M}_{f}(\bm{x},z;\bm{y}^{*},C_{\bm{N}})\coloneqq\Big{\lVert}\bm{y}^{*}-\bm{\mathcal{D}}\big{(}\mathcal{O}(\bm{x},z)\big{)}\Big{\rVert}_{C_{\bm{N}}^{-1}}^{2}\,,$
(13)

to compare the synthetic data $\bm{y}^{*}$ that is generated from the fixed
latent inputs $(\bm{x}^{*},z^{*})$ via (12), with model-predicted patterns at
other inputs $(\bm{x},z)$. We note that the data misfit is proportional to the
negative log-likelihood generated by the sensor noise in (12).

In the first numerical example, the data misfit is evaluated at $\bm{x}^{*}$
and $10^{3}$ i.i.d. samples of the auxiliary variable,
$z^{i}\sim\nu_{Z|\bm{X}}(\cdot|\bm{x}^{*})$. In Fig. 3, a histogram of the
squared $l^{2}$ data misfit evaluations is shown along with the synthetic data
$\bm{y}^{*}$, the most likely observable corresponding to the lowest $l^{2}$
misfit, and the least likely observable corresponding to the highest $l^{2}$
misfit. In the second numerical example, the data misfit is evaluated at
$z^{*}$ and $10^{3}$ i.i.d. samples of the model parameter
$\bm{x}^{i}\sim\mathcal{N}(\bm{x}^{*},\delta^{2}\mathrm{Id}_{x})$, where
$\delta$ is chosen to be small141414We choose $\delta=0.05$ for the same
computational domain defined in B. The choice of $\delta$ is small comparing
to the admissible parameter region defined in Section 7.5.. Similar plots for
the second example are shown in Fig. 4.

The synthetic data $\bm{y}^{*}$
---

The most likely observable | The least likely observable
---|---
|

Figure 3: Observables $o^{i}=\mathcal{O}(\bm{x}^{*},z^{i})$, generated with
fixed model parameters $\bm{x}^{*}$ and random samples of the auxiliary
variable, $z^{i}\sim\nu_{Z|\bm{X}}(\cdot|\bm{x}^{*})$, have large misfit
values (_left_) with respect to a synthetic top-down microscopy
characterization image $\bm{y}^{*}$ (_right top_), generated by the data model
at $\bm{x}^{*}$ and $z^{*}$. However, the most likely and the least likely
observables (_right bottom_) are visually similar in terms of phase, spot
radius, and spot periodicity.

The synthetic data $\bm{y}^{*}$
---

The most likely observable | The least likely observable
---|---
|

Figure 4: Observables $o^{i}=\mathcal{O}(\bm{x}^{i},z^{*})$, generated with
small perturbations in the model parameters,
$\bm{x}^{i}\sim\mathcal{N}(\bm{x}^{*},\delta^{2}\mathrm{Id}_{x})$, and fixed
auxiliary variable $z^{*}$, have large misfit values (_left_) with respect to
a synthetic top-down microscopy characterization image $\bm{y}^{*}$ (_right
top_), generated by the data model at $\bm{x}^{*}$ and $z^{*}$. However, the
most likely and the least likely observables (_right bottom_) are visually
similar in terms of phase, spot radius, and spot periodicity. The expected
misfit for only the sensor noise $\bm{N}$ is included (_left, dashed line_) as
a reference to contrast with the large misfit values that arise from
perturbing the model parameters.

In both figures, we observe that, while the true parameters reproduce
observables identical to the image up to the additive noise, most realizations
of the nearby parameters or auxiliary variables produce observables that are
drastically different from the image data according to the $l^{2}$ data
misfit. As a result, the misfit values for most of the samples are very large,
relative to misfit values for the latent input that generated the synthetic
data. This leads to a likelihood concentration phenomenon that results in LFI
methods requiring many joint samples to recover the true parameters.

Furthermore, similar features in characterization images $\bm{y}^{*}$ of BCP
self-assembly indicate that certain characteristics of their latent BCP
equilibrium structures are similar. This should, in general, imply certain
closeness in their underlying model parameters. However, the likelihood
function $\pi_{\bm{y}^{*}|\bm{X}}(\bm{x})$ measures the probability of the
model parameterized by $\bm{x}$ matching a _specific_ image. As a result, the
observables that are visually similar to the image data, but are mismatched in
entry-wise values, will not be weighted heavily by the likelihood function, as
demonstrated through the two numerical examples above. This leads to
unnecessarily ineffective parameter inference.

Although using entire images for inference ensures that all of the available
information about the latent state is preserved, some information may be
irrelevant and/or redundant for the purpose of model calibration and
prediction. For example, the specific orientation of phase interfaces in a
top-down microscope image of unguided lamella phase BCP thin films (i.e. the
‘figer-print’ patterns) are irrelevant (when ignoring boundary effects) if our
interest is to predict the natural periodicity length of the BCP material. In
other words, compressing irrelevant information in the data should not affect
the posterior density for the parameters.

The notion of data reduction is very common for image data. While the $l^{2}$
distance can be easily computed, it does not capture the visual similarity
between images. For instance, the distance can be arbitrarily large when the
two images are identical, up to a rotation or reflection. Therefore, it is
often recommended to measure the discrepancy between images in a more natural
way. Different metrics that mitigate the issues for high-dimensional images
have been considered by many authors [84, 85, 86] and some image processing
techniques have been applied to image data for BCP self-assembly in the past
[39, 41, 26]. These metrics and the corresponding likelihood functions are
designed to be relatively insensitive to certain transformations that do _not_
change the fundamental characteristics of the image related to the model
parameters of interest.

In this work we advocate for constructing low-dimensional summary statistics
of BCP characterization images that quantitatively describe image features
that are related to characteristics of their latent BCP equilibrium
structures. Ideally, the summary statistics $\bm{Q}$ are chosen to be
sufficient statistics. This means that they do not introduce any bias in the
inference procedure, i.e., they satisfy
$\pi_{\bm{X}|\bm{Y}}(\bm{x}|\bm{y})=\pi_{\bm{X}|\mathbf{Q}}(\bm{x}|\mathcal{Q}(\bm{y}))$
for all $\bm{x},\bm{y}$. Furthermore, we can also look for minimal sufficient
statistics, i.e., the summary statistics that most reduce the dimension of the
data without discarding any relevant information. Unfortunately, sufficient
statistics are not readily available outside of exponential family models, and
consequently many methods have been proposed to select summary (but not
necessarily sufficient) statistics for LFI in practice [87]. These summary
statistics may discard some information during the inference, as a trade-off
for suffering less from the concentration phenomenon of the likelihood for
high-dimensional images [88].

It is important to note that working with summary statistics of the data
remains feasible with LFI methods. While the conditional distributions for the
summary statistics $\pi_{\bm{Q}|\bm{X}}(\cdot|\bm{x})$ may not be analytically
available or feasible to evaluate, LFI methods only require joint samples of
the parameters $\bm{x}^{i}$ and summary statistics $\bm{q}^{i}$. These samples
can be acquired by an additional step of computing the summary statistics
$\bm{q}^{i}=\mathcal{Q}(\bm{y}^{i})$ for each data sample $\bm{y}^{i}$.

##  5 Likelihood-free inference via triangular transport maps

In this section, we present a measure transport approach for characterizing
the family of conditional densities for the parameters given the observations
$\pi_{\bm{X}|\bm{Q}}$ using only samples from the joint distribution
$\nu_{\bm{X},\bm{Q}}$. In Section 5.1 and 5.2 we introduce the notion of
measure transport and triangular transport maps, respectively. In Section 5.3
we show how to estimate the maps from samples. Then, in Section 5.4, we
demonstrate how these maps can be used for likelihood–free inference. A
numerical implementation of this approach is available online in a MATLAB
package for transport map estimation151515https://github.com/baptistar/ATM.

###  5.1 Measure transport

Let $\bm{W}\in\mathbb{R}^{w}$ be a _target_ random variable with measure
$\nu_{\bm{W}}\in\mathcal{P}(\mathbb{R}^{w})$ we would like to characterize,
and let $\bm{R}\in\mathbb{R}^{w}$ be a _reference_ random variable with
measure $\nu_{\bm{R}}\in\mathcal{P}(\mathbb{R}^{w})$ that we can easily sample
from, for example, a standard Gaussian. Our goal is to find a measurable
function $S\colon\mathbb{R}^{w}\rightarrow\mathbb{R}^{w}$ that defines a
deterministic coupling between these two variables such that the law of
$S(\bm{W})$ is equal to the law of $\bm{R}$. For such a map $S$, we say that
$S$ _pushes forward_ $\nu_{\bm{W}}$ to $\nu_{\bm{R}}$ and we denote this by
$S_{\sharp}\nu_{\bm{W}}=\nu_{\bm{R}}$. Equivalently, we say that $S$ _pulls
back_ $\nu_{\bm{R}}$ to $\nu_{\bm{W}}$ and we denote this by
$S^{\sharp}\nu_{\bm{R}}=\nu_{\bm{W}}$. For measures $\nu_{\bm{W}}$ and
$\nu_{\bm{R}}$ with densities $\pi_{\bm{W}}$ and $\pi_{\bm{R}}$, respectively,
similar relationships apply, $S_{\sharp}\pi_{\bm{W}}=\pi_{\bm{R}}$ and
$S^{\sharp}\pi_{\bm{R}}=\pi_{\bm{W}}$. If $S$ is a diffeomorphism (i.e., a
differentiable bijection with a differentiable inverse), we can evaluate the
pullback density $S^{\sharp}\nu_{\bm{R}}$ using the change-of-variables
formula as

$S^{\sharp}\pi_{\bm{R}}(\bm{w})\coloneqq\pi_{\bm{R}}\circ S(\bm{w})|\det\nabla
S(\bm{w})|.$ (14)

Therefore, having the transport map $S$ that couples $\bm{W}$ and $\bm{R}$ and
the reference density $\pi_{\bm{R}}$ enables us to evaluate the target density
$\pi_{\bm{W}}$. Furthermore, we can use this map to generate i.i.d. samples
from the target distribution by sampling from the reference measure
$\bm{r}^{i}\sim\nu_{\bm{R}}$ and applying the inverse map to generate samples
$S^{-1}(\bm{r}^{i})\sim\nu_{\bm{W}}$. In Fig. 5, we provide an example of
using $S$ to characterize the pullback density and generate samples from a
2-dimensional “banana” target distribution.

$S^{\sharp}\pi_{\mathbf{R}}$$S_{\sharp}\pi_{\mathbf{W}}$$S^{-1}(\bm{r}^{i})$$S(\bm{w}^{i})$

Figure 5: The target density $\pi_{\bm{W}}$ is represented as the pullback of
a standard Gaussian reference density $\pi_{\bm{R}}$ through the transport map
$S$ (_left_). Samples from $\nu_{\bm{R}}$ are generated by mapping reference
samples through the inverse of $S$ (_right_). Similarly we can represent the
reference density as the pushforward of the target density through $S$ and
sample from $\nu_{\bm{W}}$ by mapping target samples through $S$.

The existence of a map $S$ that pushes forward $\pi_{\bm{W}}$ to
$\pi_{\bm{R}}$ is guaranteed under very mild conditions such as $\nu_{\bm{W}}$
and $\nu_{\bm{R}}$ being absolutely continuous with respect to the Lebesgue
measure, i.e., the distributions do not contain atoms and their densities
exist. We refer the reader to [89] for approaches that relax these conditions.

###  5.2 The Knothe–Rosenblatt rearrangement

In general, many maps $S$ can satisfy the measure transport conditions above
[90]. For instance, the optimal transport literature looks for maps $S$ that
couple two distributions and also minimize an integrated transportation
distance, such as the Wasserstein metric [91]. In this work, we consider the
Knothe–Rosenblatt (KR) rearrangement for $S$, i.e., the unique monotone and
lower triangular map that pushes forward $\pi_{\bm{W}}$ to $\pi_{\bm{R}}$ [92,
93, 94]. A monotone and lower triangular map is a multivariate function of the
form

$S(\bm{w})=\begin{bmatrix}[l]S^{1}(w_{1})\\\ S^{2}(w_{1},w_{2})\\\ \vdots\\\
S^{w}(w_{1},\dots,w_{w})\end{bmatrix},$ (15)

where $\xi\rightarrow
S^{k}(w_{1},\dots,w_{k-1},\xi)\colon\mathbb{R}\to\mathbb{R}$, i.e., the
restriction of the map component $S^{k}$ to its last variable $w_{k}$, by
fixing its first $k-1$ variables, is monotone increasing for all
$(w_{1},\dots,w_{k-1})\in\mathbb{R}^{k-1}$. We denote the space of triangular
and monotone maps on $\mathbb{R}^{w}$ by
$\mathcal{S}_{\triangle}(\mathbb{R}^{w})$.

Working with the KR rearrangement, rather than other (non-triangular)
transport maps, has several computational benefits [38]. First, monotone
triangular maps $S\in\mathcal{S}_{\triangle}(\mathbb{R}^{w})$ can be easily
inverted to sample from $\nu_{\bm{W}}$ by solving the nonlinear system
$S(\bm{w}^{i})=\bm{r}^{i},\;\bm{r}^{i}\sim\nu_{\bm{R}}$ via a sequence of
1-dimensional nonlinear solves. Furthermore, each nonlinear solve has a unique
solution due to the injectivity of the map. Second, the Jacobian determinant
of a triangular map $\det\nabla S(\bm{w})$ is tractable to compute as the
product of the partial derivatives of each component with respect to the last
variables, i.e., $\det\nabla
S(\bm{w})=\prod_{k=1}^{w}\partial_{k}S^{k}(\bm{w}_{1:k})$. Third, these maps
are easier to parameterize than maps derived via optimal transport and can be
computed quickly using gradient-based optimization techniques [95, 96].

Most importantly in our context, these maps provide direct access to the
conditional densities of the target distribution. For a reference distribution
with independent components, such as the standard Gaussian, the $k$ map
component $S^{k}$ pushes forward the $k$th marginal conditional density
$\pi_{\bm{W}_{k}|\bm{W}_{1},\dots,\bm{W}_{k-1}}$ in the factorization of the
target density
$\pi_{\bm{W}}=\pi_{\bm{W}_{1}}\pi_{\bm{W}_{2}|\bm{W}_{1}}\dots\pi_{\bm{W}_{w}|\bm{W}_{1},\dots,\bm{W}_{w-1}}$,
to the $k$th marginal of the reference density $\pi_{\bm{R}}$ [97, App. A]. In
Section 5.4, we use this property to extract the posterior density, i.e., the
conditional density of the parameters given the data, from their joint
density.

###  5.3 Estimating the Knothe–Rosenblatt rearrangement

The _exact_ KR rearrangement $S_{\textrm{KR}}$ that couples two densities is,
in general, analytically unavailable. However, we consider approximations of
the KR map in this work. We follow the variational approach in [38, 97, 96] of
minimizing the KL divergence,161616While we choose the KL divergence, one can
also consider other divergences or metrics to measure the discrepancy between
densities. from the pullback density $S^{\sharp}\pi_{\bm{R}}$ to the target
density $\pi_{\bm{W}}$, which has the form

$\mathcal{D}_{\textrm{KL}}\left(\pi_{\bm{W}}||S^{\sharp}\pi_{\bm{R}}\right)=\mathbb{E}_{\nu_{\bm{W}}}\left[\log\pi_{\bm{W}}-\log\pi_{\bm{R}}\circ
S(\bm{W})-\log\det\nabla S(\bm{W})\right].$

The KR rearrangement is given by the global minimizer of the KL divergence
over the space of lower triangular and monotone maps on $\mathbb{R}^{w}$. That
is,

$\displaystyle S_{\textrm{KR}}$
$\displaystyle=\underset{S\in\mathcal{S}_{\triangle}(\mathbb{R}^{w})}{\operatorname*{arg\,min}}\;\mathcal{D}_{\textrm{KL}}(\pi_{\bm{W}}||S^{\sharp}\pi_{\bm{R}})$
$\displaystyle=\underset{S\in\mathcal{S}_{\triangle}(\mathbb{R}^{w})}{\operatorname*{arg\,min}}\;\mathbb{E}_{\nu_{\bm{W}}}\left[-\log\pi_{\bm{R}}\circ
S(\bm{W})-\log\det\nabla S(\bm{W})\right].$ (16)

In all of our numerical experiments, we let the reference distribution be
standard Gaussian. Substituting the form of the reference density
$\pi_{\bm{R}}=\mathcal{N}(\mathbf{0},\mathrm{Id}_{w})$ in (5.3), the
optimization problem can be written as

$S_{\textrm{KR}}=\underset{S\in\mathcal{S}_{\triangle}(\mathbb{R}^{w})}{\operatorname*{arg\,min}}\;\sum_{k=1}^{w}\mathbb{E}_{\nu_{\bm{W}}}\left[\frac{1}{2}S^{k}(\bm{W}_{1:k})^{2}-\log\partial_{k}S^{k}(\bm{W}_{1:k})\right],$
(17)

where each term in the sum only depends on one map component. Thus, optimizing
for $S$ can be decoupled into $w$ optimization problems for each map
component171717In fact, this holds for any reference density with independent
components., which can each be solved independently and in parallel to
determine each $S^{k}$. More information on the parametrization and
computation of monotone and triangular maps is provided in C.

###  5.4 Conditional density estimation via triangular transport maps

To estimate the conditional densities $\pi_{\bm{X}|\bm{Q}}$ and any posterior
$\pi_{\bm{X}|\bm{q}^{*}}$ in particular, we begin by considering the KR
rearrangement $S_{\textrm{KR}}$ that pushes forward the target joint density
$\pi_{\bm{X},\bm{Q}}$ to a reference density with independent components
$\pi_{\bm{R}_{1},\bm{R}_{2}}=\pi_{\bm{R}_{1}}\pi_{\bm{R}_{2}}$, that are both
defined on $\mathbb{R}^{x+q}$. Given that $S_{\textrm{KR}}$ is lower
triangular, we can partition the map into two blocks as

$S_{\textrm{KR}}(\bm{q},\bm{x})=\begin{bmatrix}[l]S^{\mathcal{Q}}(\bm{q})\\\
S^{\mathcal{X}}(\bm{q},\bm{x})\end{bmatrix},$ (18)

where $S^{\mathcal{Q}}\colon\mathbb{R}^{q}\rightarrow\mathbb{R}^{q}$ and
$S^{\mathcal{X}}(\bm{q},\cdot)\colon\mathbb{R}^{x}\rightarrow\mathbb{R}^{x}$
are monotone functions for all $\bm{q}\in\mathbb{R}^{q}$.From the measure
transport condition, then $S^{\mathcal{Q}}$ pulls back the marginal
$\pi_{\bm{R}_{1}}$ to the marginal for the summary statistics $\pi_{\bm{Q}}$,
and $S^{\mathcal{X}}(\bm{q},\cdot)$ pulls back the marginal $\pi_{\bm{R}_{2}}$
to the conditional $\pi_{\bm{X}|\bm{Q}}(\cdot|\bm{q})$ for each
$\bm{q}\in\mathbb{R}^{q}$.

From the decoupling of the optimization problem over the map components, we
only need to solve $x$ independent optimization problems for the components in
$S^{\mathcal{X}}$ to approximate the conditional density
$\pi_{\bm{X}|\bm{Q}}$. Moreover, the conditional density can be estimated
given $N$ i.i.d. samples from the joint distribution,
$\\{(\bm{x}^{i},\bm{q}^{i})\\}_{i=1}^{N}\sim\nu_{\bm{X},\bm{Q}}$. We use these
samples to define a Monte Carlo approximation of the expectation in (17) and
to find the $k$th component of $S^{\mathcal{X}}$ by solving

$(\widehat{S}^{\mathcal{X}})^{q+k}=\underset{S^{q+k}}{\operatorname*{arg\,min}}\;\frac{1}{N}\sum_{i=1}^{N}\left[\frac{1}{2}S^{q+k}(\bm{q}^{i},\bm{x}^{i}_{1:k})^{2}-\log\partial_{k}S^{q+k}(\bm{q}^{i},\bm{x}^{i}_{1:k})\right],$
(19)

under the monotonicity constraint
$\partial_{k}(S^{\mathcal{X}})^{k}(\bm{q},\bm{x}_{1:k})>0$ for all
$(\bm{q},\bm{x}_{1:k})\in\mathbb{R}^{q+k}$. Then, collecting the components
for $k=1,\dots,x$ in $\widehat{S}^{\mathcal{X}}$ produces an estimator of the
conditional density $\pi_{\bm{X}|\bm{Q}}(\cdot|\bm{q})$ for any
$\bm{q}\in\mathbb{R}^{q}$ as

$\displaystyle\bm{x}\mapsto\widehat{\pi}_{\bm{X}|\bm{Q}}(\bm{x}|\bm{q})$
$\displaystyle\coloneqq\widehat{S}^{\mathcal{X}}(\bm{q},\bm{x})^{\sharp}\pi_{\bm{R}_{2}}(\bm{x})$
$\displaystyle=\pi_{\bm{R}_{2}}\left(\widehat{S}^{\mathcal{X}}(\bm{q},\bm{x})\right)\det\nabla_{\bm{x}}\widehat{S}^{\mathcal{X}}(\bm{q},\bm{x}).$
(20)

An example of the estimated conditional densities
$\widehat{\pi}_{\bm{X}|\bm{Q}}(\cdot|\bm{q})$ extracted from the transport map
$\widehat{S}^{\mathcal{X}}(\bm{q},\cdot)$ is shown in Fig. 6. We note that a
single map $S^{\mathcal{X}}$ that is a function of the summary statistics
$\bm{q}$ can characterize conditional densities with different properties
across the data space, such as unimodal and bimodal behavior.

Figure 6: Conditional density estimation for the posterior densities with the
quadratic likelihood model: $Y^{2}=X^{2}+N$ for $X\sim\mathcal{N}(0,0.25)$ and
$N\sim\mathcal{N}(0,0.01)$. Contours of the joint density of the parameter $x$
and the observable $y$ are on the left, and the posterior densities for three
values of the observable $y^{*}\in\\{0,0.25,0.5\\}$ and the monotone maps that
pull back a standard normal density to each posterior are plotted on the
right.

In addition to having a functional form for the conditional density
$\widehat{\pi}_{\bm{X}|\bm{Q}}(\cdot|\bm{q})$, we can also use the estimated
map to generate i.i.d. samples $\bm{x}^{i}$ from the conditional distribution
$\widehat{\nu}_{\bm{X}|\bm{Q}}(\cdot|\bm{q})$ by sampling $\bm{r}^{i}$ from
the reference measure and solving for the parameter $\bm{x}^{i}$ that is
mapped to $\bm{r}^{i}$ by the function
$\widehat{S}^{\mathcal{X}}(\bm{q},\cdot)$, i.e.,

$\begin{rcases}\bm{r}^{i}\sim\nu_{\bm{R}_{2}}\\\
\widehat{S}^{\mathcal{X}}(\bm{q},\bm{x}^{i})=\bm{r}^{i}\end{rcases}\implies\bm{x}^{i}\sim\widehat{\nu}_{\bm{X}|\bm{Q}}(\cdot|\bm{q}).$
(21)

To summarize the approach outlined in this section, we use a sample
representation $(\bm{x}^{i},\bm{q}^{i})$ of the joint distribution
$\nu_{\bm{X},\bm{Q}}$ to compute an approximate map
$\widehat{S}^{\mathcal{X}}$ that couples a reference density to the
conditional densities $\pi_{\bm{X}|\bm{Q}}$; this map enables both evaluation
of approximate posterior densities $\widehat{\pi}_{\bm{X}|\bm{q}^{*}}$ for any
and all summary statistics $\bm{q}^{*}$, via (20), and simulation of i.i.d.
samples from their respective distributions
$\widehat{\nu}_{\bm{X}|\bm{Q}}(\cdot|\bm{q}^{*})$, via (21).

##  6 Estimating expected information gain via triangular transport maps

A core advantage of estimating the conditional densities for _all values_ of
$\bm{q}$ is that it enables the computation of expectations over the data
random variables. One of these quantities is the EIG defined in Section 3.3.
An example of using EIG to compare the informativeness of two summary
statistics for parameter inference is shown in Fig. 7. In Section 6.1 we
describe existing estimators for EIG and in Section 6.2 we present a lower
bound for the EIG based on the transport maps computed in Section 5.

Figure 7: Posterior densities for the Gaussian $X$ parameter in Fig. 6 given
two summary statistics $Q_{1}=-X+N$ (_left_) and $Q_{2}=X^{2}+N$ (_right_) for
$N\sim\mathcal{N}(0,0.01)$. The posteriors for $Q_{1}$ are concentrated as
compared to the bi-modal posteriors for $Q_{2}$, indicating that $Q_{1}$ is
more informative for predicting $X$. Integrating the KL divergence between the
prior for $X$ and the posteriors over the different values of $Q_{i}$ results
in EIGs of $\widehat{I}(X;Q_{1})=1.64\pm 0.01$ and
$\widehat{I}(X;Q_{2})=0.86\pm 0.01$. This provides a quantitative approach to
measure and rank the informativeness of different summary statistics.

###  6.1 Related work

Given conditional densities $\pi_{\bm{Q}|\bm{X}}$ that are tractable to
evaluate, several nested Monte Carlo (NMC) estimators have been proposed to
estimate the EIG using samples from the joint distribution of parameters and
observables [98, 99, 100]. The EIG in (10a) and (10b) contains an expectation
over $\nu_{\bm{X},\bm{Q}}$ and an unknown marginal density for the summary
statistics that is typically unavailable in closed form. To estimate the EIG,
NMC uses $N_{1}$ i.i.d samples
$\\{(\bm{x}^{i},\bm{q}^{i})\\}_{i=1}^{N_{1}}\sim\nu_{\bm{X},\bm{Q}}$ to
estimate the outer expectation as
$I(\bm{X};\bm{Q})\approx\frac{1}{N_{1}}\sum_{i=1}^{N_{1}}\log\pi_{\bm{Q}|\bm{X}}(\bm{q}^{i}|\bm{x}^{i})-\log\pi_{\bm{Q}}(\bm{q}^{i})$
and an independent set of $N_{2}$ i.i.d samples
$\\{\bm{x}^{j}\\}_{j=1}^{N_{2}}\sim\nu_{\bm{X}}$ to approximate the integral
for the marginal density of the summary statistics at $\bm{q}^{i}$ as
$\pi_{\bm{Q}}(\bm{q}^{i})\approx\frac{1}{N_{2}}\sum_{j=1}^{N_{2}}\pi_{\bm{Q}|\bm{X}}(\bm{q}^{i}|\bm{x}^{j})$.
Let us note that estimates for the marginal density are often numerically zero
from arithmetic underflow when evaluating the likelihood function in high
dimensions.

To address this, Beck et. al [101] proposed a Laplace approximation–based
importance sampling scheme to estimate the marginal likelihood. Although NMC
is asymptotically unbiased for estimating EIG as
$(N_{1},N_{2})\rightarrow\infty$, NMC is often considered too computationally
costly. Even with an optimal allocation of samples to the inner and outer
estimators, the mean squared error of a NMC estimator for EIG converges at a
rate of $\mathcal{O}(N^{-2/3})$, where $N$ is the total number of samples.
However, in order to achieve this rate, one needs to have
$N_{2}=c\sqrt{N_{1}}$, and finding what $c$ is involves some tuning [102].

To alleviate this cost, Wu et. al [103] considered EIG estimation based on
Laplace approximations that exploit likelihood-informed subspaces and data
dimension reduction. The authors numerically demonstrate that various biased
estimators are often sufficient to perform BOED with unimodal posteriors that
arise from affine or mildly nonlinear forward operators, unimodal priors, and
additive Gaussian observation errors with small variance. The downside,
however, is that accurate approximations or bounds on EIG are no longer
guaranteed, and (gradient) information of the forward operator is required.

In a likelihood-free setting, estimating EIG requires us to characterize two
densities, the prior and posterior, or the likelihood and marginal density of
the data, which can be computationally expensive. One of the most well-known
methods proposed in [104] is to estimate mutual information (MI) using
k-nearest neighbor entropy estimators. Another method is to estimate MI using
kernel density estimators (KDE) [105]. While the former requires the number of
samples to scale exponentially with the MI, the latter KDE often performs
poorly in high dimensions. To circumvent these challenges, Suzuki et. al [106]
proposed to approximate EIG by directly estimating a ratios of densities.
Thomas et. al [82] proposed an approach using logistic regression to classify
prior and posterior samples in order to estimate the ratio of the prior and
posterior densities. Although this method allows us to estimate the density
ratio directly, it requires solving an optimization problem for each
observation. Alternatively, several approaches have been proposed in the
machine learning literature for constructing variational bounds for MI as
proposed in [44, 107]. In [44], the approximate posterior is obtained by
optimizing over a parameterized family of distributions, e.g., the mean-field
approximation. The choice can introduce bias in the EIG estimation.

###  6.2 Likelihood-free estimators for expectation information gain

We now describe an approach to estimate the EIG that uses the variational
approximation $\widehat{\pi}_{\bm{X}|\bm{Q}}$ in (20) for the posterior
density $\pi_{\bm{X}|\bm{Q}}$. A simple calculation [44, Appx. A] shows that a
lower bound for the EIG in (10a) is:

$\displaystyle
I(\bm{X};\bm{Q})\geq\widehat{I}(\bm{X};\bm{Q})\coloneqq\mathbb{E}_{\nu_{\bm{X},\bm{Q}}}\left[\log\frac{\widehat{\pi}_{\bm{X}|\bm{Q}}}{\pi_{\bm{X}}}\right].$
(22)

Moreover, the gap between the EIG and the lower bound is equal to the
expectation of the KL divergence from the variational approximation to the
true posterior, i.e.,
$\mathbb{E}_{\nu_{\bm{Q}}}[\mathcal{D}_{\textrm{KL}}(\pi_{\bm{X}|\bm{Q}}(\cdot|\bm{q})||\widehat{\pi}_{\bm{X}|\bm{Q}}(\cdot|\bm{q}))]$.
Thus, the lower bound for the EIG provided by the variational approximation
approaches the true EIG as $\widehat{\pi}_{\bm{X}|\bm{Q}}(\cdot|\bm{q})$
better approximates the true posterior density.

Given $N$ i.i.d. samples from the joint density of parameters and summary
statistics, an unbiased estimate for the lower bound in (22) is

$\mathbb{E}_{\nu_{\bm{X},\bm{Q}}}\left[\log\frac{\widehat{\pi}_{\bm{X}|\bm{Q}}}{\pi_{\bm{X}}}\right]\approx\frac{1}{N}\sum_{i=1}^{N}\log\widehat{\pi}_{\bm{X}|\bm{Q}}(\bm{x}^{i}|\bm{q}^{i})-\mathbb{E}_{\nu_{\bm{X}}}[\log\pi_{\bm{X}}],$
(23)

where $\\{(\bm{x}^{i},\bm{q}^{i})\\}_{i=1}^{N}\sim\nu_{\bm{X},\bm{Q}}$. The
second term in (23) is the entropy of the prior density for the parameters.
When the prior density is specified in a parametric family, this quantity is
often available in closed form. As an example, for a standard Gaussian prior
density $\pi_{\bm{X}}=\mathcal{N}(0,\mathrm{Id}_{x})$ where
$\bm{X}\in\mathbb{R}^{x}$, the entropy is given by
$\mathbb{E}_{\nu_{\bm{X}}}[-\log\pi_{\bm{X}}]=\frac{1}{2}\log(2\pi ex)$.

We note that in the general likelihood–free setting when the prior density
$\pi_{\bm{X}}$ is also analytically unavailable, we can also find a
variational approximation $\widehat{\pi}_{\bm{X}}$ for $\pi_{\bm{X}}$ and
define an estimator for the EIG as

$I(\bm{X};\bm{Q})\approx\mathbb{E}_{\nu_{\bm{X},\bm{Q}}}\left[\log\frac{\widehat{\pi}_{\bm{X}|\bm{Q}}}{\widehat{\pi}_{\bm{X}}}\right]\approx\frac{1}{N}\sum_{i=1}^{N}\left[\log\widehat{\pi}_{\bm{X}|\bm{Q}}(\bm{x}^{i}|\bm{q}^{i})-\log\widehat{\pi}_{\bm{X}}(\bm{x}^{i})\right],$
(24)

where $\\{(\bm{x}^{i},\bm{q}^{i})\\}_{i=1}^{N}\sim\nu_{\bm{X},\bm{Q}}$.
However, as shown in [44], the estimator in (24) does not yield a bound on
EIG. This is problematic if the goal is to maximize the EIG with respect to
different design variables in the context of optimal experimental design.

An important property of MI, and hence EIG, is its invariance under marginal
transformations of the random variables [108]. If we have two bijective
functions $f$ and $g$, then $I(\bm{X};\bm{Q})=I(f(\bm{X});g(\bm{Q}))$. The
property allows us to compute the EIG by finding variational approximations on
transformed random variables so that the resulting EIG estimator has better
statistical properties. In our numerical experiments, we transform the model
parameters $\bm{X}$ defined on a bounded subset of $\mathbb{R}^{3}_{+}$ to the
unbounded domain $\mathbb{R}^{3}$ using the map described in Section 7.5. The
transformed random variables have a standard Gaussian prior distribution
allowing us to use the EIG estimator in (23). Moreover, the distributions for
the transformed parameters satisfy the absolute continuity assumptions
required by our approach to compute the transport maps in Section 5. A similar
transformation can be defined for the summary statistics. For instance, the
function $g(\bm{Q})=\log(\bm{Q})$ transforms variables fully supported on the
positive numbers $\mathbb{R}_{+}$ to be fully supported on $\mathbb{R}$.

##  7 Case study: Model definitions

In the next two sections, we perform Bayesian calibration of the OK model for
the unguided self-assembly of Di-BCP thin films with top-down microscopy image
data. In this section, we focus on defining the data model and the prior
distribution for the model parameters used by the LFI and EIG computation
methodologies described above. In the next section, we present numerical
results of model calibration and EIG computation.

We first introduce the OK model for unguided Di-BCP thin film self-assembly
and the model parameters in Section 7.1. The associated forward operator and
auxiliary variable are defined in Section 7.2. We then define a model for top-
down microscopy characterization that contains restriction operators that
mimic microscope magnification, a blurring convolution, and additive sensor
noise in Section 7.3. Several low-dimensional summary statistics that
quantitatively describe relevant features in the image data are proposed in
Section 4.4, which are categorized into energy-based statistics and Fourier-
based statistics. Lastly, we construct a model-specific prior for the
parameters based on assumptions of experimental conditions in Section 7.5.

###  7.1 The Ohta–Kawasaki model

Consider a Di-BCP melt, composed of two types of monomers that are labeled A
and B respectively, in a physical domain $\mathcal{D}\subset\mathbb{R}^{3}$.
Let $u_{A}$ and $u_{B}\colon\mathcal{D}\to[0,1]$, denote the normalized
segment densities for the two monomers. These segment densities satisfy the
incompressibility constraint $u_{A}+u_{B}=1$. Then, our state variable
$u\coloneqq u_{A}-u_{B}\in[-1,1]$, known as the order parameter, represents
the difference between the normalized local fractions of the two phases. Under
the thin film approximation, we may only consider the variation of the order
parameter parallel to the substrate; thus $\mathcal{D}\subset\mathbb{R}^{2}$
is assumed in our case study181818If one is interested in model calibration
with respect to the effects of film thickness and/or polymer-substrate
interaction, then the combination of X-ray scattering, which has the ability
to probe the internal structure of BCP melts [40, 109], and the OK model in
the full 3D domain or $x_{1}$-$x_{3}$ cross-section domain should be
considered. We refer to works discussed in Section 1.4 for setups of such a
calibration task..

The free energy of a Di-BCP melt, $\mathcal{E}^{\textrm{OK}}$, given by the OK
model [34, 21] is assumed to be the sum of three contributing energies: the
_double well energy_ $\mathcal{D}^{\textrm{OK}}$, the _interfacial energy_
$\mathcal{I}^{\textrm{OK}}$, and the _non-local energy_
$\mathcal{N}^{\textrm{OK}}$, scaled by the model parameters
$(\kappa,\epsilon,\sigma)\in\mathbb{R}_{+}^{3}$, respectively. That is,

$\displaystyle\mathcal{E}^{\textrm{OK}}(u;\kappa,m,\epsilon,\sigma)\coloneqq\kappa\mathcal{D}^{\textrm{OK}}(u)+\epsilon^{2}\mathcal{I}^{\textrm{OK}}(u)+\sigma\mathcal{N}^{\textrm{OK}}(u;m).$
(25)

The double well energy $\mathcal{D}^{\textrm{OK}}$ quantifies the free energy
associated with the mixing of monomers and favors pure monomer phases
represented by values $\pm 1$. The interfacial energy
$\mathcal{I}^{\textrm{OK}}$ is a $H^{1}(\mathcal{D})$ seminorm that
characterizes the short-range interaction between monomers and controls the
interfacial length between different monomer phases. Lastly, the nonlocal
energy $\mathcal{N}^{\textrm{OK}}$ characterizes the long-range interaction
between monomers and controls the expansion of phase boundaries. The
competition between the three energies during free-energy minimization forms a
length scale selection mechanism that leads to various structures of
equilibrium Di-BCP melts depicted by the experimental phase diagram in Fig. 1
(_left_). It has been shown numerically that most periodic structures observed
in experiments are minimizers of the OK free energy [110, 111].

The form of these contributing energies is given by:

$\begin{gathered}\mathcal{D}^{\textrm{OK}}(u)=\int_{\mathcal{D}}\mathcal{W}(u(\bm{s}))\textrm{d}\bm{s},\quad\mathcal{I}^{\textrm{OK}}(u)=\frac{1}{2}\int_{\mathcal{D}}|\nabla
u(\bm{s})|^{2}\textrm{d}\bm{s},\\\
\mathcal{N}^{\textrm{OK}}(u;m)=\frac{1}{2}\int_{\mathcal{D}}(u(\bm{s})-m)(-\triangle_{N})^{-1}(u(\bm{s})-m)\textrm{d}\bm{s}\;,\end{gathered}$
(26)

where $m\in[-1,1]$ is the spatial mean $\mathcal{M}$ of the order parameter,
i.e.,

$m=\mathcal{M}(u)=|\mathcal{D}|^{-1}\int_{\mathcal{D}}u(\bm{s})\textrm{d}\bm{s},$
(27)

and $\mathcal{W}(u)=(1-u^{2})^{2}/4$ is the double well potential that is
plotted in Fig. 8. The action of the inverse Laplacian is defined by the
solution of a Poisson equation with homogeneous Neumann boundary
conditions191919The periodic boundary condition is often used for square or
cubic domains, mainly for computational convenience and the minimization of
edge effects [19, p. 96]. The Neumann boundary additionally models the effect
of a bounded (confined) domain by assuming that polymer chains reflect at the
boundary [112].:

$w=(-\triangle_{N})^{-1}(u-m)\iff\begin{cases}-\triangle
w(\bm{s})=u(\bm{s})-m&\text{in }\mathcal{D}\;;\\\ \nabla
w(\bm{s})\cdot\bm{n}=0&\text{on }\partial\mathcal{D}\;,\end{cases}$

where $\bm{n}$ is the unit vector normal to $\partial\mathcal{D}$.

Figure 8: The quartic double well potential, which has $-1$ and $1$ as its
minima. The dashed lines delineate physical boundaries. Since $-1$ and $1$
refer to normalized pure densities of $u_{A}$ and $u_{B}$, respectively, when
the order parameter is outside of $(-1,1)$, the model does not represent a
physically meaningful quantity. In other words, this double well potential
$\mathcal{W}$ is _non-degenerate_. A _degenerate_ , or _obstacle_ , double
well potential, such as described in [113], produces an infinite wall at these
boundaries.

The parameters in the OK model are connected to the commonly used material
parameters of diblock copolymers. In particular, an appeal to the higher-
fidelity SCFT model in [35] reveals the relationships:

$\kappa=1\;,\quad\epsilon^{2}=\frac{l^{2}}{3f(1-f)\chi}\;,\quad\sigma=\frac{36}{f^{2}(1-f)^{2}l^{2}\chi
N^{2}}\;,\quad m=2f-1,$ (28)

where $N$ is the degree of polymerization of the polymer, $\chi$ is the Flory-
Huggins parameter describing the strength of repulsion between the two
monomers, $l$ is the statistical segment length of the polymer202020The
statistical segment length is assumed to be the same for the two types of
monomers in a Di-BCP, thus we refer to the OK model as a model of
conformationally symmetric Di-BCP melts., and $f\in(0,1)$ represents the
volume fraction of monomer $A$. In this study, we take the common assumption
of $\kappa=1$ and calibrate the collection of model parameters
$\bm{x}=(m,\epsilon,\sigma)$. We refer to the three parameters as the _spatial
average parameter_ , the _interfacial parameter_ , and the _nonlocal
parameter_ , respectively.

###  7.2 The forward operator of the Ohta–Kawasaki model

The forward operator $\mathcal{F}$ is defined as the map from parameters
$\bm{x}$ and a random initial state $z$ to a local minimizer of the OK energy,
$\mathcal{E}^{\textrm{OK}}_{\bm{x}}(u)\coloneqq\mathcal{E}^{\textrm{OK}}(u;1,m,\epsilon,\sigma)$,
for the given parameters. That is,

$u=\mathcal{F}(\bm{x},z)\in\operatorname*{arg\,min}_{u\in
V^{u}}\mathcal{E}^{\textrm{OK}}_{\bm{x}}(u)$ (29)

where the state space is defined as

$V^{u}\coloneqq\left\\{v\in
H^{1}(\mathcal{D}):\int_{\mathcal{D}}\left(v(\bm{s})-m\right)\textrm{d}\bm{s}=0\right\\}.$
(30)

As a result of the non-convexity of the double well potential $\mathcal{W}$,
the OK energy is non-convex and does not have a unique minimizer. Thus, the
local minimizers in (29) consist of an ensemble of numerous achievable
equilibrium patterns for Di-BCP thin film self-assembly. We note that the
states satisfying the first order condition for the minimization problem in
(29) can be found by solving the following mixed problem [114].

$\displaystyle\text{Find }(u,\mu)\in\left(H^{1}(\mathcal{D})\right)^{2}\text{
such that:
}\begin{cases}(\nabla\mu,\nabla\tilde{u})_{L^{2}}+\left(\sigma(u-m),\tilde{u}\right)_{L^{2}}=0,&\forall\tilde{u}\in
H^{1}(\mathcal{D})\;;\\\
(\mu,\tilde{\mu})_{L^{2}}-\left(\mathcal{W}^{\prime}(u),\tilde{\mu}\right)_{L^{2}}-(\epsilon^{2}\nabla
u,\nabla\tilde{\mu})_{L^{2}}=0,&\forall\tilde{\mu}\in
H^{1}(\mathcal{D})\;.\end{cases}$ (31)

where $(\cdot,\cdot)_{L^{2}}$ denotes $L^{2}(\mathcal{D})$ or
$(L^{2}(\mathcal{D}))^{2}$ inner product.

The evaluation of the forward operator $\mathcal{F}(\bm{x},z)$ is accomplished
by minimizing the energy $\mathcal{E}^{\textrm{OK}}_{\bm{x}}$ via a
variational procedure, starting from an initial state $z$. Conventionally, an
$H^{-1}$ gradient flow approach is used, which involves mimicking a time-
dependent process by solving the nonlocal Cahn–Hilliard equation [110, 111,
23], a fourth-order time-evolving partial differential equation (PDE), whose
steady-state solutions satisfy (31). According to the dynamical DFT [115], the
gradient flow approach may faithfully represent the kinetic pathway of the
self-assembly if an appropriate mobility operator [116] and/or a stochastic
noise that satisfies the fluctuation–dissipation theorem [20] are/is
incorporated. Instead, the OK model is more commonly treated as an equilibrium
model, i.e., the thermal fluctuations and dynamics of Di-BCPs are not modeled.

Several authors have considered a diverse set of strategies to evaluate the
forward operator via more efficient optimization schemes212121We note that,
the forward operator defined by these optimization schemes are often _non-
continuous_ —the resulting solutions are not guaranteed to vary smoothly with
changes in the initial state. This makes the forward operator strictly
generative.. Although there are few studies on the direct optimization of the
OK energy, quasi-Newton methods, such as Anderson acceleration [117] and semi-
implicit Seidel schemes [118], are widely used in context of the SCFT model,
which is also formulated as a free energy minimization problem. Similarly, a
modified Newton method for the direct minimization of the OK energy was
recently proposed [114] and shown to be typically three orders of magnitude
faster than the gradient flow approach. The latter Newton method is adopted as
the optimization procedure in our numerical case study.

Given the presence of many local minimizers for OK energy, an approach is
needed to model the random initial state. A commonly adopted initial state is
the realization of a uniform white noise random field [23, 110, 111]. A more
rigorous approach that uses transformed Gaussian random fields is proposed in
[114]. These random fields play the role of the  _auxiliary variable_ $Z$ in
the forward operator. When coupled with the optimization procedure, the random
fields produce a distribution of equilibrium patterns. As described in Section
2.2, a physically consistent model should satisfy the equality in distribution
in (3), meaning that the probability associated with each model solution
matches the probability of realizing the corresponding equilibrium pattern in
an experiment. The validation of these model assumptions underlying the
forward operator is an important topic that we leave for future work.

We provide a description of our model for the auxiliary variable and finite
element discretization of the state space in A and B, respectively. In
conjunction with the variational minimization procedure for the OK model and
its numerical implementation detailed in [114], this completes the definition
of the forward operator in our case study. See Fig. 9 for examples of Di-BCP
patterns generated by the forward operator.

Figure 9: Examples of Di-BCP pattern simulated via various choices of model
parameters ($m,\epsilon,\sigma)$ of the OK model. The model can characterize a
range of morphological phases, interface thicknesses, and strengths of phase
separation.

###  7.3 A model of top-down microscopy characterization

In this subsection, we define the state-to-observable map $\mathcal{J}$ for
top-down microscopy characterization of Di-BCP thin films. We present the
restriction operator $\mathcal{R}$, the imaging operator $\mathcal{I}$, and
the statistical sensor operator $\mathcal{S}$ used in our case study. Recall
that $\mathcal{J}=\mathcal{I}\circ\mathcal{R}$, as defined in Section 2.1.

We assume a set of experimental scenarios, in which microscope images are
taken in fields of view of various sizes, mimicking the effect of varying
levels of magnification. In terms of the characterization model, we first
define a reference field of view $\mathcal{D}=[0,L]^{2}$, where $L$ denotes
its side length, corresponds to a relatively large region of the specimen top
surface. Then, we define the restrictions to its centered subdomains
$R_{j}(\mathcal{D})$, where $j=2,4,\dots,10,20,\dots 90$ denotes the
percentage of the area of $\mathcal{D}$ retained by the restriction. The
restriction operator is thus defined by $\mathcal{R}_{j}:H^{1}(\mathcal{D})\to
H^{1}(R_{j}(\mathcal{D}))$. The effect of these restriction operators, i.e.,
varying levels of magnification, on model calibration will be quantitatively
studied in the next section. The numerical implementation of the restriction
operators is discussed B and a numerical example of the synthetic data
generated by these restriction operators is provided in Fig. 10.

reference field of view | $70\%$ | $40\%$ | $10\%$
---|---|---|---
|  |  |
Figure 10: A set of synthetic top-down microscope images of a Di-BCP thin film
at different levels of magnification, specified by the restriction operator
$\mathcal{R}_{j}$ to the centered subdomain $R_{j}(\mathcal{D})$ with $j\%$,
$j=100,70,40,10$ (from _left_ to _right_), of the area of the reference field
of view $\mathcal{D}$.

For the microscopy imaging operator $\mathcal{I}$, we only consider the
blurring effect, which is often modeled as a point-spread function (PSF) in
the context of electron or optical microscopy. The PSF is a linear integral
operator which degrades the fine-scale features in the observed specimen to
resemble lens-aberrations. Numerous PSFs for microscopy imaging have been
considered in the literature [119, 120, 121], such as Airy disks, Debye
diffraction integral, and Gaussian convolution integral. For our numerical
case study, we employ the convolution with a Gaussian kernel222222We apply
zero-padding to extend $u$ beyond the boundary of the reference field of view
$\partial\mathcal{D}$.as a representative prototype:

$\mathcal{I}(u)=\int_{\mathcal{R}(\mathcal{D})}\frac{1}{\sqrt{2}\sigma_{\textrm{blur}}}\exp{\bigg{(}-\frac{1}{2\sigma_{\textrm{blur}}^{2}}\|\bm{s}_{1}-\bm{s}_{2}\|^{2}\bigg{)}}u(\bm{s}_{1})\textrm{d}\bm{s}_{1},$
(32)

where $\mathcal{I}\colon H^{1}(\mathcal{R}(\mathcal{D}))\to
H^{1}(\mathcal{R}(\mathcal{D}))$ and $\sigma_{\textrm{blur}}$ is a blurring
strength parameter.

Since the sensor noise in microscope images are generally ‘fine-grained’, we
assume there is i.i.d. additive noise for each pixel. Other authors have
considered Poisson, Gaussian or mixed Poisson–Gaussian noise random variable
in a combination of additive and multiplicative form [122, 54, 120]; we employ
additive Gaussian noise as a prototype for simplicity, although more complex
and accurate models could have been employed within the introduced inference
framework without additional difficulty.

In our study, we assume a $256\times 256$ grid array of sensors. The sensor
signal is mathematically defined as local averages in grid cells, denoted as
$\mathcal{D}_{jk}$, centered at the sensor positions:

$\left(\bm{\mathcal{D}}(o)\right)_{jk}=\left|\mathcal{D}_{jk}\right|^{-1}\int_{\mathcal{D}_{jk}}o(\bm{s})\>\textrm{d}\bm{s},\quad
j,k\in(1,2,\ldots,256).$ (33)

We thus arrive at the following characterization model:

$\mathcal{S}(\mathcal{J}(u),\bm{N})=\bm{\mathcal{D}}(\mathcal{J}(u))+\bm{N},\quad\bm{N}\overset{\textrm{i.i.d.}}{\sim}\mathcal{N}_{jk}(0,\sigma^{2}_{\bm{N}}),$
(34)

where $\mathcal{J}=\mathcal{I}\circ\mathcal{R}$ as defined above.

In the next section, we examine, jointly and separately, the effect of various
levels of blurring, $\sigma_{\text{blur}}=(\text{None},5L\times
10^{-3},L\times 10^{-2})$, and/or the effect of various sensor noise variances
$\sigma^{2}_{\bm{N}}$, based on a signal-to-noise (SNR) ratio of
$(\infty,5.0,3.0,2.5)$, on the informativeness of the data. SNR is related to
the sensor noise variance $\sigma^{2}_{\bm{N}}$ by, naively,
$\sigma_{\bm{N}}^{-1}=\textrm{SNR}$.

###  7.4 Summary statistics for top-down microscope images of diblock
copolymer thin films

In this section, several summary statistics for top-down microscopy
characterization images of Di-BCP thin films, i.e.
$\bm{Q}=\mathcal{Q}(\bm{Y})\in\mathbb{R}^{q}$, are proposed. First, in Section
7.4.1 we introduce summary statistics defined by the composition of energy
functionals and an $H^{1}$-interpolation operator. Then, in Section 7.4.2, we
present summary statistics based on the magnitude spectrum of these images
obtained by discrete Fourier transform. For both cases, we discuss connections
between these summary statistics and features in the image data, and provide
reasoning for why we expect these summary statistics to be informative of the
model parameters.

####  7.4.1 Energy-based summary statistics

We introduce five summary statistics, $Q_{i}=\mathcal{Q}_{i}(\bm{Y})$,
$i=0,1,\dots,4$, where $\mathcal{Q}_{i}$ is defined by the composition of an
energy functional and an $H^{1}$ interpolation operator. The $H^{1}$
interpolation operator, denoted as $\Pi$, interpolate the image
$\bm{y}\in\mathbb{R}^{y}$ into the function space $H^{1}(\mathcal{D}_{u})$
using a linear finite element basis, where $\mathcal{D}_{u}$ is a unit square
domain232323We note that scalings of the summary statistics do not change the
EIG between $\bm{X}$ and $Q_{i}$, since EIG is invariant under invertible
transformations of $\bm{Y}$. The scalings can be helpful to adjust the
magnitude of summary statistics to the same order.. The details of the
interpolation operator for our numerical case study are given in B. The energy
functionals act on the $H^{1}$-interpolated image data. Their names,
notations, and mathematical forms are listed in Tab. 1.

Notation | Name | Functional
---|---|---
$Q_{0},\widehat{M}$ | Empirical spatial average | $\displaystyle\mathcal{M}(\widehat{u})=\int_{\mathcal{D}_{u}}\widehat{u}\;\textrm{d}\bm{s}$
$Q_{1}$ | Double well | $\displaystyle\mathcal{D}^{\textrm{OK}}(\widehat{u})=\int_{\mathcal{D}_{u}}(1-\widehat{u}^{2})^{2}/4\;\textrm{d}\bm{s}$
$Q_{2}$ | Interfacial | $\displaystyle\mathcal{I}^{\textrm{OK}}(\widehat{u})=\int_{\mathcal{D}_{u}}|\nabla\widehat{u}|^{2}/2\;\textrm{d}\bm{s}$
$Q_{3}$ | Nonlocal | $\displaystyle\mathcal{N}^{\textrm{OK}}(\widehat{u})=\int_{\mathcal{D}_{u}}(\widehat{u}-\mathcal{M}(\widehat{u}))(-\triangle_{N})^{-1}(\widehat{u}-\mathcal{M}(\widehat{u}))/2\;\textrm{d}\bm{s}$
$Q_{4}$ | Total variation | $\displaystyle\mathcal{TV}(\widehat{u})=2(1-|\mathcal{M}(\widehat{u})|)^{-1}\int_{\mathcal{D}_{u}}(\nabla\widehat{u}\cdot\nabla\widehat{u})^{1/2}\;\textrm{d}\bm{s}$
Table 1: The notations, names, and corresponding energy functionals of the
proposed energy-based summary statistics.

Three summary statistics, $Q_{1},Q_{2},Q_{3}$, are inspired by the OK free
energy. They are defined using the functionals in the OK energy, i.e.,
$\mathcal{D}^{\textrm{OK}}$, $\mathcal{I}^{\textrm{OK}}$, and
$\mathcal{N}^{\textrm{OK}}$, evaluated at the interpolated image data random
variable $\Pi(\bm{Y})$. The motivation behind choosing these summary
statistics are our perceptions that (i) the magnitude of the three energies
are sensitive to the model parameters and relatively insensitive to the
auxiliary variable, and (ii) they characterize informative geometric features
of Di-BCP patterns, such the total area of phase interface, degree of phase
separation, etc. Since the microscopy characterization visually preserves the
geometry of phase interfaces, we suspect that the negative impact of sensor
blurring and noise on the informativeness of these summary statistics for
inferring the model parameters is relatively minor. Therefore, we hypothesize
that the parameters are well-informed by these energy-based summary
statistics.

We also define two additional energy-based summary statistics, given by the
empirical spatial average,
$Q_{0}=\mathcal{Q}_{0}(\bm{Y})\coloneqq(\mathcal{M}\circ\Pi)(\bm{Y})$, and a
total-variation-like functional,
$Q_{4}=\mathcal{Q}_{4}(\bm{Y})\coloneqq(\mathcal{TV}\circ\Pi)(\bm{Y})$. We
expect the empirical spatial average to be informative of the material
parameter $f$ that represents the segment number ratio of monomer A. For the
OK model, $f$ is directly related to the model parameter $m$ through (28). The
total-variation functional $\mathcal{TV}$ formally captures the average ratio
between the interfacial length and the inclusion domain size, such as the
radius of spots for spot patterns, of Di-BCP patterns. It is also relatively
insensitive to the sensor noise with edge preserving property, which makes it
popular in image denoising problems [123].

We visualize integrands of the functionals defining the energy-based summary
statistics for synthetic data with lamella and spot latent patterns in Fig.
11.

Synthetic image $\bm{y}^{*}$ | $\mathcal{Q}_{1}(\bm{y}^{*})$ | $\mathcal{Q}_{2}(\bm{y}^{*})$ |  $\mathcal{Q}_{3}(\bm{y}^{*})$ (log scale) | $\mathcal{Q}_{4}(\bm{y}^{*})$
---|---|---|---|---
|  |  |  |
|  |  |  |
Figure 11: Two synthetic top-down microscopy characterization images,
$\bm{y}^{*}$, of Di-BCP thin films and the integrands for functionals,
$\mathcal{Q}_{i}$, $i=1,\dots,4$, evaluated at the images. These functionals
are the energy-based summary statistics shown in Tab. 1.

####  7.4.2 Fourier-based summary statistics

Di-BCP thin film equilibrium patterns are well characterized by their
intrinsic length scales. For example, the _interfacial thickness_ and the
_periodicity length_ of Di-BCP patterns are often extracted from microscopy
images for quantitative analysis [124, 41]. In particular, the azimuthally-
averaged magnitude spectrum is widely used to process and analyze the
microscopy characterization images of BCP equilibrium patterns [125, 126]. In
the absence of external guidance or boundary geometry effects, there is no
orientational preference in the latent patterns, which suggests that the
orientational averaging of the magnitude spectrum preserves most of the
information for inferring the model parameters. We thus hypothesize that the
azimuthally-averaged magnitude spectrum of microscopy image data and its
summary statistics may be effective for model parameter inference.

We note that the azimuthally-averaged magnitude spectrum is similar to
_scattering intensity profiles_ in small angle X-ray scattering. The
scattering intensity profiles are modeled as noisy Fourier domain
representations of the latent states [127, 39]. In this sense, the use of
azimuthally-averaged magnitude spectrum and its related summary statistics for
parameter inference and EIG computation can be seen as a rough prototype for
the same task performed with small angle X-ray scattering data242424However,
we acknowledge several important differences. First, the X-ray scattering
imaging operator involves a Fourier transform and the sensor noise is applied
on top of the observables, whereas the magnitude spectrum of microscope images
already includes the sensor noise and the blurring effect. Secondly, the
imaging operator for the X-ray scattering experiment depends highly on the
experimental set up. Changes in the experimental design variables, such as
beam stop positions and X-ray incident angles, can lead to dramatic changes in
the observables..

We now define the procedure for obtaining the azimuthally-averaged magnitude
spectrum from an image $\bm{y}^{*}\in\mathbb{R}^{N\times M}$. The discrete
Fourier transform of the image outputs a complex Fourier coefficient matrix
$\widehat{\bm{y}^{*}}\in\mathbb{C}^{N\times M}$. We first translate the
entries to center the zero frequency and examine the magnitude spectrum:

$\mathcal{P}(\bm{y}^{*})=|\widehat{\bm{y}^{*}}|\in\mathbb{R}^{N\times M}.$
(35)

For azimuthal averaging of the magnitude spectrum, a change-of-variables to
radial coordinates centered at the zero frequency is performed. In the radial
coordinate, the zero-centered magnitude spectrum has $N_{r}$ unique radial
frequencies, $(r_{1},\dots,r_{N_{r}})\in{(0,r_{\textrm{max}}]}^{N_{r}}$. We
only study the magnitude spectrum for frequencies no larger than some right-
closed intervals, as higher radial frequencies are not well-resolved and
typically uninformative. Finally, we obtain the azimuthally-averaged magnitude
spectrum, $\mathcal{P}_{a}\colon\mathbb{R}^{N\times M}\to\mathbb{R}^{N_{r}}$,
by averaging over orientations $\theta$. Specifically, the average magnitude
$\mathcal{P}_{a}(\bm{y}^{*})(r_{k})$ is computed for $N_{\theta}(r_{k})$
number of entries in $\mathcal{P}(\bm{y}^{*})$ with the identical radial
frequency $r_{k}$:

$\mathcal{P}_{a}(\bm{y}^{*})(r_{k})=\frac{1}{N_{\theta}(r_{k})}\sum_{j=1}^{N_{\theta}(r_{k})}\mathcal{P}(\bm{y}^{*})(r_{k},\theta_{j}),\quad
k=1,2,\dots,N_{r}.$ (36)

We hypothesize that not much, if any, information that is relevant toward
inferring the length scale dependent parameters is lost after the azimuthal
averaging.

The number of unique radial frequencies $N_{r}$ is typically large. To further
reduce the dimension of the azimuthally-averaged magnitude spectrum data, we
also apply averaging over each unit interval of the radial frequency to obtain
$\lfloor r_{\text{max}}\rfloor$-dimensional averaged spectrum
$\tilde{\mathcal{P}}_{a}(\bm{y}^{*})$:

$\left(\tilde{\mathcal{P}}_{a}(\bm{y}^{*})\right)_{j}=\frac{1}{\sum_{j-1<r_{k}\leq
j}N_{\theta}(r_{k})}\sum_{j-1<r_{k}\leq
j}\sum_{i=1}^{N_{\theta}(r_{k})}\mathcal{P}(\bm{y}^{*})(r_{k},\theta_{i})\;,\quad
j=1,2,\dots,\lfloor r_{\text{max}}\rfloor.$ (37)

As shown in Fig. 12, the zero-centered magnitude spectrum of synthetic
microscope images have ‘rings’ of frequencies with dominant magnitude values.
They represent the periodicity length of the latent Di-BCP pattern.
Furthermore, the overall shape of the magnitude spectrum distribution is
sensitive to the change in model parameters and relatively insensitive to the
auxiliary variable. Therefore, we further extract summary statistics such as
dominant radial frequency, or some quantitative descriptions of the shape of
the azimuthally-averaged magnitude spectrum, such as moments of various
orders. Moreover, using $\bm{Q}_{a}=\tilde{\mathcal{P}}_{a}(\bm{Y})$ as
summary statistics of microscopy data allows us to further investigate and
compare the informativeness of the averaged magnitude spectrum in different
ranges of the radial frequency for learning the model parameters.

Synthetic image $\bm{y}^{*}$ | $\ln(\mathcal{P}(\bm{y}^{*}))$ | $\mathcal{P}_{a}(\bm{y}^{*})$ | $\tilde{\mathcal{P}}_{a}(\bm{y}^{*})$
---|---|---|---
|  |  |
|  |  |
Figure 12: Two synthetic top-down microscopy characterization images,
$\bm{y}^{*}$, of Di-BCP thin films, their log-scale magnitude spectrum,
$\ln(\mathcal{P}(\bm{y}^{*}))$, their azimuthally-averaged magnitude spectrum,
$\mathcal{P}_{a}(\bm{y}^{*})$, with $r_{\text{max}}=50$, and their interval-
averaged magnitude spectrum, $\tilde{\mathcal{P}}_{a}(\bm{y}^{*})$. The
synthetic images are generated by the OK model at the same model parameter.

###  7.5 Prior distribution

We now define the prior distribution $\nu_{\bm{X}}$, with density
$\pi_{\bm{X}}$, for the parameters used in our case study. The parameters of
interest, again, are $\bm{x}=(m,\epsilon,\sigma)$, where $m\in(-1,1)$,
$\epsilon\in\mathbb{R}_{+}$, and $\sigma\in\mathbb{R}_{+}$. We label the
individual parameter random variables as
$\bm{X}=(M,E,\Sigma)\sim\nu_{\bm{X}}$.

$\nu_{\bm{X}}$ is constructed based on hypotheses on the experimental setting
related to the parameters of the OK model. These assumptions put restrictions
on the feasible or admissible portions in parameter space. We denote this
admissible set of parameters by $\mathcal{A}\subsetneq\mathbb{R}_{+}^{3}$.
Then, under the principal of maximum entropy [128], we consider a uniform
prior distribution252525We note that if one’s experimental setup provides
additional information about the model parameters, the prior distribution may
be much more concentrated in regions of parameter space. over this admissible
set, whose density is given by:

$\pi_{\bm{X}}(\bm{x})=\begin{cases}1/|\mathcal{A}|\quad&\bm{x}\in\mathcal{A}\\\
0\quad&\bm{x}\not\in\mathcal{A}.\end{cases}$ (38)

Our uniform prior is defined _hierarchically_ , that is, the parameters are
defined conditionally on the others given the variable ordering
$(M,E,\Sigma)$. This hierarchical distribution implies a factorization of the
prior density:

$\displaystyle\pi_{\bm{X}}(m,\epsilon,\sigma)$
$\displaystyle=\pi_{M}(m)\pi_{E|M}(\epsilon|m)\pi_{\Sigma|M,E}(\sigma|m,\epsilon).$
(39)

The marginal density $\pi_{M}$, for the spatial average parameter $M$, known
as the _hyper-parameter prior_ or _hyper-prior_ , is uniform and is
independent of the other two parameter random variables $(E,\Sigma)$; the
conditional densities $\pi_{E|M}(\cdot|m)$ and
$\pi_{\Sigma|M,E}(\cdot|m,\epsilon)$ are also uniform. We implicitly define
the factors in (39) by describing the constraints that form the boundary
$\partial\mathcal{A}$ of the admissible set $\mathcal{A}$.

First, we consider constraints on the spatial average parameter $M$. In our
experimental setting, we assume that the latent material states of observed
images exhibit phase separation with a high probability. Note that the
homogeneously mixed state are trivial critical points of the OK energy. In
this case, the length scale related parameters $E$ and $\Sigma$ do not play a
role in the forward operator, and thus not identifiable from data. To define
constraint boundaries for the hyper-prior $\pi_{M}$ based on this assumption,
we consider a result from the analysis of the OK energy functional:

######  Observation 1 (Constant state: Global minimum or metastable).

The constant state $u=m$ with no phase seperation is the global minimum of
$\mathcal{E}_{\bm{x}}^{\textrm{OK}}$ (25) whenever $|m|>1/\sqrt{3}$.
Otherwise, it is a local minimum, i.e. a metasable state. [110, 111, 114].

It is often assumed that global minimum states of
$\mathcal{E}_{\bm{x}}^{\textrm{OK}}$ are relatively more stable
configurations, i.e., more probable. Therefore, we exclude portion of
parameter space in which the constant states are global minima, and define the
admissible set for the hyper-prior for $M$ as $(0,1/\sqrt{3})$, and the
marginal prior density for $M$ as

$\pi_{M}=\mathcal{U}(0,1/\sqrt{3}),$ (40)

where we consider only positive numbers in our study, without loss of
generality, due to symmetry of the forward operator about $m=0$.

To define admissibility constraints for the conditional density of $E$ and
$\Sigma$ given the parameter $M$, we employ two assumptions on the equilibrium
patterns encountered in our experimental scenario. These assumptions put
restrictions on the size of the field of view of the microscope relative to
different length scales of latent Di-BCP patterns. We note that these
assumptions are based on a rescaled computational domain
$\mathcal{D}_{c}=[0,L]^{2}$, $L=40$ , with the material and model parameters
rescaled according to the relations in (28). See a discussion on the
computational domain in B.

######  Assumption 1 (Resolving phase interface).

We assume the phase interface of the latent Di-BCP patterns can be
sufficiently resolved if given little or no image corruptions, which
correspond to $\epsilon>\epsilon_{\textrm{min}}=0.01$.

######  Assumption 2 (Observing relatively large surfaces).

We assume the periodicity of Di-BCP patterns can be sufficiently captured in
the image data [110, 111], if given little or no image corruptions:

$L\gg O\left(\left(\epsilon/\sigma\right)^{1/3}\right)\implies\sigma>\epsilon
C^{3}/L^{3},\quad C\gg 1,$ (41)

and we thus assume $\sigma>\sigma_{\textrm{min}}=0.01$.

To complete the admissibility constraints, we emulate the work of [110, 111],
in which the authors characterize a _phase diagram_ for the OK model that
describes the morphology of equilibrium states for given model parameters
$(\epsilon,\sigma)$. By analyzing the (linear) stability of global minimizers
of the OK energy under periodic boundary conditions through Fourier
expansions, they derive “order/disorder transition (ODT) curves” that
partition the parameter space of $(\epsilon,\sigma)$ into: (i) a region in
which the OK energy functional possess the constant states as the global
minimizers, (ii) a region in which the constant states are _not_ the global
minimizers, and (iii) a transition region where the analysis is not
applicable. It is hypothesized in [111], based on numerical evidence, that the
_actual_ region in parameter space for which the constant states are _not_ the
global minimizers is actually larger, therefore expanding the ODT curve
further into the transition region.

We follow these authors in producing a _numerically discovered empirical
boundary_ for the ODT curve based on the OK model.262626While we consider the
homogeneous Neumann boundary conditions (7.1), for which the analysis in [111]
does not directly apply—since the Fourier expansion arguments only apply to
the case of periodic boundary conditions—our empirical ODT curves have similar
qualitative behavior as in the authors’ work. The result of the numerical
procedure is two parameterized curved boundary surfaces
$\sigma_{\textrm{max}}(\epsilon;m)$ and $\sigma_{\textrm{min}}(\epsilon;m)$
which define the upper and lower boundary for the parameter $\sigma$
conditioned on the values for m and $\epsilon$. A visualization of the curved
line “slices" of the surfaces, for three values of $m$ are provided in Fig.
13. Intersecting these surfaces with the fixed bounds
$\epsilon_{\textrm{min}}=0.01$, $\sigma_{\textrm{min}}=0.01$ from Assumptions
1, 2, and $M\in(0,1/\sqrt{3})$ from Observation 1, produces the boundary
$\partial\mathcal{A}$.

Figure 13: (_left_) The admissible set $\mathcal{A}\subset\mathbb{R}^{3}$.
(_Right_) The boundary $\partial\mathcal{A}$ of the admissible set for
$(\epsilon,\sigma)$ for three choices of the order parameter spatial average
$m=(0,0.3,1/\sqrt{3})$. The prior distribution for the parameters $(E,\Sigma)$
conditioned on the hyper-parameter $M$, $\nu_{E,\Sigma|M}(\cdot,\cdot|m)$, is
uniformly distributed over the admissible set. We also plot 10,000
independently drawn samples from each conditional distribution.

##  8 Case study: Numerical results

In this section, we present the numerical results for our case study. These
results demonstrate the insights one can extract from performing calibration
tasks via the measure transport approach. In Section 8.1, we provide
representative posterior samples to visually demonstrate the reduction of the
model parameter uncertainty after performing the model calibration with the
energy-based summary statistics. Samples of the observables, i.e., the latent
Di-BCP patterns, at these parameter values are also generated and plotted to
qualitatively validate the inference results. In Section 8.2, we present
results on estimating EIG of model parameters from various combinations of the
energy-based summary statistics with varying levels of noise and blurring
corruptions. The similar EIG estimation results from the Fourier-based summary
statistics are provided in Section 8.3. In Section 8.4, we study the impact of
various levels of magnification for microscopy imaging on the calibration
results. Finally, in Section 8.5, we study the EIG estimation for the
interfacial and nonlocal parameters $(E,\Sigma)$ when the spatial average
parameter $M=m$ is known.

To perform our numerical experiments, we first collect $n=10^{5}$ samples from
the joint distribution of parameters and the model-generated data
$\\{(\bm{x}^{i},\bm{y}^{i})\\}_{i=1}^{n}$, as described in Section 4.2. Then,
the samples of model-generated data are used to compute samples of the
proposed summary statistics
$\\{(\bm{x}^{i},\bm{q}^{i}=\mathcal{Q}(\bm{y}^{i}))\\}_{i=1}^{n}$. We use
$60\%$ of the samples as a training set for estimating the posterior
probability densities using the triangular transport map approach introduced
in Section 5 and use the remaining samples to estimate the EIG using the
likelihood-free estimator introduced in Section 6.

A very small subset ($\sim\thinspace 200$) of the samples are numerically
problematic due to (i) failure to resolve phase interface geometries in the
observable and (ii) failure to evaluate the forward operator within a
reasonable residual tolerance. In both cases, the numerical anomalies
primarily occur near the boundary of the admissible set $\partial\mathcal{A}$
and are caused by relatively coarse spatial discretization of the state space.
In our numerical case study, we choose to discard these samples. To detect
these numerically problematic samples, an empirical approach is adopted based
on the value of the KL divergence between
$\tilde{\mathcal{P}}_{a}(\bm{y}^{i})$ and a reference set of values. We
observe that the numerically problematic samples invariably possess a large
value of such a KL-divergence. This ansatz is used an automated criterion to
discard a majority of these numerically problematic samples. Additionally, to
deal with outlier samples, we compute the $0.05$ and $0.995$ quantiles of the
log posterior evaluations
$\widehat{\pi}_{\bm{X}|\bm{Q}}(\bm{x}^{i}|\bm{q}^{i})$ and discard samples
that are outside of these quantiles for the EIG estimation.

###  8.1 Posterior and posterior predictives

To visualize posterior predictive samples, we first generate representative
posterior quadrature points. In particular, we take Gauss–Hermite quadrature
points $\bm{r}_{\textrm{GH}}$ which are used for accurately integrating
functions with respect to the standard normal distribution, i.e, the reference
distribution $\nu_{\bm{R}_{2}}$, and transform them through the inverse of the
posterior transport map, i.e.
$\widehat{S}^{\mathcal{X}}(\bm{q}^{*},\cdot)^{-1}(\bm{r}_{\textrm{GH}})$, as
described in (20). These points are depicted within the admissible set for the
prior in Tab. 2. While these points, paired with appropriate weights, can be
used for both visualizing and accurately integrating functions with respect to
the posterior distribution, we only present numerical results on the
visualization of the model predictions at posterior samples in this study. We
note that the points include the inverse transformation of the center
quadrature point,
$\bm{x}_{\bm{0}}=\widehat{S}^{\mathcal{X}}(\bm{q}^{*},\cdot)^{-1}(\bm{0})$, as
a representative sample that is typically in the region with highest posterior
probability2727footnotemark: 27. In addition, we include the transformation of
quadrature points in the tails of each reference marginal $\mathcal{N}(0,1)$,
in order to study the ‘tails’, or regions of low probability, of the
posterior. These points are labeled as $\bm{x}_{\pm\delta
m}\coloneqq\widehat{S}^{\mathcal{X}}(\bm{q}^{*},\cdot)^{-1}(\bm{0}\pm\delta
m)$,
$\bm{x}_{\pm\delta\epsilon}:=\widehat{S}^{\mathcal{X}}(\bm{q}^{*},\cdot)^{-1}(\bm{0}\pm\delta\epsilon)$,
and