smmfstream#3dict/Type/Filespec/F (#3)/UF (#3)/EF ¡¡/F
¿¿dict/Type/RichMediaInstance/Subtype/Video/Asset /Params
¡¡/FlashVars(source=#3&skin=SkinOverAllNoFullNoCaption.swfskinAutoHide=trueskinBackgroundColor=0x5F5F5FskinBackgroundAlpha=0)¿¿dict/Type/RichMediaConfiguration/Subtype/Video/Instances
[]dict/Type/RichMediaContent/Assets ¡¡/Names [(#3) ]¿¿/Configurations
[]dict/Activation ¡¡/Condition/#1PVXA/Presentation
¡¡/Style/Embedded¿¿¿¿/Deactivation
¡¡/Condition/PI¿¿1113.33337pt6.94444pt1.94443pt/Subtype/RichMedia/BS ¡¡/W
0/S/S¿¿/Contents (embedded video file:#3)/NM (rma:#3)/AP ¡¡/N
¿¿/RichMediaSettings /RichMediaContent #2

# Anisotropies in the spatial distribution and kinematics of dwarf galaxies in
the Local Group and beyond

Isabel M.E. Santos-Santos,1 Julio F. Navarro,2 Alan McConnachie3,2
1Institute for Computational Cosmology, Department of Physics, Durham
University, South Road, Durham, DH1 3LE, UK
2Department of Physics and Astronomy, University of Victoria, Victoria, BC V8P
5C2, Canada
3NRC Herzberg Astronomy and Astrophysics, 5071 West Saanich Road, Victoria,
BC, V9E 2E7, Canada E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

The Local Group (LG) of galaxies is dominated by the Andromeda (M31) and Milky
Way (MW) pair, a configuration which suggests that the mass distribution in
the LG and its surroundings should be highly anisotropic. We use the APOSTLE
suite of cosmological hydrodynamical simulations to examine how this
anisotropy manifests on the spatial distribution and kinematics of dwarf
galaxies out to a distance of $\sim 3$ Mpc from the MW. The simulations
indicate a clear preference for dwarfs to be located close to the axis defined
by the MW-M31 direction, even for dwarfs in the LG periphery (LGP), defined as
those at distances $1.25<d_{\rm MW}/$Mpc$<3$ from the MW. The “Hubble flow” in
the periphery is also affected; at fixed distance from the MW the mean
recession speed, $\langle V_{\rm rad}\rangle$, varies with angular distance to
M31, peaking in the anti-M31 direction and reaching a minimum behind M31. The
combined M31-MW mass decelerates the local expansion; the LG “turnaround
radius” (i.e., where $\langle V_{\rm rad}\rangle=0$) in APOSTLE is located at
$r\sim 1.25$ Mpc from the LG barycentre and the pure Hubble flow (i.e., where
$\langle V_{\rm rad}\rangle$ becomes comparable to $H_{0}*d$) is not reached
out to at least $r\sim 3$ Mpc. The predicted flow is very cold, with a
barycentric dispersion of less than $\sim 40$ km/s. A comparison of these
predicted features with existing observations gives mixed results. There is
clear observational evidence for an angular anisotropy in recession velocities
around the LGP, but there is little evidence in the spatial distribution of
LGP dwarfs for a preferred direction along the MW-M31 direction. The observed
local Hubble flow is also somewhat peculiar. Although the “coldness” of the
flow is consistent with the simulations, it is significantly less decelerated:
relative to the MW, on average, all galaxies beyond $d_{\rm MW}\sim 1.25$ Mpc
seem to be already on a pure Hubble flow. We argue that these oddities may
result at least in part from incompleteness and inhomogeneous sky coverage,
but a full explanation may need to await the completion of deep all-sky
surveys able to fill the gaps in our current inventory of nearby dwarfs.

###### keywords:

galaxies: dwarf – galaxies: Local Group – cosmology: dark energy

††pubyear: 2023††pagerange: Anisotropies in the spatial distribution and
kinematics of dwarf galaxies in the Local Group and beyond–Anisotropies in the
spatial distribution and kinematics of dwarf galaxies in the Local Group and
beyond

##  1 Introduction

In the Lambda Cold Dark Matter cosmogony (LCDM), the current paradigm for
structure formation, the Local Group (LG) of galaxies is thought to arise as
two relatively isolated massive dark matter halos, the hosts of the Milky Way
(MW) and the Andromeda (M31) galaxies, detach from the universal expansion
under the influence of their mutual gravity, turn around, and start heading
towards each other on a nearly radial orbit. Achieving, on first approach, the
observed relative radial velocity ($\sim-109$ km/s) at the current MW-M31
separation ($\sim 774$ kpc) in roughly $14$ Gyr (the age of the Universe)
suggests that the MW-M31 system turned around a few Gyrs ago after reaching a
maximum separation of $\sim 1.1$ Mpc (Fattahi et al., 2016), and that their
combined mass is at least a few times $10^{12}\,M_{\odot}$ (the “timing
argument”, e.g., Kahn & Woltjer, 1959; Li & White, 2008).

Although the combined mass estimate is in reasonable agreement with current
estimates of the MW and M31 virial111We define the virial boundary of a system
as the radius where the mean enclosed density is 200$\times$ the critical
density for closure, and refer to virial quantities with the subscript ‘200’.
masses (see, e.g., Cautun et al., 2020; Patel & Mandel, 2023, and references
therein), other features of the LG formation scenario outlined above are
seemingly at odds with LG observations.

For example, relative to the Milky Way, the current LG turnaround radius
(i.e., where the mean recession velocity is $\langle V_{\rm rad}\rangle\approx
0$) should be substantially farther than $1.1$ Mpc, the expected turnaround
radius of M31. Fattahi et al. (2016), for example, estimate that the
turnaround radius at present could be as large as $1.7$ Mpc from the MW.
Galaxies just inside that radius should have already turned around, and have
today mainly negative radial velocities. Just outside turnaround, on the other
hand, galaxies should still be receding on average, but with a substantially
decelerated Hubble flow. These two robust predictions are apparently in
contrast with observations: all known dwarf galaxies beyond $d_{\rm MW}\sim
1.25$ Mpc from the Milky Way are receding from us, following an apparently
undecelerated, “pure” Hubble flow.

Before analyzing this further, we note that the true boundaries of the Local
Group are somewhat ill-defined. A common definition, which we follow here,
defines LG galaxies as those within the current turnaround radius, estimated
empirically at around $\sim 1$ Mpc from the MW-M31 barycenter (see; e.g.,
McConnachie, 2012). Galaxies just outside this radius, although receding from
the LG barycenter, may still be bound to the LG. We shall hereafter refer to
galaxies in the LG periphery, i.e., those at $1.25<d_{\rm MW}/$Mpc$<3$, as
“LGP dwarfs”.

Returning to the apparent disagreements mentioned above, one reason may be
incompleteness in the inventory of dwarfs in the LG and its periphery, as well
as their patchy distribution across the volume. Indeed, our inventory of
nearby dwarfs is likely woefully incomplete, as discussed by Fattahi et al.
(2020), who conclude that as many as $\sim 50$ dwarfs as massive as the Draco
dwarf spheroidal could be missing from our current inventory of LG and LGP
members.

This incompleteness should be taken carefully into account when examining the
Hubble flow around the Milky Way, as well as its dispersion. Because the Milky
Way is offset from the LG barycentre, recession velocities at given distance
are expected to depend on sky position, reaching a maximum in the anti-M31
direction and a minimum behind M31. (This assumes that the LG effect on the
local recession velocity field is more or less symmetric relative to the LG
barycentre, which should lie somewhere midway between MW and M31.)

Of all $31$ known LGP dwarfs only $2$ are within $45$ deg of M31 and $4$ in
the opposite anti-M31 direction. Such patchy coverage may therefore have a
strong effect not only on how decelerated the local Hubble flow may appear,
but also on estimates of its dynamical “coldness”, an issue that has been
discussed quite extensively in the literature, with conflicting claims
(Sandage et al., 1972; Schlegel et al., 1994; Macciò et al., 2005, and
references therein).

Cosmological simulations that capture the particular dynamical configuration
of the MW and M31 may also offer guidance regarding where LGP dwarfs missing
from our current inventory might be located. Fattahi et al. (2020), for
example, noted that many of them should be located behind and around M31. This
anisotropy suggests that distant dwarf galaxy searches may be substantially
more fruitful in some regions of the sky relative to others. Indeed, prior
work has suggested that LGP dwarfs should be preferentially aligned with the
MW–M31 axis, where the effects of the quadrupole of the mass distribution are
maximized (Peñarrubia et al., 2014).

Such guidance has already proven useful in the case of the search for new MW
satellites in the vicinity of the Magellanic Clouds. Indeed, Sales et al.
(2011) used cosmological N-body simulations to predict that the surroundings
of the Clouds should “prove a fertile hunting ground for faint, previously
unnoticed MW satellites”, a prediction that became spectacularly true with the
discovery of a number of nearby dwarfs in the DES survey, one of the first to
target large fractions of the southern sky to magnitudes deep enough to allow
for the identification of new MW satellites and LG members (Bechtol et al.,
2015).

We revisit these issues here, using simulations from the APOSTLE222APOSTLE
stands for “A Project Of Simulating The Local Environment” (Sawala et al.,
2016; Fattahi et al., 2016). project. Our main goal is to characterize
expected anisotropies in the spatial distribution of dwarfs in the Local
Group, with particular emphasis on the spatial and kinematic properties of the
LGP dwarf population, and to compare them with current observations.

We compare simulations and observations mainly in the MW-centric frame, since
the lack of proper motions for most nearby dwarfs means that it is not
possible to transform accurately their velocities to an LG-centric frame. This
approach differs from that commonly adopted in earlier work, where LG-centric
velocities are estimated by simply projecting the observed radial velocities
onto a frame where the dispersion in the Hubble flow of distant galaxies is
minimized (see; e.g., Karachentsev et al., 2009; Peñarrubia et al., 2014, and
references therein). This approach may introduce biases not only because it
neglects tangential velocities, but also because of the aforementioned lack of
homogeneous spatial coverage in our current inventory of nearby dwarfs, which
may compromise the minimization procedure.

This paper is organized as follows. Sec. 2 introduces the numerical
simulations and observational dataset of galaxies used. Sec. 3.1 quantifies
the anisotropies in the spatial distribution of dwarfs in the APOSTLE
simulations, while Sec. 3.2 examines anisotropies in the recession speeds of
dwarfs, with emphasis on the deceleration and dispersion of the local Hubble
flow. Finally, we summarize our results and discuss their implications in Sec.
4.

*figures/video3DLG_longfast.mp4

Figure 1: 3D positions of $M_{*}>10^{5}\,M_{\odot}$ Local Group dwarf galaxies
within a radius of $\sim 3.5$ Mpc from the LG midpoint. Left: APOSTLE,
centered on the MW analog: To enhance visually the shape of the spatial
distribution, we stack the data from the 4 APOSTLE-L1 volumes. Right:
Observational data, centered on the MW. In both panels, the X-axis is aligned
with the MW-M31 direction. Galaxies are colored according to $\alpha$, the
angular distance between a system and the direction to M31; see colorbar. Gray
spherical cones delineate an aperture of $45^{\circ}$ around this axis for
reference.

##  2 Methods

###  2.1 Numerical simulations

We use the APOSTLE project, a suite of “zoom-in” cosmological hydrodynamical
simulations of LG-like environments that include two primary halos with
masses, relative distances, and relative radial and tangential velocities
chosen to be roughly consistent with current observational constraints for the
MW-M31 pair (Fattahi et al., 2016).

In this work we have used the $z=0$ outputs of 4 different APOSTLE volumes
from the highest “L1” resolution level, with initial dark matter and gas
particle masses of $m_{\rm dm}\sim 5\times 10^{4}$ M⊙ and $m_{\rm gas}\sim
1\times 10^{4}$ M⊙, respectively, and a gravitational softening length of
$134$ pc at $z=0$. The average combined virial mass of the MW and M31 pair is
$2.7\times 10^{12}$ M⊙, with mass ratios in the range $0.65$-$0.97$. In each
APOSTLE volume, the zoomed-in region includes a sphere of radius $r\sim 3.5$
Mpc around the midpoint of the MW-M31 pair, which is fully enclosed within the
highest resolution volume.

APOSTLE has been run with the EAGLE galaxy formation code (Schaye et al.,
2015; Crain et al., 2015), which includes subgrid physics prescriptions for
radiative cooling of gas, star formation in gas particles exceeding a
metallicity-dependent density threshold, stellar feedback in the form of
stellar winds, radiation pressure and supernovae, as well as an homogeneous
X-ray/UV background radiation. The model also accounts for supermassive black
hole growth and AGN feedback, but we note that these have negligible effects
in the APOSTLE volume as it is dominated by low-mass, dwarf galaxies.

The APOSTLE project adopts a flat $\Lambda$CDM cosmological model with WMAP-7
parameters (Komatsu et al., 2011): $\Omega_{\rm m}=0.272$;
$\Omega_{\Lambda}=0.728$; $\Omega_{\rm bar}=0.0455$; $H_{0}=100\,h$ km s-1
Mpc-1; $\sigma=0.81$; $h=0.704$.

####  2.1.1 Simulated galaxies

APOSTLE halos were identified using the friends-of-friends (FoF) groupfinding
algorithm (Davis et al., 1985) assuming a linking length of 0.2 times the mean
interparticle separation. Self-bound substructures within FoF groups were then
identified using SUBFIND (Springel et al., 2001).

Luminous galaxies form in APOSTLE at the centre of halos that exceed a
redshift-dependent “critical mass”, set by the UV-ionizing background
(Benitez-Llambay & Frenk, 2020). At $z=0$ this threshold corresponds to a
virial mass of $M_{200}\sim 10^{9}$ M⊙, resolved with $>2\times 10^{4}$ dark
matter particles in APOSTLE-L1 (Pereira-Wilson et al., 2023).

Systems with lower-than-critical mass that host a luminous galaxy are either
halos which were over the critical boundary in the past but whose recent mass
accretion history has been uncharacteristically slow, or the result of tidal
stripping, which may reduce the total dark matter mass of a system that orbits
a more massive host.

In this work we are mainly interested in nearby “field” dwarfs that are not
satellites of either the MW or M31; i.e., those found at $z=0$ outside the
virial radii of the two massive APOSTLE primaries, and within a sphere of
radius $r\sim 3.5$ Mpc from the MW-M31 midpoint.

We shall hereafter refer to luminous galaxies within the virial radius of
either primary as “satellites”. To minimize numerical resolution effects, we
shall only use for our analysis galaxies with at least $\sim 10$ star
particles, or stellar masses $M_{*}>10^{5}\,M_{\odot}$.

###  2.2 Observational data

In this work we consider all currently known dwarf galaxies within $3$ Mpc of
the midpoint between the MW and M31. We use position (RA, dec), distance
modulus $(m-M)$ and line-of-sight velocity data from the latest update of
McConnachie (2012)’s Nearby Dwarf Galaxy Database333See https://www.cadc-
ccda.hia-iha.nrc-cnrc.gc.ca/en/community/nearby/, and references therein..

We consider systems within 300 kpc of the MW or M31 as “satellites” of that
primary, and those further away as “field” dwarfs. Our total sample consists
of $142$ dwarfs of which $48$ are field galaxies. To be consistent with the
simulation limitations, we also impose a minimum stellar mass of
$M_{*}=10^{5}\,M_{\odot}$. This cut removes from our sample only MW
satellites, which are not the main focus of our study, leaving a total of $96$
dwarfs.

From the catalogued data, we have computed Galactocentric positions and radial
velocities assuming a Galactocentric distance for the Sun of $R_{\odot}=8.29$
kpc, a circular velocity for the local standard of rest (LSR) of $V_{0}=239$
km/s (McMillan, 2011), and a peculiar velocity with respect to the LSR of
$(U_{\odot},V_{\odot},W_{\odot})=(11.1,12.24,7.25)$ km/s (Schönrich et al.,
2010).

Figure 2: Distribution of the (cosine of) angular distances between dwarf
galaxies and the direction to the M31 analog, in concentric spherical shells
with Galactocentric distances as indicated in the upper legend of each panel.
Light green histograms show all galaxies (field + satellites) while teal
histograms show only field galaxies. The number of dwarfs considered in each
spherical shell is quoted in the panels. Top: APOSTLE dwarf galaxies. We
indicate average number of galaxies, i.e., total normalized by 8, as we
consider 8 different LG configurations alternating the MW and M31
identification (see text). The black histogram indicates the result of using,
on average, only as many systems as are actually observed (bottom panels) in
each spherical shell. Bottom: Observational data. The first shell includes
only MW satellites, whereas the third shell includes mainly M31 satellites.
The simulations show an excess of objects along the MW-M31 direction (at
$\cos(\alpha)\approx 1$ and $-1$); a trend that is not readily seen in the
observational data.

Figure 3: Aitoff projections of nearby dwarf galaxies in Galactic coordinates.
Top: One eighth of randomly selected APOSTLE galaxies are shown after stacking
8 different configurations, rotated so in each the direction to M31 matches
the observed position of M31 in the sky (see text). Bottom: Observed LG and
periphery. Satellite galaxies are shown with star symbols (gray for MW
satellites and blue for M31 satellites). Field dwarfs are shown as circles
colored according to $\alpha$. Filled circles correspond to LGP dwarfs and
open circles to the rest of field dwarfs at distances $<1.25$ Mpc. Black-
dashed circles mark an area of $45^{\circ}$ around the MW-M31 direction. In
the bottom panel, patches of different colors illustrate the footprints of the
observational surveys indicated in the legend. The effect of Galactic disk
obscuration is shown in yellow, and the Supergalactic Plane is marked with a
brown line.  Figure 4: Radial velocity versus distance for nearby dwarf
galaxies in the APOSTLE simulations. We stack 8 different LG configurations
and plot only dwarfs with $M*>10^{5}$ M⊙ within 3.5 Mpc of the LG midpoint.
Top: Centered on each of the 8 APOSTLE LG primaries. Bottom: Centered on the
barycenter of each LG. Satellites are shown with star symbols while field
galaxies are shown as circles. All dwarfs are colored according to $\alpha$ as
in previous figures (see legend). Black circles indicate the mean radial
velocity computed in radial bins within $1.25<d_{\rm MW}$/Mpc$<3$, and the
gray shaded areas show the $\pm 1\sigma$ deviation from the mean. The gray
dotted line shows a pure Hubble law with $H_{0}=73$ km/s/Mpc (Riess et al.,
2022). Tab. 1 gives linear fits to the recession velocity of LGP APOSTLE
dwarfs in the MW-centric reference frame (upper panel), expressed in the form
$V_{\rm rad}=\mathcal{H}(d_{\rm MW}/$Mpc$-1.25)+V_{1.25}$, for each angular
bin in $\alpha$.

##  3 Results

###  3.1 The spatial distribution of Local Group dwarf galaxies

The left panel of Fig. 1 shows the 3D-positions of APOSTLE nearby dwarfs
(field and satellite galaxies) with stellar masses $M_{*}>10^{5}$ M⊙. The
coordinate system is centered on the smaller of the two primaries in each
volume (we shall refer to it as the “MW analog”444In APOSTLE the less massive
of the two halo primaries is usually considered the Milky Way, although for
some of the analysis we shall drop the distinction between “MW” and “M31”
analogs.), with the “$X$” axis being coincident with the direction connecting
the two primaries. (M31 is thus located on the $X$-axis at roughly $X\sim 800$
kpc.) The plot includes all simulated dwarfs within a radius of $3.5$ Mpc from
the LG midpoint of each of the 4 APOSTLE-L1 volumes, stacked together (i.e.,
the plot shows four times as many dwarfs as would be expected in the true LG;
the increased number helps to visualize the 3D structure of the LG).

Galaxies are colored by the angle, $\alpha$, between the position of a galaxy
and the direction to the M31 analog; i.e., the angular distance between the
dwarf and the MW-M31 axis. Systems colored purple are those closer than
$45^{\circ}$ to the direction towards M31, while objects that are colored red
are those “behind” the MW along the same axis
($135^{\circ}<\alpha<180^{\circ}$). For reference, the gray spherical cones in
Fig. 1 highlight the $\alpha=45^{\circ}$ and $\alpha=135^{\circ}$ boundaries
in the volume.

The spatial distribution of APOSTLE nearby dwarfs is not isotropic, but
clearly elongated, with the majority of objects located close to the $X$-axis
joining the MW and M31 analogs. This seems true at all radii, not only for
relatively nearby galaxies, which include the satellites of the main halos,
but also for the more distant LGP dwarfs far outside their virial radii.

The anisotropic spatial distribution of APOSTLE nearby dwarfs is quantified in
Fig. 2, where we plot the angular distribution of galaxies on the sky in four
different spherical shells centred on the MW analog. Each panel shows a
histogram of $\cos(\alpha)$ for all $M_{*}>10^{5}\,M_{\odot}$ galaxies in a
shell.

The top row in Fig. 2 corresponds to the average of $8$ APOSTLE
configurations, two for each APOSTLE volume, alternating the designation of MW
or M31 analogs. This helps to reduce noise and to characterize more robustly
the anisotropic distribution in each shell.

Light-green histograms in Fig. 2 show results for all galaxies (i.e.,
satellites + field), while teal histograms show results considering only field
dwarfs. The first radial bin (shown in the leftmost panels) is, as expected,
dominated by the satellite population of each primary chosen as centre.

The second radial bin ($300<d_{\rm MW}$/kpc$<600$) is chosen to exclude most
satellites, whereas the third bin includes the majority of the satellites of
the second primary. Note as well that the second bin includes a high fraction
of so-called “backsplash” dwarfs ($>50\%$ within $2\times r_{200}$), i.e.,
galaxies found presently outside the virial radii of a primary but which were
in the past inside $r_{200}$ (Santos-Santos et al., 2023). The final,
outermost radial bin (rightmost panels in Fig. 2, $1.25<d_{\rm MW}$/Mpc$<3$)
includes only LGP field dwarfs.

A spatially uniform galaxy distribution would show in Fig. 2 as a constant
number of systems as a function of $\cos(\alpha)$. This is clearly not the
case for APOSTLE dwarfs, which show a clear excess at $\cos(\alpha)\approx 1$
and $\cos(\alpha)\approx-1$ for all radial shells. The “U-shaped”
distributions indicate a clear preference for dwarfs to align with the MW-M31
axis, a preference which persists even for LGP dwarfs located as far as $3$
Mpc away from either primary.

Is the same predicted anisotropy observed in our Local Group? The right-hand
panel of Fig. 1 shows the positions of observed dwarfs in a reference frame
with the Milky Way at the origin. A total of $96$ galaxies with
$M_{*}>10^{5}\,M_{\odot}$ are known within $3$ Mpc from the Milky Way,
whereas, on average, each of the APOSTLE volumes has $176$ such galaxies
within the same volume. In other words, APOSTLE predicts that roughly $80$
galaxies as massive as $M_{*}>10^{5}\,M_{\odot}$ may be currently missing from
our inventory of nearby dwarfs (see also Fattahi et al., 2020). Many of these
missing dwarfs are predicted to be in the outermost radial shell presented in
Fig. 2. In this figure, the number of such galaxies observed in each shell is
listed in each panel of the bottom row of that figure, and may be compared
with the average number of APOSTLE dwarfs quoted in the top row.

Intriguingly, the observed dwarfs in the LG and periphery do not seem to
follow the same anisotropic distribution predicted by APOSTLE. Indeed, there
is little evidence for a preferred alignment of dwarfs with the axis defined
by M31 and the MW; if anything, the opposite trend appears to prevail, with
clear peaks at $\cos(\alpha)\approx 0$, both in the satellites of the Milky
Way (bottom left panel of Fig. 2) and in the LGP dwarfs (bottom-right panel of
Fig. 2).

Could the difference be due to small-number statistics, or caused by
incompleteness in our inventory of currently known nearby dwarfs? The former
seems unlikely, as shown by the solid black histograms in the upper panels of
Fig. 2, which indicate the result of choosing at random in each APOSTLE radial
bin only as many galaxies as are available in the observational sample. The
“U-shaped” anisotropy is still clearly noticeable for APOSTLE galaxies but not
in the observations. A K-S test rejects the hypothesis that the distributions
in the upper and lower panels of Fig. 2 are drawn from the same parent
distribution at high confidence in all cases where numbers permit (see
p-values in the legends).

If the lack of “U-shaped” anisotropy in the observed dwarfs is caused by
incompleteness, then there should be a number of undetected dwarfs both behind
M31 and in the anti-M31 direction. We show this in Fig. 3, where the colored
circles in the upper Aitoff diagram shows the on-sky distribution of APOSTLE
field dwarfs randomly sampled from 8 possible LG configurations, obtained by
alternating the MW and M31 analogs. Filled circles represent LGP field dwarfs
while open circles indicate the rest of field dwarfs at distances $<1.25$ Mpc
from the primary. The satellites of each primary are shown by stars but are
de-emphasized in this plot.

The high density of red and purple LG dwarfs is quite clear in APOSTLE,
compared to that of galaxies further away from the MW-M31 direction: indeed,
while the circles around M31 and the anti-M31 direction cover only $30\%$ of
the sky, they contain roughly $57\%$ of all dwarfs within $3.5$ Mpc from the
LG barycenter (or $54\%$ of all dwarfs in the LGP). In terms of numbers,
according to APOSTLE we would expect to find, on average, $8$ LG and $13$ LGP
dwarfs in the red circle, and $14$ LG and $25$ LGP dwarfs in the purple one,
and a total of about $12$ LG and $33$ LGP dwarfs in other regions of the sky
(note this count is for dwarfs within $3.5$ Mpc from the LG barycentre).

In contrast, there are very few observed LGP dwarfs in the red and purple
areas, as shown in the bottom panel of Fig. 3: only $2$ in the M31 direction,
and $4$ in the anti-M31 direction. One reason for this may be that imaging
surveys have covered the sky unevenly, reducing the possibility of detecting
dwarfs in the M31 and anti-M31 directions. The “zone of avoidance” (i.e.,
$b<15^{\circ}$) caused by Galactic disk obscuration may also play a role:
given M31’s low Galactic latitude, roughly $\sim 39\%$ of the sky with
$\alpha<45^{\circ}$ or $\alpha>135^{\circ}$ is obscured by the disk.

The bottom panel of Fig. 3 also suggests that the anti-M31 direction, in
particular, might not have been surveyed as deeply as other parts of the sky.
Indeed, most LGP dwarfs have been found by visual inspection of photographic
plates (see; e.g., Whiting et al., 1999; Karachentseva & Karachentsev, 1998),
which suggests that there is scope for new discoveries with the advent of
digital surveys of the whole sky. The anti-M31 direction is also outside the
footprint of both the Sloan Digital Sky Survey, which covers most of the
northern sky (Willman et al., 2005, SDSS; gray-shaded area), and of the DES
survey (purple contour, Drlica-Wagner et al., 2015), which has imaged a large
fraction of the southern sky. There have been several surveys in the general
direction of M31 (most notably the PAndAS survey, shown with a cyan contour,
McConnachie et al., 2009) but they only cover a very small fraction of the
$\alpha<45^{\circ}$ area of the sky around M31.

Recently, McNanna et al. (2023) have searched for field dwarfs in the DES
footprint with distances between $300$ kpc and $2$ Mpc from the MW. Although
their search should detect all $M_{*}>10^{5}\,M_{\odot}$ galaxies like the
ones we study here, they report no new discoveries aside from the $7$ already
known dwarfs in that volume. For comparison, averaging all volumes and various
orientations, APOSTLE predicts a median number of $3$ in that region, with a
$10$th and $90$th percentile of $0$ and $8$, respectively (and a full range
that goes from $0$ to $37$ dwarfs). Given the large volume-to-volume scatter
in APOSTLE, we find no obvious conflict between the result of the latest
search in the DES footprint and the results from the simulations.

To summarize, the APOSTLE runs predict a clear anisotropic distribution for
nearby galaxies, which are expected to align strongly with the M31 and
anti-M31 directions. This prediction seems at odds with currently available
data on LG and LGP dwarfs. The reason for the disagreement is unclear, but, if
due to incompleteness, then the simulations suggest that those directions
could prove fruitful targets for future searches of LGP dwarfs.

Table 1: Parameters of linear fits to the recession velocity of LGP dwarfs in APOSTLE (in the MW-centric reference frame) of the form $V_{\rm rad}=\mathcal{H}(d_{\rm MW}/$Mpc$-1.25)+V_{1.25}$, for each of the angular bins $\alpha$. These are illustrated as colored lines in the inset to Fig. 5. APOSTLE |  $\mathcal{H}$ [km/s] |  $V_{1.25}$ [km/s]
---|---|---
Mean MW-centric | 134.3 | -17.3
$\alpha<45^{\circ}$ | 120 | -61.2
$45^{\circ}>\alpha<90^{\circ}$ | 102.4 | -13.4
$90^{\circ}>\alpha<135^{\circ}$ | 104.6 | 26.5
$\alpha>135^{\circ}$ | 97.3 | 86.0
Mean LG-centric | 108.3 | 7.8

###  3.2 The Local Hubble Flow in APOSTLE

The anisotropic distribution of mass predicted by APOSTLE should also have
consequences on the velocity field of galaxies around the Milky Way. As
discussed in Sec. 1, galaxies beyond the LG turnaround radius should be
expanding away with a decelerated Hubble flow. The velocity field is also
expected to be fairly symmetric relative to the LG barycentre, and to show a
clear asymmetry when expressed in the Galactocentric (MW) reference frame.
From that perspective, distant galaxies at given distance from the MW, $d_{\rm
MW}$, should have recession velocities which depend on the angular distance to
the direction of M31, peaking in the anti-M31 direction and having a minimum
behind M31.

We show this in the upper panel of Fig. 4, which shows the radial velocity
relative to the MW analog in APOSTLE, as a function of MW-centric distance.
Satellite galaxies are shown with star symbols, while field galaxies are shown
as circles, all colored according to $\alpha$ as in previous figures. Black
squares indicate the radial distance and velocity of the “M31 analog” in each
volume, and thin gray lines indicate the virial radii of all primaries, for
reference. Results are stacked for 8 LG configurations, obtained by
alternating the MW and M31 analogs.

Galaxies in Fig. 4 are, as in previous plots, colored by their angular
distance to the M31 analog. In particular, systems in purple are in the
general direction of M31, and those in red are located close to the anti-M31
direction. The black points and gray shade represent the mean and $\pm
1\sigma$ standard deviation of $V_{\rm rad}$ computed in radial bins in the
range $1.25$ to $3$ Mpc. The average rms about the mean velocity is quoted in
each panel.

As may be seen from the top panel of Fig. 4, the recession velocities of LGP
dwarfs, as seen from the MW, are clearly modulated by their angular distance
to M31. The difference is not subtle. Behind M31, galaxies at $d_{\rm MW}\sim
2$ Mpc have a mean recession velocity of just $\sim 30$ km/s. At the same
distance, galaxies in the anti-M31 direction are receding on average at
roughly $\sim 160$ km/s.

This angular dependence hinders a proper characterization of the local
recession velocity field, including a precise determination of the LG
turnaround radius (where the mean radial velocity vanishes), the dispersion
about the mean flow (the “coldness” of the local Hubble flow), and how
decelerated the local velocity field is relative to the pure Hubble flow.

As expected, the angular dependence disappears when referring velocities and
distances to the barycentre of the M31-MW system, as shown in the bottom panel
of Fig. 4. This makes it easy to estimate the average APOSTLE LG turnaround
radius ($r_{\rm ta}\sim 1.2$ Mpc from the LG barycentre), and the velocity
dispersion about the mean flow ($\sim 40$ km/s). The flow is clearly
decelerated relative to a pure Hubble flow (indicated by the gray dotted line
in the bottom panel of Fig. 4): even galaxies as far away as $\sim 3$ Mpc from
the LG barycentre have not yet reached the pure Hubble flow.

Figure 5: Galactocentric radial velocity versus radial distance for observed
dwarf galaxies. Symbols and color-coding are the same as in Fig. 4. A thin
gray line and shaded area indicate the linear fit and standard deviation
obtained for galaxies in the range $1.25<d_{\rm MW}$/Mpc$<3$. For comparison,
we overplot the black circles shown in the top panel of Fig. 4, i.e., the mean
distance-velocity relation for APOSTLE. The dotted line marks a pure Hubble
law with $H_{0}=73$ km/s/Mpc. The inset panel shows the same observational
data with the addition of 4 lines representing linear fits to APOSTLE LGP
dwarfs in each of the four $\alpha$ bins (see Tab. 1).

How do these results compare with observational data for the Local Group? In
this case, we can only compute the Galactocentric flow, because proper motions
for all galaxies, which are unavailable, or a detailed 3D velocity model,
would be needed to refer recession velocities to the LG barycentre (see; e.g.,
Karachentsev et al., 2009; Peñarrubia et al., 2014).

Fig. 5 shows the distance-recession velocity relation for known dwarfs,
measured with respect to the MW. The thin solid gray line indicates a linear
fit to the data in the distance range $1.25<d_{\rm MW}$/Mpc$<3$. The shaded
area around that fit shows the corresponding $1\sigma$ standard deviation, of
order $\sim 65$ km/s. Note that this is in excellent agreement with the
coldness of the Galactocentric local Hubble flow in APOSTLE (see upper panel
of Fig. 4), which implies that the APOSTLE runs have no problem accounting for
the observed coldness of the local Hubble flow.

There is also good qualitative agreement between APOSTLE and observations
regarding the angular dependence of the recession velocity. As visual
inspection shows, at given distance, the recession velocities are highest in
the anti-M31 direction and lowest behind M31, with a velocity difference
between antipodal directions exceeding $100$ km/s. We note that the reflex
motion of the Milky Way caused by the recent infall of the Magellanic Clouds
system could add a similar dipole-like effect along a similar axis, but the
amplitude in that case is likely to be much smaller, roughly of only $\sim 30$
km/s according to Petersen & Peñarrubia (2020). The velocity anisotropy in the
MW frame is therefore mainly due to the offset between the Milky Way and the
LG barycentre, as shown for APOSTLE in Fig. 4.

The main difference between observations and simulations is that the observed
local Hubble flow seems much less decelerated than that of APOSTLE. At given
distance form the MW, APOSTLE’s Galactocentric recession velocities (shown by
the black connected circles in Fig. 5) are well below observed ones (shown by
the solid gray line; the dotted gray line indicates a pure Hubble flow with
$H_{0}=73$ km/s/Mpc, Riess et al. 2022).

One way of reconciling this difference would be to assume that the combined
MW+M31 mass is much lower than assumed in APOSTLE, but this seems unlikely.
Indeed, as discussed by Fattahi et al. (2016), a total combined mass of order
$10^{11.5}$ M⊙ would be needed, at least an order of magnitude lower than
current estimates.

Could the disagreement be due instead to the patchy sky coverage of the
observational sample? We explore this in the inset to Fig. 5, where we compare
linear fits to the MW-centric APOSTLE Hubble flow on different parts of the
sky (colored lines) with the observational data. The disagreement is worse for
galaxies with $\alpha>90^{\circ}$ (i.e., the hemisphere in the anti-M31
direction, shown in red and green), where the observed recession speeds
systematically exceed the APOSTLE predictions. Although the disagreement seems
clear, it is important to keep in mind how sparse the sky coverage of LGP
dwarfs is: only $4$ “red” distant dwarfs are known in the anti-M31 direction,
and only $2$ (shown in purple) are known behind M31.

It is therefore certainly possible that the comparison may change if even a
few more systems are added in each of these directions. APOSTLE predicts that
about $\sim 60$ field dwarfs with $M_{*}>10^{5}\,M_{\odot}$ are likely be
missing from our nearby dwarf galaxy inventory, $41$ of them LGP dwarfs (see
Sec. 3.1). In particular, the lack of deceleration in the observed MW-centric
Hubble flow could be due to the paucity of LGP dwarfs behind M31 (i.e., purple
circles in Fig. 5). Finding some of the “missing dwarfs” predicted by APOSTLE
in the MW-M31 direction could certainly impact our understanding of the Hubble
flow in the outskirts of the Local Group.

If not due to incompleteness, it is possible that the apparent lack of
deceleration in the velocity field of LGP dwarfs may be due to the large-scale
distribution of matter around the Local Group. This could in principle be
investigated using simulations tailored to reproduce not only the LG
environment but also that of its surrounding large-scale structure (see; e.g.,
Carlesi et al., 2017; Libeskind et al., 2020; Sawala et al., 2022).

##  4 Summary and Conclusions

We have used the APOSTLE suite of cosmological hydrodynamical simulations to
study anisotropies in the spatial distribution and kinematics of dwarf
galaxies in the Local Group and periphery, out to a distance of $\sim 3$ Mpc
from the MW. The simulations show that the anisotropy induced by the presence
of two massive primaries on first approach is reflected in the spatial
distribution of nearby dwarf galaxies. At all distances from the MW, the
simulations predict a strong preference for dwarfs to be located close to the
axis defined by the MW-M31 direction, from the satellites of the primary
galaxy to even the “distant” LGP dwarfs, defined as those in the LG periphery,
i.e., at distances $1.25<d_{\rm MW}/$Mpc$<3$ from the MW.

The local “Hubble flow” of the LGP dwarfs is also expected to be anisotropic
if measured in the Galactocentric rest frame. At fixed distance from the MW
the mean recession speed, $\langle V_{\rm rad}\rangle$, varies with angular
distance to M31, peaking in the anti-M31 direction and reaching a minimum
behind M31, mainly as a result of the offset between the MW and the LG
barycentre, which lies somewhere between MW and M31.

The combined M31-MW mass also decelerates the local Hubble flow of LGP dwarfs;
the LG “turnaround radius” (i.e., where $\langle V_{\rm rad}\rangle=0$) in
APOSTLE is located at $r_{\rm ta}\sim 1.2$ Mpc from the LG barycentre and the
pure Hubble flow (i.e., $\langle V_{\rm rad}\rangle=H_{0}*r$) is not reached
out to at least $r\sim 3$ Mpc. The predicted flow is very cold, with a
barycentric dispersion of less than $\sim 40$ km/s.

A comparison of these features with existing observations raises interesting
questions. Although there is agreement with the predicted angular anisotropy
in recession velocities around the Milky Way, there is little evidence in the
spatial distribution of LGP dwarfs for a preferred direction along the MW-M31
direction.

The “coldness” of the local Hubble flow also seems consistent with the
simulations, but it is significantly less decelerated. Indeed, in the
Galactocentric frame, all dwarfs beyond $r\sim 1.25$ Mpc seem to be receding
with velocities consistent with a pure, undecelerated Hubble flow. Although
the reason for these differences is so far unclear, APOSTLE also predicts that
the true number of LGP dwarfs should be substantially higher than observed,
suggesting that our local inventory of dwarfs is rather incomplete.

Another reason for the disagreement may be that the APOSTLE volume selection
made no attempt to account for structures beyond $\sim 3$ Mpc from the LG
barycentre. The presence of large galaxies just outside that volume, like M81
or NGC 5128, as well as the influence of the Virgo cluster or the Local Void,
may all have an influence over the spatial distribution and velocity field of
LGP dwarfs and clearly need to be taken into account in future, higher
fidelity simulations of the LG volume.

It is thus possible that the oddities described above may result at least in
part from incompleteness and inhomogeneous sky coverage, but a full
explanation will need to await the completion of deep all-sky surveys able to
fill the gaps in our current inventory of the Local Group, and of simulations
able to fully reproduce the configuration of the Local Group within the
larger-scale distribution of matter in the local Universe.

## Acknowledgements

ISS acknowledges support from the European Research Council (ERC) through
Advanced Investigator grant to C.S. Frenk, DMIDAS (GA 786910). We wish to
acknowledge the generous contributions of all those who made possible the
Virgo Consortium’s EAGLE/APOSTLE simulation projects. This work used the
DiRAC@Durham facility managed by the Institute for Computational Cosmology on
behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The equipment was
funded by BEIS capital funding via STFC capital grants ST/K00042X/1,
ST/P002293/1, ST/R002371/1 and ST/S002502/1, Durham University and STFC
operations grant ST/R000832/1. DiRAC is part of the National e-
Infrastructure. This research made use of Astropy (http://www.astropy.org) a
community-developed core Python package for Astronomy.

## Data Availability

The simulation data underlying this article can be shared on reasonable
request to the corresponding author. The references for the observational data
for nearby dwarfs used in this article are listed in Sec. 2.2.

## References

  * Bechtol et al. (2015) Bechtol K., et al., 2015, ApJ, 807, 50
  * Benitez-Llambay & Frenk (2020) Benitez-Llambay A., Frenk C., 2020, MNRAS, 498, 4887
  * Carlesi et al. (2017) Carlesi E., Hoffman Y., Sorce J. G., Gottlöber S., 2017, MNRAS, 465, 4886
  * Cautun et al. (2020) Cautun M., et al., 2020, MNRAS, 494, 4291
  * Crain et al. (2015) Crain R. A., et al., 2015, MNRAS, 450, 1937
  * Davis et al. (1985) Davis M., Efstathiou G., Frenk C. S., White S. D. M., 1985, ApJ, 292, 371
  * Drlica-Wagner et al. (2015) Drlica-Wagner A., et al., 2015, ApJ, 813, 109
  * Fattahi et al. (2016) Fattahi A., et al., 2016, MNRAS, 457, 844
  * Fattahi et al. (2020) Fattahi A., Navarro J. F., Frenk C. S., 2020, MNRAS, 493, 2596
  * Kahn & Woltjer (1959) Kahn F. D., Woltjer L., 1959, ApJ, 130, 705
  * Karachentsev et al. (2009) Karachentsev I. D., Kashibadze O. G., Makarov D. I., Tully R. B., 2009, MNRAS, 393, 1265
  * Karachentseva & Karachentsev (1998) Karachentseva V. E., Karachentsev I. D., 1998, A&AS, 127, 409
  * Komatsu et al. (2011) Komatsu E., et al., 2011, ApJS, 192, 18
  * Li & White (2008) Li Y.-S., White S. D. M., 2008, MNRAS, 384, 1459
  * Libeskind et al. (2020) Libeskind N. I., et al., 2020, MNRAS, 498, 2968
  * Macciò et al. (2005) Macciò A. V., Governato F., Horellou C., 2005, MNRAS, 359, 941
  * McConnachie (2012) McConnachie A. W., 2012, AJ, 144, 4
  * McConnachie et al. (2009) McConnachie A. W., et al., 2009, Nature, 461, 66
  * McMillan (2011) McMillan P. J., 2011, MNRAS, 414, 2446
  * McNanna et al. (2023) McNanna M., et al., 2023, arXiv e-prints, p. arXiv:2309.04467
  * Patel & Mandel (2023) Patel E., Mandel K. S., 2023, ApJ, 948, 104
  * Peñarrubia et al. (2014) Peñarrubia J., Ma Y.-Z., Walker M. G., McConnachie A., 2014, MNRAS, 443, 2204
  * Pereira-Wilson et al. (2023) Pereira-Wilson M., Navarro J. F., Benítez-Llambay A., Santos-Santos I., 2023, MNRAS, 519, 1425
  * Petersen & Peñarrubia (2020) Petersen M. S., Peñarrubia J., 2020, MNRAS, 494, L11
  * Riess et al. (2022) Riess A. G., et al., 2022, ApJ, 934, L7
  * Sales et al. (2011) Sales L. V., Navarro J. F., Cooper A. P., White S. D. M., Frenk C. S., Helmi A., 2011, MNRAS, 418, 648
  * Sandage et al. (1972) Sandage A., Tammann G. A., Hardy E., 1972, ApJ, 172, 253
  * Santos-Santos et al. (2023) Santos-Santos I. M. E., Navarro J. F., McConnachie A., 2023, MNRAS, 520, 55
  * Sawala et al. (2016) Sawala T., et al., 2016, MNRAS, 457, 1931
  * Sawala et al. (2022) Sawala T., McAlpine S., Jasche J., Lavaux G., Jenkins A., Johansson P. H., Frenk C. S., 2022, MNRAS, 509, 1432
  * Schaye et al. (2015) Schaye J., et al., 2015, MNRAS, 446, 521
  * Schlegel et al. (1994) Schlegel D., Davis M., Summers F., Holtzman J. A., 1994, ApJ, 427, 527
  * Schönrich et al. (2010) Schönrich R., Binney J., Dehnen W., 2010, MNRAS, 403, 1829
  * Springel et al. (2001) Springel V., Yoshida N., White S. D. M., 2001, New A, 6, 79
  * Whiting et al. (1999) Whiting A. B., Hau G. K. T., Irwin M., 1999, AJ, 118, 2767
  * Willman et al. (2005) Willman B., et al., 2005, ApJ, 626, L85


# Unsupervised physics-informed disentanglement of multimodal data
for high-throughput scientific discovery

Nathaniel Trask  Carianne Martinez  Kookjin Lee  Brad Boyce

###### Abstract

We introduce physics-informed multimodal autoencoders (PIMA) - a variational
inference framework for discovering shared information in multimodal
scientific datasets representative of high-throughput testing. Individual
modalities are embedded into a shared latent space and fused through a product
of experts formulation, enabling a Gaussian mixture prior to identify shared
features. Sampling from clusters allows cross-modal generative modeling, with
a mixture of expert decoder imposing inductive biases encoding prior
scientific knowledge and imparting structured disentanglement of the latent
space. This approach enables discovery of fingerprints which may be detected
in high-dimensional heterogeneous datasets, avoiding traditional bottlenecks
related to high-fidelity measurement and characterization. Motivated by
accelerated co-design and optimization of materials manufacturing processes, a
dataset of lattice metamaterials from metal additive manufacturing
demonstrates accurate cross modal inference between images of mesoscale
topology and mechanical stress-strain response.

Multimodal machine learning, physics-informed machine learning, variational
inference, variational autoencoders, fingerprinting, mixture of experts

##  1 Motivation

Many scientific and engineering datasets are multimodal, necessitating the
fusion of disparate sources and datatypes for informed analysis. For example,
in the realm of process optimization for materials manufacturing, processes
ranging from microelectronic fabrication to metal additive manufacturing
involve a myriad of process settings along with in-process and post-process
measurements (Liu et al., 2012; Sochol et al., 2018). Moreover, automated
high-throughput characterization methods are increasingly generating large,
rich, multimodal datasets, fueled by advances in robotics and automation
(Boyce & Uchic, 2019). Many scientific datasets admit fingerprints: easily
measurable signals which correlate with a difficult to measure underlying
physical process. The hunt for exploitable fingerprints extends beyond
material science (Isayev et al., 2015), and can be found in all science and
engineering domains ranging from quantum mechanics (Chakraborty et al., 2021)
to climate change (Hasselmann, 1997; Hegerl et al., 2007). Rapid datasets
designed to detect fingerprints may potentially serve as a surrogate for, or
in conjunction with, bespoke experiments capturing high-fidelity modalities.
Accordingly, we aim to discover comprehensive fingerprints constructed from
the weighted integration of several disparate data sources, each with unique
fidelity, sparsity, and spatiotemporal resolution. The multimodal scientific
data under consideration includes both physical and simulated data and differs
from the text/audio/video modalities commonly considered in the multimodal
literature (Baltrušaitis et al., 2018), providing the opportunity to impose
physics-based inductive biases and move beyond purely data-driven linear
techniques such as principle component analysis typically used for fingerprint
detection.

Herein, we present a novel variational inference framework for synthesizing
multimodal scientific data with the aim of cross-modal inference; if one can
reliably perform generative modeling of a high-fidelity but slow measurement
from a low-fidelity but fast fingerprint, high-throughput experimentation and
material characterization are no longer infeasible. Such applications however
mandate an unsupervised approach, as costly human-in-the-loop data labelling
precludes high-throughput testing.

Concretely, cross-modal inference corresponds to training jointly across
modalities $X_{1},...,X_{D}$ in a manner that supports generative sampling of
individual modalities $p(X_{i}|X_{j})$ for $i\neq j$. We achieve this in a
variational inference setting by combining the following algorithmic
contributions (Figure 1): 1. encoding data into unimodal embeddings
$q(Z|X_{i})$ and applying a product of experts (PoE) model to fuse data into a
multimodal posterior $q(Z|\mathbf{X})=\Pi_{i}q(Z|X_{i})$; 2. adopting a
Gaussian mixture prior $p(Z|c)$ to determine latent clusters shared across
modalities; and 3. decoding with a physics-informed mixture of experts (MoE)
model $p(X_{i}|c,Z)$ to impose inductive biases. For scientific settings, the
expert model provides a critical new means of fusing experimental audiovisual
data with traditional scientific models; rather than considering generalized
linear models commonly used in MoE (Jordan & Jacobs, 1994), we may incorporate
parameterized physical models, surrogates or simulators for the system under
consideration. These ingredients are designed to yield an ELBO loss with
closed form expressions for requisite integrals and is amenable to a novel
expectation maximization strategy to fit clusters and experts. In concert,
this architecture produces fingerprints in the form of latent clusters
spanning modalities, with cross-modal estimators allowing inference of cluster
membership for a single modality.

Figure 1: Individual modalities are encoded into Gaussian distributions in a
shared latent space. During training the posterior is sampled from a product
of experts distribution fusing complementary information into a shared
multimodal Gaussian distribution. A Gaussian mixture prior parameterizes
clusters encoding cross-modal shared information. Sampling from mixture
components provides generative models using either black-box decoders or
expert physics models encoding prior physics knowledge. To facilitate cross-
modal generative inference, unimodal embeddings are trained to reproduce the
multimodal embedding, allowing inference of $p(c|X_{i})$. Shown here, an
expert strain-hardening plasticity model allows two types of cross-modal
inference: costly measurements of stress-strain response may be inferred from
high-throughput imaging of lattice topology, or generative microstructural
images can be provided to suggest microstructure which correlate with a given
stress/strain measurement.

For unsupervised learning, several works apply variational autoencoders (VAE)
to seek latent disentangled representations of data which admit efficient
separation into meaningful classes (Burgess et al., 2018; Chen et al., 2018;
Locatello et al., 2019; Kim & Mnih, 2018). While desirable from an
interpretability and accuracy perspective, such representations are often
challenging to reliably discover in the absence of labels. The complementary
information available in multimodal data has been shown to provide multiple
pathways to disentanglement; a human may be unable to distinguish an image of
a one and a seven, but if the digit is read aloud there is no confusion.
However, the fact that scientific data is governed by physical models
potentially allows the expert model to extract more information than purely
data-driven approaches - known physics encodes the generative process and
therefore imposing even a low-fidelity physical model as inductive bias may
provide substantial disentanglement. Our tests demonstrate that the
combination of PoE, GMM, and expert models provide not only disentangled
clusters, but also an ordering of clusters in latent space reflecting
information shared across modalities. For examples with unambiguous class
ownership, we demonstrate 100% test accuracy when classifying data into
clusters. We anticipate these physics-based disentangled representations will
enable future causal analysis of large high-throughput datasets.

###  1.1 Relationship to prior literature

This work draws from several thematic bodies of literature. The non-exhaustive
list below denotes those works which have most informed our approach as well
as provide overviews for the recent state-of-the-art.

#### Gaussian mixture embeddings

For deep unsupervised clustering, several works replace the standard normal
prior from (Kingma & Welling, 2014; Rezende et al., 2014) with a Gaussian
mixture model (GMM) to facilitate disentanglement and provide an explicit
parameterization of clusters (Dilokthanakul et al., 2016; Jiang et al., 2017;
Rao et al., 2019; Lee et al., 2020). Each modality in the multi-modal prior
distribution is expected to provide disentangled latent representations of
data which admit an explicit parameterization of class distributions. The
current work is most similar to VaDE (Jiang et al., 2017) in its use of mean-
field distributions to obtain a separable ELBO, and Bayesian estimator for
$q(c|X)$. This work builds upon VaDE by incorporating multimodal data inputs
while maintaining computational tractability of the ELBO, as well as employing
clusters to decode into physics-informed MoE models.

#### Disentanglement

Another line of research is to extract latent disentangled representations
into different factors of variations in data using VAEs. Earlier works such as
$\beta$-VAE (Higgins et al., 2017) and Annealed VAE (Burgess et al., 2018)
introduce additional weighting parameters to the KL divergence term of the
original VAE ELBO loss. In Factor VAE (Kim & Mnih, 2018) and $\beta$-TCVAE
(Chen et al., 2018) the ELBO is further decomposed to derive and penalize the
total correlation to promote disentanglement in learned representations. For
our purposes however, without an explicit parameterization of the cluster
distributions to condition off of, it is not possible to introduce a physics-
informed expert MoE model.

#### Multimodal inference

Generative modeling from multimodal data can be broadly categorized into
either conditional generative models (Sohn et al., 2015; Pu et al., 2016)
which directly learn conditional cross-modal distributions $p(X_{i}|X_{j})$,
or joint models (Suzuki et al., 2017; Vedantam et al., 2018; Wu & Goodman,
2018), which explicitly learn joint distributions that learn
$p(Z,X_{1},...,X_{D})$. We pursue the later as (Wu & Goodman, 2018) has been
shown to provide better description of the underlying data distribution. We
pursue the strategy used by works such as joint multimodal VAE (Suzuki et al.,
2017) and joint VAE (Vedantam et al., 2018), where a a joint inference network
$q(Z|X_{1},X_{2})$ is trained, followed by training of two additional unimodal
inference networks $q(Z|X_{1})$ and $q(Z|X_{2})$ which handle missing data at
test time. The unimodal inference networks are trained to either match the
joint inference network or to maximize an ELBO derived to perform unimodal
variational inference. More recently, MVAE (Wu & Goodman, 2018) and MMVAE (Shi
et al., 2019) were proposed to model the joint posterior as a product of
experts (PoEs) and a mixture of experts (MoEs). Most recently, MoPoE-VAE
(Sutter et al., 2021) proposed a new ELBO formulation, which generalizes ELBO
formulations derived from PoEs and MoEs. Our encoder bears similarities to
MMVAE, MoPoE-VAE, and PoE, while preserving a computationally tractable closed
form ELBO when combined with the GMM prior.

#### Physics-informed ML and fingerprinting

Substantial works in recent years have focused on introducing physics into
either solving partial differential equations (PDEs) or for building
surrogates, typically introducing a PDE residual regularizer in physics-
informed neural networks (Lagaris et al., 1998; Raissi et al., 2019) or by
embedding physics directly into network architecture in structure-preserving
ML (Trask et al., 2020). Such tools can be combined to provide parametric
surrogates of simulations which can perform real-time inference over a
database of parameterized PDE solutions (Lu et al., 2019; Wang et al., 2021;
Mao et al., 2021). This paper provides a framework to fuse either these
physics-informed surrogates or simpler empirical models together with
experimental data. In contrast to traditional tools for fingerprinting which
rely on purely data-driven techniques like PCA (Hasselmann, 1997; Hegerl et
al., 2007), the current framework provides a means to incorporate domain
expertise into fingerprints tailored toward a scientific task.

#### Major contributions:

  * •

Novel fusion of PoE w/ Gaussian mixture to obtain parameterized cluster
“fingerprints” for downstream data analysis and high-throughput diagnostic
tasks.

  * •

Multimodal embedding allows cross-modal inference while preserving closed form
expressions for expectations in ELBO.

  * •

Mixture of experts decoding allows incorporation of interpretable inductive
biases by assuming model form describing scientific processes. Potential for
embedding physics-informed surrogates or simulators.

  * •

Improvements over SotA unimodal unsupervised techniques.

  * •

Disentanglement of clusters into structured latent space exposing
relationships across modalities.

##  2 Framework construction

Given individual modalities $X_{i}\subset\mathbb{R}^{d_{i}}$, we partition the
set of all modalities $\mathcal{M}=\left\\{X_{1},...,X_{D}\right\\}$ into
$\mathcal{M}_{DD}$ consisting of images/videos/audio amenable to purely data-
driven modeling, and $\mathcal{M}_{S}$ consisting of scientific modalities
amenable to expert modeling. For each $X_{i}\in\mathcal{M}$ we seek an
embedding $Z\in\mathbb{R}^{l}$ in latent dimension $l<<{d_{i}}$. Assuming a
categorical variable $c$ clustering data into $C$ clusters in latent space,
our variational autoencoder amounts to introducing parameterized prior $p$ and
posterior $q$ distributions that maximize the following ELBO loss:

$\mathcal{L}=\mathbb{E}_{q(Z,c|X_{1},...,X_{D})}\left[\log\frac{p(X_{1},...,X_{D},Z)}{q(Z,c|X_{1},...,X_{D})}\right].$
(1)

We further assume separability of both prior and posterior:

$\displaystyle
p(X_{1},...,X_{D},Z,c)=\left(\Pi_{i=1}^{D}p(X_{i}|Z,c)\right)p(Z|c)p(c),$ (2)
$\displaystyle
q(Z,c|X_{1},...,X_{D})=q(Z|X_{1},...,X_{D})q(c|X_{1},...,X_{D}).$ (3)

Our framework consists of four components: 1. unimodal deep encodings with a
product of experts (PoE) multimodal fusion, 2. a mixture of Gaussians prior,
3. a mixture of experts decoding of modalities $X\in\mathcal{M}_{S}$, and 4.
unimodal encoders for cross-modal inference. We introduce each component
sequentially, derive a closed form expression for the ELBO, and introduce an
expectation maximization assignment of clusters and expert models.

###  2.1 Multi-modal embedding

Assuming the unimodal embeddings may be modeled as multivariate Gaussians with
diagonal covariance, we obtain posterior probabilities
$q(Z_{i}|X_{i})=\mathcal{N}(Z_{i};\mu_{i},\sigma^{2}_{i}\mathbf{I})$, with
mean and covariance provided by the set of neural networks

$[\mu_{i},\sigma^{2}_{i}]=F_{i}(X_{i};\theta_{i}),$ (4)

where $\theta_{i}$ denotes trainable weights and biases. For this work we
consider a simple class of 1D/2D convolutional encoders, whose architecture is
provided in Appendix A.

To estimate $q(Z|X_{1},...,X_{D})$ in the ELBO, it follows from Bayes’ rule
and pairwise independence that

$q(Z|X_{1},...,X_{D})=q(Z)^{1-D}\Pi_{i=1}^{D}q(Z|X_{i}),$ (5)

so that the posterior is a scaled product of individual modalities. To obtain
closed form expressions for the ELBO later, we assume

$q(Z|X_{1},...,X_{D})\propto\Pi_{i=1}^{D}q(Z|X_{i}).$ (6)

The product of Gaussian distributions is again Gaussian, yielding the
multimodal distribution:

$\displaystyle q(Z|X_{1},...,X_{D})=\mathcal{N}(\mu,\sigma^{2}\mathbf{I}),$
(7)
$\displaystyle\sigma^{-2}=\sum_{i=1}^{D}\sigma_{i}^{-2},\qquad\frac{\mu}{\sigma^{2}}=\sum_{i=1}^{D}\frac{\mu_{i}}{\sigma_{i}^{2}},$
(8)

which may be sampled during training using the reparameterization trick: by
sampling $\epsilon\sim\mathcal{N}(0,\mathbf{I})$ and calculating
$Z=\mu+\epsilon\odot\sigma$, we may back-propagate through the random node $Z$
into the unimodal encoders, where $\odot$ denotes the Hadamard product.

###  2.2 Gaussian mixture prior and expert decoding

We adopt a simple Gaussian mixture prior, modeling

$\displaystyle p(c)=Cat(c|\mathbf{\pi}),$ (9) $\displaystyle
p(Z|c)=\mathcal{N}(\mu_{c},\sigma^{2}_{c}\mathbf{I}),$ (10) $\displaystyle
p(X_{i}|Z,c)=\mathcal{N}(\hat{\mu}_{i},\hat{\sigma}^{2}_{i}\mathbf{I}).$ (11)

To ensure a positive $\mathbf{\pi}$ that sums to unity, we parameterize it as
the softmax of a trainable vector. To decode $X_{i}\in\mathcal{M}_{DD}$, we
employ a neural network with parameters $\hat{\theta_{i}}$

$[\hat{\mu}_{i},\hat{\sigma}^{2}_{i}]=G_{i}(X_{i};\hat{\theta_{i}}).$ (12)

For $X_{i}\in\mathcal{M}_{S}$, we assume an expert model
$p(X_{i}|c)=\mathcal{N}(\mathcal{E}(t;\hat{\theta}_{c}),\hat{\sigma}^{2}_{c}\mathbf{I})$,
where $t$ is an independent variable and $\hat{\theta}_{c}$ denotes expert
parameters associated with each cluster. The specific choice of $\mathcal{E}$
will be problem dependent and specified in the experiment section.

Because $p(X_{i}|Z,c)$ admits interpretation as a mixture of experts model
(Jordan & Jacobs, 1994), we obtain closed form expressions for the mean and
variance to facilitate postprocessing and uncertainty quantification:

$\displaystyle\mathbf{E}[X_{i}]=\sum_{c}\pi_{c}\mathcal{E}(t;\hat{\theta}_{c}),$
(13)
$\displaystyle\text{Var}[X_{i}]=\left(\sum_{c}\pi_{c}(\hat{\sigma}_{c}^{2}-\mathcal{E}(t;\hat{\theta}_{c})^{2})\right)-\mathbf{E}[X_{i}]^{2}.$
(14)

In practice, $\mathcal{E}$ may take a variety of forms and its judicious
selection imparts significant prior knowledge. We consider in this work simple
generalized linear models and, for the mechanical data, an empirical linear
strain-hardening model. In general, these could range from empirical
engineering correlations obtained from e.g. dimensionless analysis or singular
perturbation theory, to analytic parametric solutions to PDE based models, or
to parametric physics-informed ML surrogates/reduce order models (see e.g. (Lu
et al., 2019; Wang et al., 2021; Mao et al., 2021; Trask et al., 2020)).

###  2.3 ELBO loss and EM minimizer

A modification of the derivation in (Jiang et al., 2017) to account for
multimodality yields the closed form for the single sample ELBO:

$\displaystyle\mathcal{L}_{d}=-\sum_{X_{i}\in\mathcal{M}_{DD}}||X_{i}-\hat{\mu}_{i}||^{2}$
$\displaystyle-\sum_{X_{i}\in\mathcal{M}_{S}}\sum_{t_{n}}\sum_{c}\gamma_{c}\left(\log\hat{\sigma}^{2}_{c}+\frac{(X_{i,n}-\mathcal{E}(t_{n};\hat{\theta}_{c}))^{2}}{\hat{\sigma_{c}}^{2}}\right)$
$\displaystyle-\sum_{c}\sum_{j}\gamma_{c}\left(\log\sigma_{c,j}^{2}+\frac{\sigma^{2}_{j}}{\sigma^{2}_{c,j}}+\frac{(\mu_{j}-\mu_{c,j})^{2}}{\sigma^{2}_{c,j}}\right)$
$\displaystyle+2\sum_{c}\gamma_{c}\log\frac{\pi_{c}}{\gamma_{c}}+\sum_{j}\left(1+\log\sigma^{2}_{j}\right),$

where $||\cdot||$ denotes the $\ell_{2}-$norm, subscripts denote scalar
components of tensors, $\gamma_{c}$ is the posterior distribution

$\gamma_{c}=p(c|Z)=\frac{\pi_{c}p(Z|c)}{\sum_{c^{\prime}}\pi_{c^{\prime}}p(Z|c^{\prime})},$
(15)

we estimate $q(c|X_{1},...,X_{d})=p(c|Z)=\gamma_{c}$ following (Jiang et al.,
2017), and we have taken $\hat{\sigma_{c}}=1$. The derivation may be found in
Appendix B. We seek the minimizer of this loss over the entire data set:
$\mathcal{L}=-\sum_{d}\mathcal{L}_{d}$.

In standard expectation-maximization fashion, we note that for fixed value of
$\gamma_{c}$, the variation of $\mathcal{L}$ with respect to the cluster
centers $\mu_{c}$ yields the global minimizer

$\mu_{c}=\frac{\sum_{d}\gamma_{cd}\mu_{d}}{\sum_{d}\gamma_{cd}},$ (16)

where $\mu_{d}$ and $\gamma_{cd}$ denote the encoded mean and posterior
$p(c|X_{1},...,X_{D})$ of the $d^{th}$ data point, respectively.

To facilitate batched access to large datasets, this may be calculated
incrementally following the streaming algorithm outlined in Algorithm 1,
alternating between an EM update for $\mu_{c}$ followed by an Adam update
(Kingma & Ba, 2014) of the remaining variables.

A weighted least squares problem for the optimal expert model parameters may
be similarly obtained by taking the variation of the ELBO with respect to
$\hat{\theta}_{c}$:

$\hat{\theta}_{c}=\underset{\theta^{\prime}}{\text{argmin}}\sum_{d}\sum_{t_{n}}\gamma_{cd}\left(X_{i,n}-\mathcal{E}(t_{n};\hat{\theta}^{\prime})\right)^{2}.$
(17)

Efficient solution of this nonlinear least squares problem at each epoch will
be dependent upon the problem-specific expert model and data stream, and to
perform batching may require a streaming technique such as recursive least
squares or Kalman filtering (Cioffi & Kailath, 1984). For simplicity, we
update $\hat{\theta}_{c}$ with Adam in this work but note that solving this at
each epoch to ensure the expert model provides a best fit to the current
partitions is likely to provide substantial improvement.

###  2.4 Cross-modal inference

The primary objective of this work is to perform cross-modal inference:
sampling from $q(Z|X_{i})$ allows: 1. generative modeling by decoding
$p(X_{j}|Z)$ for $i\neq j$, and 2. an estimate of $p(c|Z)$ via Eqn. (15).
Unfortunately, sampling from the unimodal encoders $q(Z|X_{i})$ provides poor
embeddings far from the multimodal embedding $q(Z|X_{1},...,X_{d})$. To remedy
this, we introduce a second set of unimodal encoders
$\tilde{q}(Z|X_{i})\sim\mathcal{N}(Z;\tilde{\mu}_{i},\tilde{\sigma}_{i}^{2}\mathbf{I})$
with identical architecture to $q(Z|X_{i})$. After the multimodal network is
trained, we minimize the KL-divergence between $\tilde{q}(Z|X_{i})$ and
${q}(Z|X_{1},...,X_{D})$ so that the unimodal embeddings reproduce the
multimodal one. In this sense, the unsupervised multimodal training provides
labels which allows supervised training of unimodal embeddings. The KL loss
admits the closed form expression

$\mathcal{L}_{i}=\frac{1}{2}\sum_{d,k}\left[\log\frac{\tilde{\sigma}^{2}_{i,kd}}{\sigma^{2}_{kd}}-l+\frac{\sigma^{2}_{kd}}{\tilde{\sigma}^{2}_{i,kd}}+\frac{\left(\tilde{\mu}_{i,kd}-\mu_{kd}\right)^{2}}{\tilde{\sigma}^{2}_{i,kd}}\right],$
(18)

which may be sequentially optimized with Adam for each modality $i$ after the
multimodal model has been fit. An additional possibility not pursued here is
to perform a Bayesian estimate to identify either the cluster most likely to
generate the data

$p(c|X_{i})=\frac{p(X_{i}|c)p(c)}{\sum_{c^{\prime}}p(X_{i}|c^{\prime})p(c^{\prime})}\,\,\text{
for }\,\,X_{i}\in\mathcal{M}_{S},$ (19)

or the cluster centroid most likely to have generated the data

$p(\mu_{c}|X_{i})=\frac{p(X_{i}|\mu_{c})p(\mu_{c})}{\sum_{c^{\prime}}p(X_{i}|\mu_{c^{\prime}})p(\mu_{c^{\prime}})}\text{
for }\,\,X_{i}\in\mathcal{M}_{DD}.$ (20)

Algorithm 1 Training with streaming EM for cluster centers

Input: data $\mathbf{x}=\\{X_{1},...,X_{d}\\}$ in batches $\mathcal{B}$

for $i=1$ to $N_{epochs}$ do

Calculate $\gamma=p(c|\mathbf{x})$ for all $\mathbf{x}$

Let $W^{new}=0$, $M=0$

for $\mathbf{b}\in\mathcal{B}$ do

$W^{old}\leftarrow W^{new}$

$W^{new}\leftarrow
W^{new}+\sum_{\mathbf{x}\in\mathbf{b}}\gamma_{c}(\mathbf{x})$

$M\leftarrow W^{old}M+\frac{\gamma_{c}\mu}{W^{new}}$

end for

$\mu_{c}\leftarrow M$

for $\mathbf{b}\in\mathcal{B}$ do

Calculate Adam update on ELBO

end for

end for

##  3 Experiments

Hyper-parameters for both training and architecture were selected using the
Weights and Biases experiment tracking tool (Biewald, 2020) and were based on
the MNIST dataset (LeCun et al., 1998) using a 90/10 train/validation split of
the training data. MNIST results are reported on the standard 10,000 held out
test examples (LeCun et al., 2010). All parameters used in the study, code and
data to reproduce experiments, and discussion of hardware and training time
for each experiment, may be found in Appendix A. We consider a non-physical
and physical pair of experiments described in Figure 2.

Figure 2: Experimental setup for unsupervised multimodal MNIST (left) and
multimodal mechanical test (right). For MNIST, we replace the labels on digits
$c\in\left\\{0,...,9\right\\}$ by a sample of the function
$X_{2}=ct+\epsilon$, for $t\in[0,1]$ and Gaussian noise $\epsilon$; c is
therefore encoded either as an image in $X_{1}$ or as the slope of a 1D
function in $X_{2}$. For the mechanical problem, a high-throughput compression
test is performed on two populations of additively manufactured lattices
corresponding to gyroid and octet microstructure (Garland et al., 2020). We
supplement the resulting stress-strain curves $X_{2}$ with images of the
microstructure $X_{1}$ to obtain low- and high-throughput modalities,
respectively.

Encoding Dim. | Mean Acc. | Min Acc. | Max Acc.
---|---|---|---
1 | 0.793 +/- 0.147 | 0.606 | 0.999
2 | 0.854 +/- 0.121 | 0.617 | 1.00
3 | 0.967 +/- 0.047 | 0.894 | 1.00
4 | 0.921 +/- 0.614 | 0.810 | 1.00

Figure 3: MNIST clusters and resulting accuracies for multimodal clean
dataset, varying encoding dimension of size 1, 2, 3, and 4 and selecting model
from hyperparameter training with lowest validation loss. Test accuracy
achieved by clusters shown: 99.9%, 100.0%, 99.9%, and 100.0%, respectively.
Statistics correspond to 10 runs using optimal hyperparameters. 3 of the 4
encoding dimensions from the 4D model are shown.

Method | CNN SotA∗ | VAE+GMM† | DEC† | VaDE† | GMVAE†† | GMVAE††
---|---|---|---|---|---|---
Notes | Supervised |  |  |  | 10 clusters | 16 clusters
Accuracy (max) | 99.91 | 72.94 | 84.30 | 94.46 | 88.54 | 96.92
Accuracy (mean±stdev) | n/a | n/a | n/a | n/a | 82.31 (3.75) | 87.82 (5.33)

Method | PIMA | PIMA | PIMA | PIMA | PIMA | PIMA
---|---|---|---|---|---|---
Notes |  | multimodal
---
low noise
10 clusters
| multimodal
---
high noise
10 clusters
| unimodal $X_{1}$
---
low noise
10 clusters
| unimodal $X_{2}$
---
low noise
10 clusters
| unimodal $X_{1}$
---
high noise
10 clusters
| unimodal $X_{2}$
---
high noise
10 clusters
Accuracy (max) | 100.0 | 98.53 | 79.62 | 99.98 | 84.85 | 78.33
Accuracy (mean±stdev) | 96.74 (4.72) | 93.15 (7.07) | - | - | - | -

Table 1: Unsupervised classification accuracy for MNIST. Results gathered from
(An et al., 2020), (Jiang et al., 2017) and (Dilokthanakul et al., 2016)
denoted by $*$,$\dagger$ and $\dagger\dagger$, respectively. If statistics
were not provided we assume maximum accuracy was reported. While the data
augmentation offered by $X_{2}$ is not incorporated in comparisons to unimodal
unsupervised benchmarks, a comparison to the supervised setting is valid. For
all experiments we do not overparameterize and keep clusters equal to the
number of digits.

Figure 4: Generative modeling for multimodal MNIST experiment. Clustering in
latent space with generative digit images $X_{1}$ at cluster centroids (left)
and expert fit to $X_{2}$ (center) for clean (top) and noisy (bottom) $X_{2}$
data. In both cases, data is embedded into disentangled clusters adjacent to
digits to obtain ordered clusters; for noisy data clusters overlap to reflect
class confusion in $X_{2}$ but preserve ordering. In this noisy case, this is
reflected when sampling $X_{1}$ from
$Z_{i}\in\left\\{\mu_{c}\text{\textpm}3\sigma_{c}\right\\}_{i}$ for clusters
1-9 (right): at the periphery of each cluster a digit ±1 is generated,
reflecting the fuzzy class boundary of $X_{2}$.

Figure 5: Confusion matrices denoting frequency with which each cluster
reproduces labels during unsupervised training. The approximate banding of the
matrix illustrates that the sequential embedding of clusters limits
misclassified digits to numbers with similar values in $X_{2}$. Left/right:
clean/noisy, top: Multimodal inference center/bottom: inference from
$X_{1}$/$X_{2}$, respectively.

###  3.1 Unsupervised multimodal MNIST

In the first experiment, we take the traditional MNIST dataset consisting of
images $X_{1}$ and labels $c\in\left\\{0,...,9\right\\}$ and manufacture a
synthetic 1D modality $X_{2}=ct+\epsilon$, where $t\in[0,1]$ and
$\epsilon\sim\mathcal{N}(0,0.01)$ for a clean dataset or
$\epsilon\sim\mathcal{N}(0,0.5)$ for a noisy data set. We adopt the affine
expert model $\mathcal{E}(t;\theta_{c})=\theta_{c}t$, and perform unsupervised
clustering of the multimodal dataset $(X_{1},X_{2})$, and additionally perform
cross-modal inference. For this artificial problem, the labels are thinly
veiled as the slope of second modality, and so we expect that if we
successfully perform multimodal inference we should obtain accuracy comparable
to a supervised MNIST benchmark.

We define unsupervised clustering accuracy ($acc$) as in (Xie et al., 2016),
(Jiang et al., 2017):

$acc=\underset{m\in
M}{\text{max}}\frac{\sum_{i=1}^{N}\mathbbm{1}\\{l_{i}=m(c_{i})\\}}{N},$ (21)

where $N$ is the number of examples, $M$ is the set of all possible mappings
from a cluster to a label assignment, $l_{i}$ is the true label and $c_{i}$ is
the cluster assignment by the model for example $i$. Calibration of the latent
dimension (Figure 3) revealed $l=3$ to be an optimal embedding dimension. In
Table 1 we provide a comparison to classification accuracy against state of
the art supervised and unsupervised models trained on images only. For
multimodal inference with low noise we obtain perfect classification
outperforming the current state-of-the-art in supervised classification and
outperforming the state-of-the-art in unsupervised learning.

In addition to observing improved accuracy, we investigate in Figure 4 the
disentangled representation of clusters and provide examples of generative
models. Surprisingly, the sequential ordering of clusters in $X_{2}$ induces
an sequential embedding of corresponding clusters in latent space. As noise is
increased, the data in $X_{2}$ is sufficiently large to prevent distinct
clustering into disentangled classes, however the ordering of digits is
roughly preserved. When generative modeling is performed for $X_{2}$, this is
reflected by samples from the tail of the Gaussian mixture which generate
images of adjacent digits. For example, the cluster of 2’s contains images of
3’s in its periphery. This suggests that the sequential ordering of data in
$X_{2}$ may induce generative models for $X_{1}$ which reflect cross-modal
ordering. Confusion matrices for multimodal and cross-modal classification in
Figure 5 reveal an approximately banded structure, whereby misclassified
modalities primarily occur between adjacent digits. While this single
experiment is insufficient to remark on the generality of this result, it
suggests the potential for learning generative models of images which reflect
information from the expert model, a particularly exciting prospect for
scientific datasets.

###  3.2 Metal additive lattice fingerprinting

Figure 6: Generative modeling for metal additive dataset. Top left: Data is
clustered into two gaussian distributions with mean visualized by black dot.
Shown here is test data held out during training. Bottom left: Stress-strain
curves from two populations split naturally between two populations, and
linear-strain hardening expert model provides provides best fit (orange and
blue lines). Right: Generative sampling from the means of each cluster
provides a prediction for microstructure associated with either octet (bottom)
or gyroid (top) lattice topologies.

The lattice dataset consists of $X_{1}$ images of 3D printed lattices split
between two types of printed metamaterials and corresponding $X_{2}$
stress/strain curves performed in a high-throughput uniaxial compression
machine. Even with high-throughput testing, only 91 pairs of images were able
to be generated, highlighting the utility of high-throughput photographs
$X_{1}$ as surrogates to infer $X_{2}$. We select as expert model a linear
strain-hardening model partitioning the stress-strain response into two
piecewise linear regions (Jones, 2009). Even a simple model like this
succinctly encodes a large number of interpretable quantities of interest: a
yield stress, together with elastic and plastic moduli. We gather in Figure 6
results from generative modeling and include in the appendix additional
results demonstrating 94.74%/94.74%/94.74% classification accuracy of the two
clusters for multimodal/$X_{1}$/$X_{2}$-cross-modal inferences, respectively.

##  4 Conclusions and future work

The present approach provides an abstract variational inference for
discovering fingerprints in an unsupervised manner while incorporating
physical model biases. This framework is widely applicable to a range of
scientific disciplines where detection of fingerprints are crucial for tasks
ranging from predicting and attributing climate change to designing
biochemical pathways at a molecular level. In addition to the application
focus on fingerprint generation, this framework may be used for a variety of
general purpose downstream tasks based on multimodal processing of scientific
data. For this work we have focused on a simple MNIST example to probe
dynamics for an easily replicable and understandable dataset, along with a
simple high-throughput manufacturing example to illustrate feasibility for
facilitating high-throughput experimentation.

In future work, we will employ more sophisticated physics-informed surrogates
as expert models for processes in additive manufacturing and semiconductor
device design bridging multiscale and multifidelity information from data
sources corresponding to both physical experiment (transmission electron
microscopy, atom probe tomography, micro x-ray computerized tomography, or
synchotron x-ray diffraction) and high-fidelity simulations (quantum density
functional theory, molecular dynamics, crystal plasticity, and continuum
mechanics). To date, these costly but rich sources of information are
antagonistic to high-throughput testing and simulation. This framework
provides an exciting platform for discovering data-driven scientific
fingerprints which may be combined with advances in automated experimentation
to accelerate scientific discovery.

## Software and Data

Pending acceptance, all data and code used to generate results will be hosted
on a Github page. For now, we include a subset of the code to reproduce the
MNIST examples through the anonymous
Github111https://anonymous.4open.science/r/PIMA-5D6B/; unfortunately anonymous
github does not offer enough storage to host the accompanying dataset.

## Acknowledgements

The authors thank Warren Davis, Anthony Garland and Lekha Patel for providing
guidance on the variational inference framework and review of the manuscript
and Kat Reiner and Greg Geller for providing computing support. All authors
acknowledge funding under the Beyond Fingerprinting Sandia Grand Challenge
Laboratory Directed Research and Developement program. N. Trask acknowledges
funding under the Collaboratory on Mathematics and Physics-Informed Learning
Machines for Multiscale and Multiphysics Problems (PhILMs) project funded by
DOE Office of Science (Grant number DE-SC001924) and the DOE Early Career
program. Sandia National Laboratories is a multi-mission laboratory managed
and operated by National Technology and Engineering Solutions of Sandia, LLC.,
a wholly owned subsidiary of Honeywell International, Inc., for the U.S.
Department of Energy’s National Nuclear Security Administration under contract
DE-NA0003525. This paper describes objective technical results and analysis.
Any subjective views or opinions that might be expressed in the paper do not
necessarily represent the views of the U.S. Department of Energy or the United
States Government. SAND number: SAND2022-1159 O

## References

  * Abadi et al. (2016) Abadi, M., Barham, P., Chen, J., Chen, Z., Davis, A., Dean, J., Devin, M., Ghemawat, S., Irving, G., Isard, M., et al.  Tensorflow: A system for large-scale machine learning.  In _12th $\\{$USENIX$\\}$ symposium on operating systems design and implementation_, pp. 265–283, 2016.
  * An et al. (2020) An, S., Lee, M., Park, S., Yang, H., and So, J.  An ensemble of simple convolutional neural network models for MNIST digit recognition.  _arXiv preprint arXiv:2008.10400_ , 2020.
  * Baltrušaitis et al. (2018) Baltrušaitis, T., Ahuja, C., and Morency, L.-P.  Multimodal machine learning: A survey and taxonomy.  _IEEE transactions on pattern analysis and machine intelligence_ , 41(2):423–443, 2018.
  * Biewald (2020) Biewald, L.  Experiment tracking with weights and biases, 2020.  URL https://www.wandb.com/.  Software available from wandb.com.
  * Boyce & Uchic (2019) Boyce, B. L. and Uchic, M. D.  Progress toward autonomous experimental systems for alloy development.  _MRS Bulletin_ , 44(4):273–280, 2019.  ISSN 0883-7694.
  * Burgess et al. (2018) Burgess, C. P., Higgins, I., Pal, A., Matthey, L., Watters, N., Desjardins, G., and Lerchner, A.  Understanding disentangling in $\beta$-VAE.  _arXiv preprint arXiv:1804.03599_ , 2018.
  * Chakraborty et al. (2021) Chakraborty, A., Nandi, P., and Chakraborty, B.  Fingerprints of the quantum space-time in time-dependent quantum mechanics: An emergent geometric phase.  _arXiv preprint arXiv:2110.04370_ , 2021.
  * Chen et al. (2018) Chen, R. T., Li, X., Grosse, R., and Duvenaud, D.  Isolating sources of disentanglement in vaes.  In _Proceedings of the 32nd International Conference on Neural Information Processing Systems_ , pp. 2615–2625, 2018.
  * Cioffi & Kailath (1984) Cioffi, J. and Kailath, T.  Fast, recursive-least-squares transversal filters for adaptive filtering.  _IEEE Transactions on Acoustics, Speech, and Signal Processing_ , 32(2):304–337, 1984.
  * Dilokthanakul et al. (2016) Dilokthanakul, N., Mediano, P. A., Garnelo, M., Lee, M. C., Salimbeni, H., Arulkumaran, K., and Shanahan, M.  Deep unsupervised clustering with Gaussian mixture variational autoencoders.  _arXiv preprint arXiv:1611.02648_ , 2016.
  * Garland et al. (2020) Garland, A. P., White, B. C., Jared, B. H., Heiden, M., Donahue, E., and Boyce, B. L.  Deep convolutional neural networks as a rapid screening tool for complex additively manufactured structures.  _Additive Manufacturing_ , 35:101217, 2020.
  * Hasselmann (1997) Hasselmann, K.  Multi-pattern fingerprint method for detection and attribution of climate change.  _Climate dynamics_ , 13(9):601–611, 1997.
  * Hegerl et al. (2007) Hegerl, G., Zwiers, F., Braconnot, P., Gillett, N. P., Luo, Y. M., Orsini, J. M., Nicholls, N., Penner, J. E., and Stott, P. A.  Understanding and attributing climate change, 2007.
  * Higgins et al. (2017) Higgins, I., Matthey, L., Pal, A., Burgess, C., Glorot, X., Botvinick, M., Mohamed, S., and Lerchner, A.  beta-vae: Learning basic visual concepts with a constrained variational framework.  In _5th International Conference on Learning Representations, ICLR 2017_ , 2017.
  * Hunter (2007) Hunter, J. D.  Matplotlib: A 2d graphics environment.  _Computing in science & engineering_, 9(03):90–95, 2007.
  * Isayev et al. (2015) Isayev, O., Fourches, D., Muratov, E. N., Oses, C., Rasch, K., Tropsha, A., and Curtarolo, S.  Materials cartography: representing and mining materials space using structural and electronic fingerprints.  _Chemistry of Materials_ , 27(3):735–743, 2015\.
  * Jang et al. (2016) Jang, E., Gu, S., and Poole, B.  Categorical reparameterization with gumbel-softmax.  _arXiv preprint arXiv:1611.01144_ , 2016.
  * Jiang et al. (2017) Jiang, Z., Zheng, Y., Tan, H., Tang, B., and Zhou, H.  Variational deep embedding: an unsupervised and generative approach to clustering.  In _Proceedings of the 26th International Joint Conference on Artificial Intelligence_ , pp. 1965–1972, 2017.
  * Jones (2009) Jones, R. M.  _Deformation theory of plasticity_.  Bull Ridge Corporation, 2009.
  * Jordan & Jacobs (1994) Jordan, M. I. and Jacobs, R. A.  Hierarchical mixtures of experts and the em algorithm.  _Neural computation_ , 6(2):181–214, 1994.
  * Kim & Mnih (2018) Kim, H. and Mnih, A.  Disentangling by factorising.  In _International Conference on Machine Learning_ , pp. 2649–2658. PMLR, 2018.
  * Kingma & Ba (2014) Kingma, D. P. and Ba, J.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ , 2014.
  * Kingma & Welling (2014) Kingma, D. P. and Welling, M.  Auto-Encoding Variational Bayes.  In _2nd International Conference on Learning Representations, ICLR 2014_ , 2014.
  * Kuhn (1955) Kuhn, H. W.  The hungarian method for the assignment problem.  _Naval Research Logistics Quarterly_ , 2(1-2):83–97, 1955.
  * Lagaris et al. (1998) Lagaris, I. E., Likas, A., and Fotiadis, D. I.  Artificial neural networks for solving ordinary and partial differential equations.  _IEEE transactions on neural networks_ , 9(5):987–1000, 1998.
  * LeCun et al. (1998) LeCun, Y., Bottou, L., Bengio, Y., and Haffner, P.  Gradient-based learning applied to document recognition.  _Proceedings of the IEEE_ , 86(11):2278–2324, 1998.
  * LeCun et al. (2010) LeCun, Y., Cortes, C., and Burges, C.  Mnist handwritten digit database.  _ATT Labs [Online]. Available: http://yann.lecun.com/exdb/mnist_ , 2, 2010.
  * Lee et al. (2020) Lee, D. B., Min, D., Lee, S., and Hwang, S. J.  Meta-GMVAE: Mixture of Gaussian VAE for unsupervised meta-learning.  In _International Conference on Learning Representations_ , 2020.
  * Liu et al. (2012) Liu, A., Zhu, W., Tsai, D., and Zheludev, N. I.  Micromachined tunable metamaterials: a review.  _Journal of Optics_ , 14(11):114009, 2012.
  * Locatello et al. (2019) Locatello, F., Bauer, S., Lucic, M., Raetsch, G., Gelly, S., Schölkopf, B., and Bachem, O.  Challenging common assumptions in the unsupervised learning of disentangled representations.  In _International Conference on Machine Learning_ , pp. 4114–4124. PMLR, 2019.
  * Lu et al. (2019) Lu, L., Jin, P., and Karniadakis, G. E.  Deeponet: Learning nonlinear operators for identifying differential equations based on the universal approximation theorem of operators.  _arXiv preprint arXiv:1910.03193_ , 2019.
  * Mao et al. (2021) Mao, Z., Lu, L., Marxen, O., Zaki, T. A., and Karniadakis, G. E.  Deepm&mnet for hypersonics: Predicting the coupled flow and finite-rate chemistry behind a normal shock using neural-network approximation of operators.  _Journal of Computational Physics_ , 447:110698, 2021.
  * Pedregosa et al. (2011) Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, O., Blondel, M., Prettenhofer, P., Weiss, R., Dubourg, V., et al.  Scikit-learn: Machine learning in python.  _The Journal of Machine Learning Research_ , 12:2825–2830, 2011.
  * Pu et al. (2016) Pu, Y., Gan, Z., Henao, R., Yuan, X., Li, C., Stevens, A., and Carin, L.  Variational autoencoder for deep learning of images, labels and captions.  _Advances in neural information processing systems_ , 29:2352–2360, 2016.
  * Raissi et al. (2019) Raissi, M., Perdikaris, P., and Karniadakis, G. E.  Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations.  _Journal of Computational Physics_ , 378:686–707, 2019\.
  * Rao et al. (2019) Rao, D., Visin, F., Rusu, A., Pascanu, R., Teh, Y. W., and Hadsell, R.  Continual unsupervised representation learning.  _Advances in Neural Information Processing Systems_ , 32:7647–7657, 2019.
  * Rezende et al. (2014) Rezende, D. J., Mohamed, S., and Wierstra, D.  Stochastic backpropagation and approximate inference in deep generative models.  In _International Conference on Machine Learning_ , pp. 1278–1286. PMLR, 2014.
  * Shi et al. (2019) Shi, Y., Paige, B., Torr, P., et al.  Variational mixture-of-experts autoencoders for multi-modal deep generative models.  _Advances in Neural Information Processing Systems_ , 32:15718–15729, 2019.
  * Sochol et al. (2018) Sochol, R. D., Sweet, E., Glick, C. C., Wu, S.-Y., Yang, C., Restaino, M., and Lin, L.  3d printed microfluidics and microelectronics.  _Microelectronic Engineering_ , 189:52–68, 2018.
  * Sohn et al. (2015) Sohn, K., Lee, H., and Yan, X.  Learning structured output representation using deep conditional generative models.  _Advances in neural information processing systems_ , 28:3483–3491, 2015.
  * Sutter et al. (2021) Sutter, T. M., Daunhawer, I., and Vogt, J. E.  Generalized multimodal ELBO.  In _9th International Conference on Learning Representations, ICLR_ , 2021.
  * Suzuki et al. (2017) Suzuki, M., Nakayama, K., and Matsuo, Y.  Joint multimodal learning with deep generative models.  In _5th International Conference on Learning Representations, ICLR 2017_ , 2017.
  * Trask et al. (2020) Trask, N., Huang, A., and Hu, X.  Enforcing exact physics in scientific machine learning: a data-driven exterior calculus on graphs.  _arXiv preprint arXiv:2012.11799_ , 2020.
  * Vedantam et al. (2018) Vedantam, R., Fischer, I., Huang, J., and Murphy, K.  Generative models of visually grounded imagination.  In _6th International Conference on Learning Representations, ICLR_ , 2018.
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., Haberland, M., Reddy, T., Cournapeau, D., Burovski, E., Peterson, P., Weckesser, W., Bright, J., et al.  Scipy 1.0: fundamental algorithms for scientific computing in python.  _Nature methods_ , 17(3):261–272, 2020.
  * Wang et al. (2021) Wang, S., Wang, H., and Perdikaris, P.  Learning the solution operator of parametric partial differential equations with physics-informed deeponets.  _arXiv preprint arXiv:2103.10974_ , 2021.
  * Waskom (2021) Waskom, M. L.  seaborn: statistical data visualization.  _Journal of Open Source Software_ , 6(60):3021, 2021.  doi: 10.21105/joss.03021.  URL https://doi.org/10.21105/joss.03021.
  * Wu & Goodman (2018) Wu, M. and Goodman, N.  Multimodal generative models for scalable weakly-supervised learning.  _arXiv preprint arXiv:1802.05335_ , 2018.
  * Xie et al. (2016) Xie, J., Girshick, R., and Farhadi, A.  Unsupervised deep embedding for clustering analysis.  In _International Conference on Machine Learning_ , pp. 478–487. PMLR, 2016.

##  Appendix A Architecture, hyperparameters, and implementation

Model Architectures

We employ relatively small convolutional architectures to serve as encoders
for both modalities. The image modality encoder consists of 2 2D convolutional
layers with 32 and 64 channels respectively, each with 3x3 kernels. We apply
the exponential linear unit (ELU) activation function as well as batch
normalization after each convolutional layer, then pass the output to a fully
connected layer of size $encoding_{d}im\times 2$ to enable an embedding into a
representation of the mean and standard deviation of the input in the latent
space. For the 1D modality encoder, in place of 2D convolutional layers, we
use 1D convolutions with 8 and 16 channels in the respective layers, but with
an otherwise identical architecture.

The image decoder begins with a fully connected layer of appropriate size to
be reshaped into 32 channels of 2D arrays, with each dimension having a length
$\frac{1}{4}$ of the length of the number of pixels per side of the original
image. The reshaped output of the initial dense layer is passed into a series
of 3 deconvolutional layers with 64, 32, and 1 channel respectively, each with
a kernel size of 3. The first 2 deconvolutional layers use a stride of 3 and a
Rectified Linear Unit (ReLU) activation function, and the final
deconvolutional layer uses a stride of 1. Zero padding is used to retain the
input shape while traversing these layers.

Hyperparameters

We used the MNIST dataset along with the Weights and Biases tool (Biewald,
2020) to perform a hyperparameter search over learning rates and encoding
dimension size. For the multimodal models, we found that a learning rate of
$1.97e-5$ gave the best result as measured by validation loss. For the MNIST
dataset, an encoding dimension of size 3 gave the most consistent results, and
we used a 2D encoding dimension for the lattice dataset. For the unimodal
models, we found that a learning rate of $4.398e-5$ for the image modality and
a learning rate of $4.398e-3$ for the 1D modality supported learning, but we
leave a rigorous hyperparameter tuning for unimodal models for future work.

Implementation details

Each of the input data modalities was normalized to have values in $[0,1]$.
The 1D data was sampled to generate an array of length 20 for MNIST and length
100 for the lattice stress-strain data. The lattice images were cropped and
subsampled into quadrants resulting in images of dimension 152x152. We further
augmented the dataset by flipping the images along each axis.

To train the MNIST multimodal models, we used 10% of the standard MNIST train
set for validation, and selected the model with the lowest validation loss. We
did not observe any indication of overfitting the multimodal model to the
training data. All results reported are on the held out test dataset. For the
MNIST unimodal models, we again did not observe overfitting when measuring
model performance by accuracy to the cluster labels generated from the trained
multimodal model, and we report results on the test set from the model
resulting from the final training epoch.

We split the lattice dataset into a train/test 80/20 split. Since we did not
observe overfitting during multimodal training, we selected the model with the
lowest training loss and applied that model to the held out test set, and we
report those results. We select the model resulting from the final training
epoch for the unimodal tasks.

Our models are implemented in Python using Tensorflow (Abadi et al., 2016),
and we leverage the Scikit-learn library (Pedregosa et al., 2011) for data
preparation and accuracy metrics, Scipy (Virtanen et al., 2020) for data
preparation and the linear_sum_assignment implementation of the Hungarian
method (Kuhn, 1955) for efficient computation of the unsupervised cluster
accuracy. We visualize our results using the Matplotlib (Hunter, 2007) and
Seaborn (Waskom, 2021) libraries. Training was performed on NVIDIA DGX-1 and
DGX-2 machines with each run executed on 1 GPU. A combination of P100 and V100
GPUs were used in this work. Training runs took on average approximately
(machine dependent) 10-14 hours with our longest run on lattice data taking
approximately one day. We made no attempt to optimize parallel training to
improve run time or training efficiency.

##  Appendix B Derivation of ELBO

To derive a closed form expression for the single sample ELBO

$\mathcal{L}_{d}=\mathbb{E}_{q(Z,c|X_{1},...,X_{D})}\left[\log\frac{p(X_{1},...,X_{D},Z)}{q(Z,c|X_{1},...,X_{D})}\right]$
(22)

we apply the separability assumptions

$\displaystyle
p(X_{1},...,X_{D},Z,c)=\left(\Pi_{i=1}^{D}p(X_{i}|Z,c)\right)p(Z|c)p(c),$ (23)
$\displaystyle
q(Z,c|X_{1},...,X_{D})=q(Z|X_{1},...,X_{D})q(c|X_{1},...,X_{D}).$ (24)

providing the additive decomposition

$\mathcal{L}=\mathbb{E}_{q(Z,c|X_{1},...,X_{D})}\left[\log
p(X_{1},...,X_{D},Z)\right]-\mathbb{E}_{q(Z,c|X_{1},...,X_{D})}\left[\log
q(Z,c|X_{1},...,X_{D})\right]$ (25)
$=\sum_{i=1}^{D}\mathbb{E}_{q(Z,c|X_{1},...,X_{D})}\left[\log
p(X_{i}|Z,c)\right]+\mathbb{E}_{q(Z,c|X_{1},...,X_{D})}\left[\log
p(Z|c)\right]$ $+\mathbb{E}_{q(Z,c|X_{1},...,X_{D})}\left[\log
p(c)\right]-\mathbb{E}_{q(Z,c|X_{1},...,X_{D})}\left[\log
q(Z|X_{1},...,X_{D})\right]$ $-\mathbb{E}_{q(Z,c|X_{1},...,X_{D})}\left[\log
q(c|X_{1},...,X_{D})\right].$

For convenience we denote
$\mathbb{E}_{q(Z,c|X_{1},...,X_{D})}=\mathbb{E}_{q}$. The separability
assumptions therefore decompose the ELBO into constituent expectations of the
form

$\mathbb{E}_{q}\left[\log
f(Z,c)\right]=\sum_{c}\int_{\mathbb{R}^{l}}f(Z,c)\log g(Z,c)\,dZ$ (26)

which may be integrated exactly for the Gaussian/categorical $f$ and $g$
appearing in the ELBO. The only term which may not be immediately computed is
$\mathbb{E}_{q}\left[\log q(c|X_{1},...,X_{D})\right]$. The lack of a
reparameterization trick for the categorical distribution precludes
backpropagation into the encoder, forcing us to consider an encoder which only
provides predictions for $Z$. While there are options to use e.g. a
regularized Gumbel-softmax approximation to the categorical distribution (Jang
et al., 2016), we would lose the tractability of the closed form expression
for the ELBO. Instead we follow (Jiang et al., 2017) and approximate
$q(c|X_{1},...,X_{D})=p(c|Z)$ using the following justification.

Rewriting the ELBO

$\mathcal{L}_{d}=\mathbb{E}_{q(Z,c|X_{1},...,X_{D})}\left[\log\frac{p(X_{1},...,X_{D},Z)}{q(Z,c|X_{1},...,X_{D})}\right]$
(27)
$=\mathbb{E}_{q(Z,c|X_{1},...,X_{D})}\left[\log\frac{p(X_{1},...,X_{D}|Z)p(Z)}{q(Z|X_{1},...,X_{D})}+\log\frac{p(c|Z)}{q(c|X_{1},...,X_{D})}\right]$
$=\int_{\mathbb{R}^{l}}\log\frac{p(X_{1},...,X_{D}|Z)p(Z)}{q(Z|X_{1},...,X_{D})}dZ+D_{KL}(q(c|X_{1},...,X_{D})||p(c|Z)),$

we seek extremal points with respect to $c$. The first term is independent of
$c$, and the positive second term takes zero value when
$q(c|X_{1},...,X_{D})=p(c|Z)$, providing the desired maximum. We caution
however that this holds only at local minima of the loss landscape, but
empirically has been shown to perform well as an estimator.

Finally for completeness, we gather from (Jiang et al., 2017) the various
integral formulas required to compute the expectations in closed form with
modifications for our multimodal setting.

######  Lemma B.1.

Given Gaussian distributions
$f(Z)=\mathcal{N}(Z;\mu_{1},\sigma_{1}\mathbf{I})$ and
$g(Z)=\mathcal{N}(Z;\mu_{2},\sigma_{2}\mathbf{I})$

$\int_{\mathbb{R}^{l}}f(Z)\log g(Z)\,dZ=-\frac{1}{2}\sum_{j}\log
2\pi\sigma_{2,j}^{2}+\frac{\sigma_{2,j}^{2}}{\sigma_{1,j}^{2}}+\frac{(\mu_{1,j}-\mu_{2,j})^{2}}{\sigma_{2,j}^{2}}.$
(28)

######  Lemma B.2.

$\mathbb{E}_{q(Z,c|X_{1},...,X_{D})}\left[\log
p(Z|c)\right]=\sum_{c}q(c|X_{1},...,X_{D})\int_{\mathbb{R}^{l}}q(Z|X_{1},...,X_{D})\log
p(Z|c)dZ,$ (29)

where the integrand may be computed from Lemma 28.

######  Lemma B.3.

$\mathbb{E}_{q(Z,c|X_{1},...,X_{D})}\left[\log
p(Z|c)\right]=\int_{\mathbb{R}^{l}}q(Z|X_{1},...,X_{D})dZ\sum_{c}q(c|X_{1},...,X_{D})\pi_{c},$
(30) $=\sum_{c}q(c|X_{1},...,X_{D})\pi_{c}.$

######  Lemma B.4.

$\mathbb{E}_{q(Z,c|X_{1},...,X_{D})}\left[\log
q(Z|X_{1},...,X_{D})\right]=\sum_{c}q(c|X_{1},...,X_{D})\int_{\mathbb{R}^{l}}q(Z|X_{1},...,X_{D})\log
q(Z|X_{1},...,X_{D})dZ,$ (31) $=\int_{\mathbb{R}^{l}}q(Z|X_{1},...,X_{D})\log
q(Z|X_{1},...,X_{D})dZ,$

where the integrand may be computed from Lemma 28.

######  Lemma B.5.

$\mathbb{E}_{q(Z,c|X_{1},...,X_{D})}\left[\log
q(c|X_{1},...,X_{D})\right]=\int_{\mathbb{R}^{l}}q(Z|X_{1},...,X_{D})dZ\sum_{c}q(c|X_{1},...,X_{D})\log
q(c|X_{1},...,X_{D}),$ (32) $=\sum_{c}q(c|X_{1},...,X_{D})\log
q(c|X_{1},...,X_{D}).$

For all terms, $q(c|X_{1},...,X_{D})$ is calculated via the posterior
estimator $\gamma$ given in Equation (15).