figuret

# Geometric Phase-Driven Scattering Evolutions

Pengxiang Wang  School of Optical and Electronic Information, Huazhong
University of Science and Technology, Wuhan, Hubei 430074, P. R. China
Yuntian Chen  Email<EMAIL_ADDRESS>School of Optical and Electronic
Information, Huazhong University of Science and Technology, Wuhan, Hubei
430074, P. R. China  Wuhan National Laboratory for Optoelectronics, Huazhong
University of Science and Technology, Wuhan, Hubei 430074, P. R. China  Wei
Liu  Email<EMAIL_ADDRESS>College for Advanced Interdisciplinary
Studies, National University of Defense Technology, Changsha 410073, P. R.
China.  Nanhu Laser Laboratory and Hunan Provincial Key Laboratory of Novel
Nano-Optoelectronic Information Materials and Devices, National University of
Defense Technology, Changsha 410073, P. R. China.

###### Abstract

We explore the classical topic of scattering manipulation, from a different
perspective of controlled excitations and interferences of quasi-normal modes
(QNMs). Scattered waves can be expanded as coherent additions of radiations
from the QNMs excited, and thus relative amplitudes and phases among them are
crucial factors to engineer for scattering shaping. Here relying on the
electromagnetic reciprocity, we provide full geometric representations based
on the Poincaré sphere for those factors, and identify the hidden underlying
geometric phases that drive the scattering evolutions. Further synchronous
exploitations of the incident polarization-dependent geometric phases and
excitation amplitudes enable efficient manipulations of both scattering
intensities and polarizations. Continuous geometric phase spanning $2\pi$ is
directly manifest through scattering variations, even in the rather elementary
configuration of an individual particle scattering waves of varying
polarizations. We have essentially merged three vibrant fields of geometric
phase, Mie scattering and QNM, and unlocked an extra dimension of geometric
phase for scattering manipulations, which will greatly broaden the horizons of
many disciplines associated with not only electromagnetic scatterings, but
also scatterings of waves in other forms.

The seminal topic of electromagnetic scatterings by particles has been the
cornerstone for investigations of light-matter interactions and various
scattering-related applications [1, 2, 3, 4, 5, 6, 7]. Stimulated by the
vibrantly developing fields of topological, non-hermitian, and singular
photonics [8, 9, 10, 11, 12, 13], their core and sweeping concepts of
topology, non-hermiticity, and singularity have recently been rapidly
introduced into this field, rendering unprecedented opportunities for many
associated disciplines [14, 15, 16, 17].

Despite its rather long history and the aforementioned conceptual advances,
the central mathematical and physical technique for the field of particle
scatterings remains to be spherical harmonics and electromagnetic multipolar
expansions [1, 18]. Though many breakthroughs in this field have been made
based on this technique (such as recent introductions of Poincaré-Hopf theorem
[19], electromagnetic multipolar parity [20] and duality [21] into Mie theory
to reveal its intrinsic topological and geometric structures [22, 23, 24, 25,
26]), the language of spherical harmonics and electromagnetic multipoles is
more descriptive than predictive. That is, except for some special scenarios
of particles with ideal spherical or cylindrical symmetries, this technique
describes the already known fields (either the near or scattered far fields
calculated with other numerical methods) rather than predict them. For
example, even for the elementary case of plane waves scattered by a particle
of arbitrary shape, knowing the multipolar components of one scattering
configuration barely tells anything about the scatterings of another even
neighbouring configuration with a slightly different incident direction and/or
polarization. Moreover, formulisms based on spherical harmonics are usually
cumbersome and tend to obscure rather than clarify the profound physical
picture. To circumvent those limitations of the conventional method and
further advance this seminal field, new concepts and techniques have to be
introduced.

Here we investigate the problem of electromagnetic scatterings from a
different perspective of engineered QNM [27] excitations and interferences.
Scattered waves by the particles can be expanded into QNM radiations, and thus
relative amplitudes and phases among them would decide the scattering
patterns. Relying on the principle of reciprocity [28], we manage to provide
full geometric representations for the excitation coefficients and discover
the geometric phase (Pancharatnam-Berry phase) [29, 30, 31, 32] that varies
with changing incident polarizations. We further exploit the unveiled
geometric phase and the also incident polarization-dependent excitation
amplitudes for scattering intensity and polarization manipulations, such as
eliminating total and directional scatterings and designing directions of
polarization singularities. Continuous geometric phases from $0$-$2\pi$ are
directly manifest through scattering variations, even in the rather simple
configuration of an individual particle scattering plane waves of varying
polarizations. Our work has uncovered a hidden dimension of freedom for
scattering manipulations, which can potentially accelerate both fundamental
explorations and practical applications relying on scatterings, not only of
electromagnetic waves, but also of waves of other forms in which the geometric
phase would generically emerge.

Throughout this work, we study reciprocal nonmagnetic scatterers of spatial
relative permittivity distribution $\mathbf{\epsilon}(\mathbf{r},\omega)$ in
background of refractive index $1$, where $\omega$ is the angular frequency.
The scatterer supports a set of discrete QNMs [27], and its scattered fields
can be expanded as:

$\mathbf{E}_{\rm{sca}}(\mathbf{{r}})=\sum\alpha_{m}{\mathbf{E}}_{m}(\mathbf{{r}}),$
(1)

where ${\mathbf{E}}_{m}(\mathbf{{r}})$ denotes the radiation of the $m^{th}$
QNM and $\alpha_{m}$ is the complex expanding (excitation) coefficient. In the
far field, both $\mathbf{E}_{\rm{sca}}(\mathbf{{r}})$ and
${\mathbf{E}}_{m}(\mathbf{{r}})$ are transverse, and thus Eq. (1) can be
reformulated as
$\rm{E}_{\rm{sca}}(\mathbf{\hat{r}})\mathbf{{J}}_{sca}(\mathbf{\hat{r}})=\sum\alpha_{m}(\omega){{\rm{E}}}_{m}(\mathbf{\hat{r}}){\mathbf{{J}}}_{m}(\mathbf{\hat{r}})$.
Here $\mathbf{\hat{r}}$ is unit direction vector
$\mathbf{\hat{r}}=\mathbf{r}/|\mathbf{r}|$; $\rm{E}_{\rm{sca}}$ and
$\rm{E}_{m}$ are field amplitudes; $\mathbf{{J}}_{\rm{sca}}(\mathbf{\hat{r}})$
and $\mathbf{{J}}_{m}(\mathbf{\hat{r}})$ are the corresponding unit Jones
(row) vectors [33]. Both scattering intensity and polarization distributions
are dictated by the relative amplitudes and phases among the excitation
coefficients $\alpha_{m}$, which are decided by the incident wave (initial
condition). For an incident (along $\mathbf{\hat{r}}_{\rm{inc}}$) plane wave
of Jones vector $\mathbf{{J}}_{\rm{inc}}$, based on reciprocity it is recently
revealed that $\alpha_{m}$ can be calculated directly in the far field [34]:

$\alpha_{m}\propto\mathbf{J_{\rm{inc}}}\mathbf{J}_{m}^{\dagger}(-\mathbf{\hat{r}}_{\rm{inc}}),$
(2)

where $\dagger$ denotes combined operations of complex conjugate and
transpose, and thus $\mathbf{J}_{m}^{\dagger}(-\mathbf{\hat{r}}_{\rm{inc}})$
is the corresponding column Jones vector for the QNM radiation along
$-\mathbf{\hat{r}}_{\rm{inc}}$ (opposite to the incident direction).

To reveal the core principles, we start with the simplest scenario of two-QNM
(denoted by A and B) excitations and the framework established can be
naturally generalized to deal with multi-mode cases. To compact the notations,
we simplify $\mathbf{J}_{\rm{A,B}}(-\mathbf{\hat{r}}_{\rm{inc}})$ as
$\mathbf{J}_{\rm{A,B}}$, and $\mathbf{J}_{\rm{A,B}}(\mathbf{\hat{r}})$ as
$\mathbb{J}_{\rm{A,B}}$. Then scattering intensity along $\mathbf{\hat{r}}$
can be expressed as [35, 36, 37]:

$\begin{split}I_{\rm{sca}}\propto|\mathbf{J_{\rm{inc}}}\mathbf{J}^{\dagger}_{\rm{A}}|^{2}~{}I_{\rm{A}}+|\mathbf{J_{\rm{inc}}}\mathbf{J}_{\rm{B}}^{\dagger}|^{2}~{}I_{\rm{B}}+\\\
2|\mathbf{J_{\rm{inc}}}\mathbf{J}^{\dagger}_{\rm{A}}||\mathbf{J_{\rm{inc}}}\mathbf{J}_{\rm{B}}^{\dagger}|\sqrt{I_{\rm{A}}I_{\rm{B}}}|\mathbb{J}_{\rm{A}}\mathbb{J}_{\rm{B}}^{\dagger}|\cos[\varphi(\mathbf{\hat{r}})],\end{split}$
(3)

where $\varphi(\mathbf{\hat{r}})$ denotes the total phase difference between
the directional radiations of the two modes, and $I_{\rm{A,B}}$ denotes their
respective intensity:
$I_{\rm{A,B}}=|{\mathbf{E}}_{\rm{A,B}}(\mathbf{\hat{r}})|^{2}$. To geometrize
the relative excitation amplitude and phase, we map the Jones vectors of
$\mathbf{J_{\rm{inc}}}$, $\mathbf{J}_{\rm{A}}$ and $\mathbf{J}_{\rm{B}}$ to
three points P, A, B on the Poincaré sphere [33, 35, 36, 37] of unit radius
[shown in Fig. 1(a)]. Then Eq. (3) can be written in a geometric form as [35,
36, 37]:

$\begin{split}I_{\rm{sca}}\propto\cos^{2}(\frac{1}{2}{\ooalign{\leavevmode\resizebox{13.47224pt}{0.0pt}{{\phantom{A}}}\cr
PA}})~{}I_{\rm{A}}+\cos^{2}(\frac{1}{2}{\ooalign{\leavevmode\resizebox{13.88893pt}{0.0pt}{{\phantom{A}}}\cr
PB}})~{}I_{\rm{B}}\\\
+2\cos(\frac{1}{2}{\ooalign{\leavevmode\resizebox{13.47224pt}{0.0pt}{{\phantom{A}}}\cr
PA}})\cos(\frac{1}{2}{\ooalign{\leavevmode\resizebox{13.88893pt}{0.0pt}{{\phantom{A}}}\cr
PB}})~{}\sqrt{I_{\rm{A}}I_{\rm{B}}}|\mathbb{J}_{\rm{A}}\mathbb{J}_{\rm{B}}^{\dagger}|\cos[\varphi(\mathbf{\hat{r}})],\end{split}$
(4)

where  A PA (  A PB ) denotes the length of the great arc (shorter segment
$\in[0,\pi]$) connecting PA (PB). When the incident polarization is orthogonal
to that of mode A along $-\mathbf{\hat{r}}_{\rm{inc}}$
($\mathbf{J_{\rm{inc}}}\mathbf{J}^{\dagger}_{\rm{A}}=0$), P and A are
antipodal points
(${\ooalign{\leavevmode\resizebox{13.47224pt}{0.0pt}{{\phantom{A}}}\cr
PA}}=\pi$) and thus mode A would not be excited
[$\cos(\frac{1}{2}{\ooalign{\leavevmode\resizebox{13.47224pt}{0.0pt}{{\phantom{A}}}\cr
PA}})=0$], being consistent with the special scenario of single-mode
excitations [34, 38].

Figure 1: (a) Poincaré sphere on which P, A and B represent respectively the
incident polarization and polarizations of the mode radiations opposite to the
incident direction. The geometric phase is half solid angle of the orientated
spherical triangle $\Delta\rm{PAB}$: $\varphi_{g}=\frac{1}{2}\Omega$. (b) A
gold particle with 4-fold rotation symmetry and a spherical coordinate system
parameterized by $\mathbf{r}=(r,~{}\theta,~{}\phi)$. Scattering intensity
distributions on the $\mathbf{x}$-$\mathbf{y}$ plane: (c) P overlaps with A;
(d) P overlaps with B.

After visualizing the excitation amplitude, now we turn to the phase
$\varphi(\mathbf{\hat{r}})$ which can be generally categorized into two parts:

$\varphi(\mathbf{\hat{r}})=\varphi_{g}+\varphi_{j}(\mathbf{\hat{r}}).$ (5)

Here $\varphi_{g}$ is the global scattering direction
$\mathbf{\hat{r}}$-independent geometric phase
$\varphi_{g}=\frac{1}{2}\Omega$, where $\Omega$ denotes the solid angle
enclosed by the great-arc circuit  A PABP [see Fig. 1(a); $\Omega$ is positive
(negative) for counter-clockwise (clockwise) circuit viewed above] [29, 30,
31, 32]. The other term $\varphi_{j}(\mathbf{\hat{r}})$ is the local
$\mathbf{\hat{r}}$-dependent phase difference between the directional
radiations:
$\varphi_{j}(\mathbf{\hat{r}})=\rm{Arg}[\mathbb{J}_{\rm{A}}\mathbb{J}^{\dagger}_{\rm{B}}]$.
Since the overall phase for each QNM (and thus also their Jones vectors) can
be arbitrarily chosen, $\varphi_{j}$ seems to be indeterminate due to the
gauge freedom of adding an arbitrary extra constant phase term $\Delta\varphi$
[32]:
$\varphi_{j}(\mathbf{\hat{r}})\mapsto\varphi_{j}(\mathbf{\hat{r}})+\Delta\varphi$.
In contrast, the geometric phase term is gauge-independent and has nothing to
do with the overall assigned phase to each QNM (the position of Jones vector
on the Poincaré sphere is irrelevant to its overall phase). Such phase
indeterminacy (gauge freedom) can be fixed by the optical theorem, which
imposes a stringent constraint on both the phase and amplitude of the forward
scattering $\mathbf{E}_{\rm{sca}}(\mathbf{\hat{r}}_{\rm{inc}})$ [1].
Alternatively, $\Delta\varphi$ can be quickly decided through a single fitting
of Eq. (4) with numerically calculated results. This $\Delta\varphi$ obtained
can then be used for any incident polarizations, though generally it is
dependent on the incident direction, along which the optical theorem imposes
its constraint. For a fixed incident direction, $\Delta\varphi$ can be gauged
out by assigning to modes A and B an equal constant relative phase, and then
Eq. (1) can be simplified as (for two QNMs):

$\mathbf{E}_{\rm{sca}}(\mathbf{\hat{r}})\propto\cos(\frac{1}{2}{\ooalign{\leavevmode\resizebox{13.47224pt}{0.0pt}{{\phantom{A}}}\cr
PA}})~{}{\mathbf{E}}_{\rm{A}}(\mathbf{\hat{r}})+\cos(\frac{1}{2}{\ooalign{\leavevmode\resizebox{13.88893pt}{0.0pt}{{\phantom{A}}}\cr
PB}})~{}{\mathbf{E}}_{\rm{B}}(\mathbf{\hat{r}})\exp({i\varphi_{g}}),$ (6)

through which the distributions of not only intensities but also polarizations
can be directly calculated.

Though the geometric phase has already been widely employed for various
photonic functionalities [39, 40, 41], our theoretical framework and formulism
are profoundly different, being both deeper and broader: (i) Previous studies
are mostly based on structures consisting of optical items with spatially
varying orientations, where the geometric phase emerges among the modes of
individual items (the near-field couplings between them are thus neglected)
and the formulism is valid only for a limited range of incident directions and
polarizations. Here we discuss the optical structure (could be an individual
particle or their clusters) as a whole and the geometric phase emerges among
different QNMs of the same structure. As a result, our model has automatically
considered all sorts of couplings and our formulism can accommodate arbitrary
incident directions and polarizations; (ii) As we have revealed, what
essentially decide the geometric phase are not geometric orientations of
optical structures, but rather the polarization orientations (and also
ellipticities) of the mode radiations opposite to the incident direction. The
former is a good approximation only for some special scenarios, e.g. each item
supports solely an electric dipole along its orientation direction and the
incident direction is perpendicular to it, as is indeed the case for many
previous demonstrations. As a result, the previously widely-adopted approach
to directly associate geometric phase with the structure orientation is a
reduced approximation of our model.

We now turn to a specific scattering particle shown in Fig. 1(b) to verify our
theoretical framework (numerical calculations are performed using COMSOL
Multiphysics throughout this work). The particle consists of gold (effective
permittivity fitted from data in Ref. [42]) and exhibits four-fold rotation
symmetry that secures a pair of degenerate QNMs with central resonant angular
frequency $\omega_{1}=8.0645\times 10^{14}~{}\rm{rad}/s$. We shine plane waves
along $\mathbf{+z}$ with $\omega=\omega_{1}$ and track scattering
distributions on the $\mathbf{x}$-$\mathbf{y}$ plane, unless specified
otherwise. For other off-resonance frequencies, as long as no other modes are
excited, scattering distributions would be the same as the relative amplitude
and phase are independent of frequency detuning due to mode degeneracy. The
radiations of the two QNMs along -$\mathbf{z}$ (opposite to the incident
direction) are almost linearly polarized along $\mathbf{x}$ and $\mathbf{y}$,
with the corresponding normalized Stokes parameters [33] being respectively
$S_{1}\approx\pm 1$ [see points A and B in Fig. 2(b)]. The normalized
scattering intensity distributions on the $\mathbf{x}$-$\mathbf{y}$ plane
[parameterized by the azimuthal angle $\phi$ as shown in Fig. 1(b)] are shown
in Fig. 1(c) and 1(d), for incident linear polarizations along $\mathbf{x}$
(${\ooalign{\leavevmode\resizebox{13.47224pt}{0.0pt}{{\phantom{A}}}\cr
PA}}\approx 0$,
${\ooalign{\leavevmode\resizebox{13.88893pt}{0.0pt}{{\phantom{A}}}\cr
PB}}\approx\pi$ and thus mode B is not excited: $I_{\rm{sca}}\propto
I_{\rm{A}}$) and $\mathbf{y}$
(${\ooalign{\leavevmode\resizebox{13.88893pt}{0.0pt}{{\phantom{A}}}\cr
PB}}\approx 0$,
${\ooalign{\leavevmode\resizebox{13.47224pt}{0.0pt}{{\phantom{A}}}\cr
PA}}\approx\pi$ and thus mode A is not excited: $I_{\rm{sca}}\propto
I_{\rm{B}}$), respectively. As is clearly shown in Fig. 1(c) and 1(d),
radiations of modes A and B [$I_{\rm{A,B}}(\phi)$] agree well with the
directly simulated intensity distributions without involving QNMs, confirming
the accuracy of our model and also that by properly selecting the incident
polarizations, modes can be selectively excited as long as A and B do not
overlap.

The scattering intensity distributions for two other incident orthogonal
linear polarizations [polarized along $\phi=\pi/4$ ($S_{2}=1$) and polarized
along $\phi=3\pi/4$ ($S_{2}=-1$); see positions of P in Fig. 2(b)] are
summarized in Fig. 2(a). For both scenarios
${\ooalign{\leavevmode\resizebox{13.47224pt}{0.0pt}{{\phantom{A}}}\cr
PA}}\approx{\ooalign{\leavevmode\resizebox{13.88893pt}{0.0pt}{{\phantom{A}}}\cr
PB}}\approx\pi/2$, and the core difference is that for the former
$\varphi_{g}\approx 0$ while for the latter $\varphi_{g}\approx\pi$ [see Fig.
2(b)], leading to fully distinct scattering distributions $I_{\rm{sca}}\propto
I_{\rm{A}}+I_{\rm{B}}\pm
2\sqrt{I_{\rm{A}}I_{\rm{B}}}|\mathbb{J}_{\rm{A}}\mathbb{J}_{\rm{B}}^{\dagger}|\cos[\varphi_{j}(\mathbf{\hat{r}})]$.
The geometric phase $\pi$ is directly manifest through the two distinct curves
shown in Fig. 2(a). We note that our demonstration (with an individual
scattering particle) of such classical geometric phase $\pi$ is even simpler
and more direct than the earliest ones by Fresnel-Arago and Hamilton-Lloyd
[32]. Results for two other scenarios with fixed
${\ooalign{\leavevmode\resizebox{13.47224pt}{0.0pt}{{\phantom{A}}}\cr
PA}}\approx{\ooalign{\leavevmode\resizebox{13.88893pt}{0.0pt}{{\phantom{A}}}\cr
PB}}\approx\pi/2$ but different geometric phases [Fig. 2(b)] of
$\varphi_{g}\approx\pi/2$ [left-handed circular polarization (LCP) incidence
with $S_{3}=1$] and $\varphi_{g}\approx\pi/4$ [elliptic polarization with P
locating at ($0$, $\sqrt{2}/2$, $\sqrt{2}/2$)] are summarized in Fig. 2(c),
where the influence of geometric phase is also apparent. As is shown in Fig.
2(a), for $\varphi_{g}\approx 0$ the scattering along the direction
$\phi=\phi_{0}\approx 113^{\circ}$ is zero. Along this direction, Eq. (4) can
be reduced to $I_{\rm{sca}}(\phi=\phi_{0})\propto 1-\cos(\varphi_{g})$ as long
as ${\ooalign{\leavevmode\resizebox{13.47224pt}{0.0pt}{{\phantom{A}}}\cr
PA}}\approx{\ooalign{\leavevmode\resizebox{13.88893pt}{0.0pt}{{\phantom{A}}}\cr
PB}}\approx\pi/2$. We then show the evolutions of
$I_{\rm{sca}}(\phi=\phi_{0})$ on such a polarization great circle ($S_{1}=0$)
parameterized by $\beta\in[0,2\pi]$, as shown in the inset in Fig. 2(d).
Obviously the solid angle $\Omega=2\beta$ and $\varphi_{g}=\beta$ [Fig. 2(b)],
and thus the evolution observe the relation
$I_{\rm{sca}}(\phi=\phi_{0})\propto 1-\cos(\beta)$, which agrees perfectly
with the numerical results included in Fig. 2(d).

Figure 2: Scattering intensity distributions on the $\mathbf{x}$-$\mathbf{y}$
plane with P locating at ($0$, $1$, $0$) and ($0$, $-1$, $0$) for (a); ($0$,
$0$, $1$) and ($0$, $\sqrt{2}/2$, $\sqrt{2}/2$) for (c). The geometric phase
are respectively $\varphi_{g}=0,\pi,\pi/2$ and $\pi/4$, as shown in (b). In
(a) one diretion of zero scattering $\phi=\phi_{0}\approx 113^{\circ}$ is
specified. (d) The dependence of $I_{\rm{sca}}(\phi=\phi_{0})$ on $\beta$ that
parameterizes a great circle ($S_{1}=0$) on the Poincaré sphere (see the
inset).

Figure 3: (a) $S_{3}$ distributions on the $\mathbf{x}$-$\mathbf{y}$ plane
with P locating at ($0$, $0$, $1$) and ($0$, $1$, $0$), with
$\varphi_{g}=\pi/2$ and $0$, respeively. (b) $S_{3}$ distributions on the
whole scattering momentum sphere for P locating at ($0$, $1$, $0$), on which
the predicted directions of cirular polarization singularities are marked by
crosses. (c) Distributions of
$|\mathbf{J}_{\rm{sca}}^{\rm{LCP}}(\mathbf{J}_{\rm{sca}}^{\rm{RCP}})^{\dagger}|$
for waves incident along
$\mathbf{\hat{r}}_{\rm{inc}}=-\mathbf{\hat{r}}_{\rm{A=B}}$ (left) and
$\mathbf{\hat{r}}_{\rm{inc}}=\mathbf{+z}$ (right). (d) Variations of
$C_{\rm{sca}}$ for P transversing a great circle ($S_{2}=0$), covering matched
polarization ($\beta=0$), orthogonal polarization (marked; $\beta=\pi$), LCP
and RCP ($\beta=\pi/2$ and $3\pi/2$).

Our theoretical framework can predict the distributions of not only scattering
intensities according to Eq. (4), but also of polarizations according to Eq.
(6). The polarization distributions (in terms of $S_{3}$ of the scattered
fields) on the $\mathbf{x}$-$\mathbf{y}$ plane are shown in Fig. 3(a), for
linear polarization ($S_{2}=1$ and $\varphi_{g}\approx 0$) and LCP ($S_{3}=1$
and $\varphi_{g}\approx\pi/2$) incidences, where the agreement further
confirms the accuracy of our theoretical framework. According to Eq. (6), the
locations of circularly-polarized scatterings ($S_{3}=\pm 1$; circular
polarization singularities [13]) can be directly predicted and even designed
by selecting proper incident polarization and directions. We show the
scattering polarization distributions (simulated results) on the whole
scattering momentum sphere in Fig. 3(b) for the linear polarization incidence
($S_{2}=1$) and the predicted locations of polarization singularities are also
marked by crosses, agreeing well with the numerical calculations.

According to Eq. (6), there is a rather interesting scenario of overlapped A
and B: directions where the radiation polarizations for both QNMs are the
same, with
$\cos(\frac{1}{2}{\ooalign{\leavevmode\resizebox{13.47224pt}{0.0pt}{{\phantom{A}}}\cr
PA}})=\cos(\frac{1}{2}{\ooalign{\leavevmode\resizebox{13.88893pt}{0.0pt}{{\phantom{A}}}\cr
PB}})$ and $\varphi_{g}=0$ [see Fig. 1(a)]. For waves incident opposite to
those directions, along any scattering direction the polarization would not
change, irrespective of the incident polarization. For the marked direction
($\phi=\phi_{0}$) in Fig. 2(a) where the scattering is zero, since both modes
are excited, along this direction mode radiation polarization must be the same
or they would not cancel each other. This is also the requirement of Eq. (4),
as $I_{\rm{sca}}=0$ requires
$|\mathbb{J}_{\rm{A}}\mathbb{J}_{\rm{B}}^{\dagger}=1|$ which means an
identical radiation polarization. We denote this direction as
$\mathbf{\hat{r}}_{\rm{A=B}}$ and shine opposite to it
($\mathbf{\hat{r}}_{\rm{inc}}=-\mathbf{\hat{r}}_{\rm{A=B}}$) plane waves of
LCP and RCP (left-handed circular polarization). We then track the scattering
polarization variation on the whole momentum sphere through the parameter
$|\mathbf{J}_{\rm{sca}}^{\rm{LCP}}(\mathbf{J}_{\rm{sca}}^{\rm{RCP}})^{\dagger}|$,
where
$|\mathbf{J}_{\rm{sca}}^{\rm{LCP}}(\mathbf{J}_{\rm{sca}}^{\rm{RCP}})^{\dagger}|=1$
means the scattering polarization does not change for RCP and LCP incidences
[see Fig. 3(c)]. Here $\mathbf{J}_{\rm{sca}}$ denotes the normalized Jones
vector for the scattered expressed by Eq. (6), and its superscript indicates
the incident polarization. For comparison, we also shine circularly-polarized
waves along +$\mathbf{z}$ (where A and B do not overlap) and show its
$|\mathbf{J}_{\rm{sca}}^{\rm{LCP}}(\mathbf{J}_{\rm{sca}}^{\rm{RCP}})^{\dagger}|$
distribution in Fig. 3(c), confirming that for general incident directions the
scattering polarizations are dependent on incident polarizations.

Another interesting property for overlapped A and B is that the incident
polarization can be tuned to be orthogonal to that of both QNMs
[$\cos(\frac{1}{2}{\ooalign{\leavevmode\resizebox{13.47224pt}{0.0pt}{{\phantom{A}}}\cr
PA}})=\cos(\frac{1}{2}{\ooalign{\leavevmode\resizebox{13.88893pt}{0.0pt}{{\phantom{A}}}\cr
PB}})=0$]. For such an incident polarization, neither QNM would be excited and
thus the particle would become invisible. For the same incident direction
($\mathbf{\hat{r}}_{\rm{inc}}=-\mathbf{\hat{r}}_{\rm{A=B}}$), we show the
evolutions of scattering cross sections ($C_{\rm{sca}}$) with varying incident
polarizations in Fig. 3(d). The incident polarization covers a great circle,
covering matched polarization
($\beta={\ooalign{\leavevmode\resizebox{13.47224pt}{0.0pt}{{\phantom{A}}}\cr
PA}}={\ooalign{\leavevmode\resizebox{13.88893pt}{0.0pt}{{\phantom{A}}}\cr
PB}}=0$), LCP ($\beta=\pi/2$), orthogonal polarization
($\beta={\ooalign{\leavevmode\resizebox{13.47224pt}{0.0pt}{{\phantom{A}}}\cr
PA}}={\ooalign{\leavevmode\resizebox{13.88893pt}{0.0pt}{{\phantom{A}}}\cr
PB}}=\pi$). According to Eq. (4), scattering along any direction and thus also
the scattering cross section would adhere to
$I_{\rm{sca}},C_{\rm{sca}}~{}\propto 1+\cos(\beta)$, which is verified by Fig.
3(d). We note that the scattering evolutions in Fig. 2(d) are driven by
changing $\varphi_{g}$ with
${\ooalign{\leavevmode\resizebox{13.47224pt}{0.0pt}{{\phantom{A}}}\cr
PA}}={\ooalign{\leavevmode\resizebox{13.88893pt}{0.0pt}{{\phantom{A}}}\cr
PB}}$ fixed at $\pi/2$, while those in Fig. 3(d) driven by changing
${\ooalign{\leavevmode\resizebox{13.47224pt}{0.0pt}{{\phantom{A}}}\cr
PA}}={\ooalign{\leavevmode\resizebox{13.88893pt}{0.0pt}{{\phantom{A}}}\cr
PB}}$, with $\varphi_{g}$ fixed at $0$.

Figure 4: (a) A pair of gold split-ring resonators, and two directions of
interest ($\mathbf{\hat{r}}_{\rm{sca}}$ and $\mathbf{\hat{r}}_{\rm{inc}}$) are
also marked. The corresponding A $(-0.29,-0.41,-0.86)$ and B
$(0.38,-0.81,0.45)$ for this $\mathbf{\hat{r}}_{\rm{inc}}$ are shown in (b).
The dependence of the directional scattering intensity along
$\mathbf{\hat{r}}_{\rm{sca}}$ on incident polarizations is shown in (c), with
P transversing a great circle on the Poincaré sphere ($S_{2}=0$), covering LCP
($\beta=\pi/2$) and RCP ($\beta=3\pi/2$). (d) Scattering intensity
distributions on the $\mathbf{x}$-$\mathbf{y}$ plane with LCP and RCP
incidences along the marked incident direction.

So far we have showcased only the scenario of two degenerate QNMs while our
theoretical framework does not require such degeneracy. For non-degenerate
QNMs, deviations from the central resonant positions would only shift the
overall phase and amplitude of the corresponding QNM radiations [27], the
effect of which can be taken into consideration in our model by adjusting the
corresponding excitation coefficients accordingly, with the geometric phase
term unaffected. We have further studied a scattering structure shown in Fig.
4(a), which supports a pair of non-degenerate QNMs with central resonant
frequencies close to $\omega_{2}=1.1695\times 10^{15}~{}\rm{rad}/s$, at which
the angular frequency of the incident wave is fixed ($\omega=\omega_{2}$). We
have identified a scattering direction $\mathbf{\hat{r}}_{\rm{sca}}$
($\phi=41^{\circ}$ and $\theta=90^{\circ}$; on the $\mathbf{x}$-$\mathbf{y}$
plane) along which the radiation polarizations of both QNMs are identical. As
has already been elaborated, through tuning the relative exciting amplitude
and phase, scattering along this direction can be eliminated, similar to that
marked in Fig. 2(a). We have selected an incident direction
$\mathbf{\hat{r}}_{\rm{inc}}$ ($\phi=0$ and $\theta=10^{\circ}$) and the
corresponding positions of A and B are shown in Fig. 4(b). The evolution of
the directional scattering intensity along $\mathbf{\hat{r}}_{\rm{sca}}$ with
changing incident polarizations covering a great circle ($S_{2}=0$) is shown
in Fig. 4(c), where the scattering is fully suppressed for LCP incidence
($\beta=\pi/2$). We have also tracked the scattering intensity distributions
on the $\mathbf{x}$-$\mathbf{y}$ plane for LCP and RCP incidences, with the
results summarized in Fig. 4(d).

In conclusion, we have unveiled the hidden excitation geometric phases among
QNMs of scatterers and reveal how they drive scattering evolutions on the
incident Poincaré sphere. The geometric phase can be exploited to efficiently
manipulate the scatterings, such as scattering eliminations and polarization
singularity (zero and circularly-polarized scatterings) generations. For the
general scenario of more than two QNMs being simultaneously excited, the
relative amplitude and phase among any two QNMs can be calculated using our
model and then calculations for interferences among all QNMs become routine.
This means that the theoretical framework we have constructed is generic and
widely applicable. We have essentially merged three vibrant fields of
geometric phase, Mie scattering and QNM, and unlocked an extra hidden
dimension of electromagnetic scattering. Similar dimensions might be uncovered
for waves of other forms for which geometric phase is generic and ubiquitous,
providing new flexibilities for many physics and interdisciplinary branches
that are related to wave scatterings.

_Acknowledgment_ : This work is supported by National Natural Science
Foundation of China (Grant No. 11874426 and 61405067), and several Researcher
Schemes of National University of Defense Technology. W. L. acknowledges many
illuminating correspondences with Sir Michael Berry, whose monumental paper on
geometric phase was published $40$ years ago today.

## References

  * [1] C. F. Bohren and D. R. Huffman, _Absorption and Scattering of Light by Small Particles_ (Wiley, 1983).
  * [2] A. I. Kuznetsov, A. E. Miroshnichenko, M. L. Brongersma, Y. S. Kivshar, and B. Luk'yanchuk, ``Optically resonant dielectric nanostructures,'' Science 354, aag2472 (2016).
  * [3] W. Liu and Y. S. Kivshar, ``Generalized Kerker effects in nanophotonics and meta-optics [Invited],'' Opt. Express 26, 13085–13105 (2018).
  * [4] A. Krasnok, D. Baranov, H. Li, M.-A. Miri, F. Monticone, and A. Alú, ``Anomalies in light scattering,'' Adv. Opt. Photon. 11, 892–951 (2019).
  * [5] Y. Kivshar, ``The Rise of Mie-tronics,'' Nano Lett. 22, 3513–3515 (2022).
  * [6] H. Altug, S.-H. Oh, S. A. Maier, and J. Homola, ``Advances and applications of nanophotonic biosensors,'' Nat. Nanotechnol. 17, 5–16 (2022).
  * [7] S. Fan and W. Li, ``Photonics and thermodynamics concepts in radiative cooling,'' Nat. Photon. 16, 182–190 (2022).
  * [8] S. Pancharatnam, ``The propagation of light in absorbing biaxial crystals,'' Proc. Indian. Acad. Sci. 42, 86–109 (1955).
  * [9] M. Berry, ``Pancharatnam, virtuoso of the Poincaré sphere: An appreciation,'' Curr. Sci. 67, 220–220 (1994).
  * [10] T. Ozawa, H. M. Price, A. Amo, N. Goldman, M. Hafezi, L. Lu, M. C. Rechtsman, D. Schuster, J. Simon, O. Zilberberg, and I. Carusotto, ``Topological photonics,'' Rev. Mod. Phys. 91, 015006 (2019).
  * [11] C. Wang, Z. Fu, W. Mao, J. Qie, A. D. Stone, and L. Yang, ``Non-Hermitian optics and photonics: From classical to quantum,'' Adv. Opt. Photon., AOP 15, 442–523 (2023).
  * [12] C. W. Hsu, B. Zhen, A. D. Stone, J. D. Joannopoulos, and M. Soljačić, ``Bound states in the continuum,'' Nat. Rev. Mater. 1, 16048 (2016).
  * [13] G. J. Gbur, _Singular Optics_ (CRC Press Inc, Boca Raton, 2016).
  * [14] W. Liu, W. Liu, L. Shi, and Y. Kivshar, ``Topological polarization singularities in metaphotonics,'' Nanophotonics 10, 1469–1486 (2021).
  * [15] M.-A. Miri and A. Alù, ``Exceptional points in optics and photonics,'' Science 363, eaar7709 (2019).
  * [16] K. Koshelev, A. Bogdanov, and Y. Kivshar, ``Meta-optics and bound states in the continuum,'' Science Bulletin 64, 836–842 (2019).
  * [17] M. Kang, T. Liu, C. T. Chan, and M. Xiao, ``Applications of bound states in the continuum in photonics,'' Nat Rev Phys 5, 659–678 (2023).
  * [18] A. Doicu, T. Wriedt, and Y. A. Eremin, _Light Scattering by Systems of Particles: Null-Field Method with Discrete Sources: Theory and Programs_ , vol. 124 (Springer, 2006).
  * [19] T. Needham, _Visual Differential Geometry and Forms: A Mathematical Drama in Five Acts_ (Princeton University Press, Princeton, 2021).
  * [20] W. Liu, J. Zhang, B. Lei, H. Ma, W. Xie, and H. Hu, ``Ultra-directional forward scattering by individual core-shell nanoparticles,'' Opt. Express 22, 16178 (2014).
  * [21] J. D. Jackson, _Classical Electrodynamics Third Edition_ (Wiley, New York, 1998), 3rd ed.
  * [22] W. Liu, ``Generalized magnetic mirrors,'' Phys. Rev. Lett. 119, 123902 (2017).
  * [23] W. Chen, Y. Chen, and W. Liu, ``Singularities and Poincaré indices of electromagnetic multipoles,'' Phys. Rev. Lett. 122, 153907 (2019).
  * [24] W. Chen, Y. Chen, and W. Liu, ``Line Singularities and Hopf Indices of Electromagnetic Multipoles,'' Laser Photonics Rev. 14, 2000049 (2020).
  * [25] I. Fernandez-Corbaton, X. Zambrana-Puyalto, N. Tischler, X. Vidal, M. L. Juan, and G. Molina-Terriza, ``Electromagnetic Duality Symmetry and Helicity Conservation for the Macroscopic Maxwell's Equations,'' Phys. Rev. Lett. 111, 060401 (2013).
  * [26] Q. Yang, W. Chen, Y. Chen, and W. Liu, ``Electromagnetic Duality Protected Scattering Properties of Nonmagnetic Particles,'' ACS Photonics 7, 1830–1838 (2020).
  * [27] P. Lalanne, W. Yan, K. Vynck, C. Sauvan, and J.-P. Hugonin, ``Light Interaction with Photonic and Plasmonic Resonances,'' Laser Photonics Rev. 12, 1700113 (2018).
  * [28] C. Caloz, A. Alù, S. Tretyakov, D. Sounas, K. Achouri, and Z.-L. Deck-Léger, ``Electromagnetic Nonreciprocity,'' Phys. Rev. Applied 10, 047001 (2018).
  * [29] S. Pancharatnam, ``Generalized theory of interference, and its applications,'' Proc. Indian. Acad. Sci. 44, 247–262 (1956).
  * [30] M. V. Berry, ``Quantal Phase Factors Accompanying Adiabatic Changes,'' Proc. R. Soc. A 392, 45–57 (1984).
  * [31] M. V. Berry, ``The Adiabatic Phase and Pancharatnam Phase for Polarized-Light,'' J. Mod. Opt. 34, 1401 (1987).
  * [32] M. Berry, ``A Geometric-Phase Timeline,'' Optics & Photonics News 3, 42–49 (2024).
  * [33] A. Yariv and P. Yeh, _Photonics: Optical Electronics in Modern Communications_ (Oxford University Press, New York, 2006), 6th ed.
  * [34] W. Chen, Q. Yang, Y. Chen, and W. Liu, ``Extremize Optical Chiralities through Polarization Singularities,'' Phys. Rev. Lett. 126, 253901 (2021).
  * [35] G. N. Ramachandran and S. Ramaseshan, ``Crystal Optics,'' in ``Kristalloptik $\cdot$ Beugung/Crystal Optics $\cdot$ Diffraction,'' , S. Flügge, ed. (Springer, Berlin, Heidelberg, 1961), Handbuch Der Physik / Encyclopedia of Physics, pp. 1–217.
  * [36] S. Pancharatnam, _Collected Works of S. Pancharatnam_ (Oxford University Press for the Raman Research Institute, London, 1975).
  * [37] M. V. Berry, ``The adiabatic phase and pancharatnam phase for polarized-light,'' J. Mod. Opt. 34, 1401 (1987).
  * [38] C. Wen, J. Zhang, S. Qin, Z. Zhu, and W. Liu, ``Momentum-Space Scattering Extremizations,'' Laser & Photonics Reviews 18, 2300454 (2024).
  * [39] G. Li, S. Zhang, and T. Zentgraf, ``Nonlinear photonic metasurfaces,'' Nat. Rev. Mater. 2, 17010 (2017).
  * [40] E. Cohen, H. Larocque, F. Bouchard, F. Nejadsattari, Y. Gefen, and E. Karimi, ``Geometric phase from Aharonov–Bohm to Pancharatnam–Berry and beyond,'' Nat. Rev. Phys. 1, 437–449 (2019).
  * [41] C. P. Jisha, S. Nolte, and A. Alberucci, ``Geometric Phase in Optics: From Wavefront Manipulation to Waveguiding,'' Laser & Photonics Reviews 15, 2100003 (2021).
  * [42] P. B. Johnson and R. W. Christy, ``Optical constants of the noble metals,'' Phys. Rev. B 6, 4370 (1972).


# Neural Bayesian Network Understudy

Paloma Rabaey, Cedric De Boom, Thomas Demeester
IDLab, Dept. of Information Technology
Ghent University - imec
Ghent, Belgium
<EMAIL_ADDRESS>

###### Abstract

Bayesian Networks may be appealing for clinical decision-making due to their
inclusion of causal knowledge, but their practical adoption remains limited as
a result of their inability to deal with unstructured data. While neural
networks do not have this limitation, they are not interpretable and are
inherently unable to deal with causal structure in the input space. Our goal
is to build neural networks that combine the advantages of both approaches.
Motivated by the perspective to inject causal knowledge while training such
neural networks, this work presents initial steps in that direction. We
demonstrate how a neural network can be trained to output conditional
probabilities, providing approximately the same functionality as a Bayesian
Network. Additionally, we propose two training strategies that allow encoding
the independence relations inferred from a given causal structure into the
neural network. We present initial results in a proof-of-concept setting,
showing that the neural model acts as an understudy to its Bayesian Network
counterpart, approximating its probabilistic and causal properties.

##  1 Introduction

Bayesian Networks (BNs) have been the subject of a large body of research, and
are considered an important technology for high-stakes decision making, in
particular in the healthcare domain. Kyrimi et al. [2021b] provide a clear
overview of the key properties that make BNs highly appealing in that context:
(i) their ability to model complex problems involving causal dependencies as
well as uncertainty, (ii) their ability to combine data and expert knowledge,
(iii) their interpretable graphical structure, and (iv) their ability to model
interventions and reason diagnostically as well as prognostically. Despite
efforts from the research community to support healthcare practitioners in
using BNs [Arora et al., 2019, Kyrimi et al., 2020], their practical adoption
remains very limited. The analysis of Kyrimi et al. [2021a] indicates data
inadequacies as one of several key limiting factors: BNs have a limited
capacity to address continuous data, and cannot directly deal with
unstructured data such as images or free text. According to Tayefi et al.
[2021], the value of information entailed by unstructured text fields in
electronic medical records for clinical decision support cannot be
underestimated.

In contrast to Bayesian networks, neural networks are well-suited for dealing
with unstructured data. However, they are inherently not interpretable, which
stands in the way of clinical adoption [Montani and Striani, 2019, Peiffer-
Smadja et al., 2020]. Whereas Hinton [2018] argues strongly in favour of the
adoption of neural networks in healthcare, he also recognises interpretability
as a property that is desired by practitioners.

Our envisioned research goal is to build models that combine the advantages of
both approaches: a fully neural model that can be trained on observations
including unstructured data, but that also adopts the desirable properties of
BNs as listed above. Realistically speaking, we can at best expect only
approximate probabilistic reasoning capabilities from a generic neural
network, and it will lack the direct interpretability of a BN. However, if the
model displays the causal properties of a BN, the ideas of counterfactual
explanation [Wachter et al., 2017] can be applied to explain its predictions.

This paper forms the first step along that outlined research track. We focus
on the problem of training a neural network on a given set of discrete data
samples as well as knowledge of a Directed Acyclic Graph (DAG) [Russell and
Norvig, 2009] reflecting causal relationships between the variables. We do not
yet aim to augment such models with unstructured data, but for now require
that it behaves as a BN with the given DAG, if we were to train it on the same
samples. In this sense, our model can be considered a _neural understudy_ to
the BN. We want to emphasise that our current preliminary work does not aim to
build a neural network which outperforms its Bayesian network counterpart, but
one that approximates its capabilities. After this necessary first step, we
expect to reap the benefits of the neural network understudy by incorporating
unstructured data. While the current work presents results only in a proof-of-
concept setting, the possibility to reason under uncertainty, while combining
causal structure and unstructured data, has great potential in real-world
(clinical) applications.

In short, this paper makes the following contributions: (1) we propose a
neural network model able to infer the probability of a set of target
variables conditioned on observed evidence, (2) we introduce two approaches
for injecting DAG knowledge while training the neural network, and (3) we
present empirical results (including robustness analysis) in a small proof-of-
concept setting. The necessary code to reproduce these results can be found in
our Github repository.111https://github.com/prabaey/NBN-understudy

##  2 Related work

This section focuses on recent works most related to the proposed ideas; we do
not aim to provide a broad overview in terms of related research fields. A
basic premise for the presented work is the existence of causal knowledge in
addition to observational data, to be injected into a neural network model.
Our goal is therefore not _causal discovery_ , intended to recover the causal
mechanisms underlying the distribution from which the observed data was
generated [Ke et al., 2019, Brouillard et al., 2020].

Geiger et al. [2022] present a method to align different parts of a neural
network with nodes in a causal graph, resulting in a model with improved
interpretability. Harradon et al. [2018] pursue a similar goal, summarising a
neural model into a Bayesian causal model, to provide counterfactual
explanations for the neural model. Instead of starting from a neural network
and using causal models to understand or guide its inner workings, we work the
other way round.

In Rohekar et al. [2018], conditional independencies in the input distribution
are encoded hierarchically in the network structure, as a way of performing
unsupervised neural network structure learning. Pawlowski et al. [2020] and
Parafita and Vitrià [2020] present each parent-child relation in a given
structural causal model as a deep neural network unit. While we don’t
explicitly align parts of the neural network with a causal model, we do
include structural information in the training process as a regulariser. This
is related to the work from Kyono et al. [2020], who propose a prediction task
with a regularisation objective to detect causal relationships between
features. It is not based on prior knowledge of causal relationships, as in
our work.

Attempts have been made to build neural networks that implement belief
propagation [Pearl, 1988] to improve on accuracy and efficiency [Garcia
Satorras and Welling, 2021, Kuck et al., 2020]. Our proposed neural networks
also aim at answering probabilistic queries, but can freely learn the
inference mechanics from the data.

Another key related contribution is DeepProbLog [Manhaeve et al., 2018]. This
neuro-symbolic model is able to implement our envisioned goal: its
probabilistic logic programming core can represent any BN exactly, and its so-
called neural predicates are able to encode unstructured data. However,
whereas a fully neural model can at best reason only approximately with
probabilities, it may have some benefits over the exact probabilistic
reasoning component in DeepProbLog. Indeed, in future research we aim to model
complex realistic observational data. Some of its properties a suitable model
should be able to capture, may be hard to express in a probabilistic program.

Finally, our work is not to be confused with the so-called _Bayesian neural
networks_ MacKay [1992], aiming to obtain a probability distribution over
neural networks. Instead, our goal is training a neural understudy for a
Bayesian Network, i.e., a neural model that behaves similarly.

##  3 Training a neural understudy of a Bayesian network

The following paragraphs describe our proposed neural architecture with full
flexibility in the choice of input and output variables (Section 3.1),
followed by our proposed training strategies to inject causal structure into
the model (Section 3.2).

###  3.1 Neural architecture and training

A fully specified BN can be queried by providing inputs to any selection of
variables (henceforth called the _evidence_), after which the probabilities
for each of the remaining variables (the _targets_) can be inferred,
conditioned on the evidence. We will require the same ability from our neural
counterpart, and show that this can be achieved by using complementary masks
at its inputs and outputs. This only induces restrictions on the input and
output layers but not on the internal neural architecture.

We focus on problems involving a set of $N$ discrete variables
$\mathcal{V}=\\{X_{1},X_{2},...,X_{N}\\}$, in which $X_{i}$ can take $n_{i}$
possible values in $\mathcal{X}_{i}=\\{x_{i1},...,x_{in_{i}}\\}$. Given a set
$\mathcal{E}\subset\mathcal{V}$ of $m$ evidence variables ($m<N$), each with
an assigned value $\\{E_{1}=e_{1},...,E_{m}=e_{m}\\}$ (for convenience written
as $\mathcal{E}=e$), our goal is to train a neural network to predict the
probability distribution $P(X|\mathcal{E}=e)$ for each of the remaining
variables $X$. The latter are called the target variables, aggregated in the
set $\mathcal{T}=\mathcal{V}\setminus\mathcal{E}$.

After training on an observed set of individual samples of the form
$\\{X_{1}=x_{1},\ldots,X_{N}=x_{N}\\}$, our model should be able to deal with
any selection of evidence variables. The key idea to achieve that is through
dynamic masking (i.e., a different mask for any variable split $\mathcal{E}$
vs. $\mathcal{T}$): all variables are present at the model input and output,
but target variables are filtered out of the inputs through an evidence mask,
whereas predictions for the evidence variables are filtered from the outputs
through a target mask. This is illustrated in Fig. 1.

Figure 1: Illustration: training with evidence $\mathcal{E}=\\{X_{1},X_{2}\\}$
and targets $\mathcal{T}=\\{X_{3}\\}$, for
$P(X_{3}|X_{1}=x_{12},X_{2}=x_{21})$. The mask at the input (output) indicates
selection of evidence (target) entries.

#### Model architecture:

Each input to the model is a vector $v_{\text{in}}$ of dimension
$k=\sum_{i=1}^{N}n_{i}$, and is formed by concatenating the $n_{i}$-sized one-
hot representation of the value held by each of the $N$ variables $X_{i}$
($i=1,\ldots,N$). The masked positions of the target variables in
$v_{\text{in}}$ are substituted by the corresponding components of a static
vector $v_{0}$. The input $v_{\text{in}}$ is subsequently passed through one
or more hidden layers, resulting at the output in a $k$-dimensional vector of
logits. These are then locally normalised through softmax activations, for
each variable $X_{i}$ considering its corresponding $n_{i}$ entries. The
output vector $v_{\text{out}}$ is finally constructed by replacing the entries
at the positions of evidence variables, by the corresponding one-hot
representations from $v_{\text{in}}$. The values corresponding to the target
variables can be interpreted as the predicted target distributions given the
evidence. This is illustrated in Fig. 1. Note that the $k$-dimensional vector
$v_{0}$ is obtained by applying the linear output layer with per-variable
softmax normalisation on a randomly initialised trainable vector.

#### Model training:

Whereas a BN is naturally able to answer any query once its conditional
probability tables are specified, we need to explicitly train our neural
network for that ability. To that end, we iterate over the observed samples,
randomly dividing the variables into evidence set $\mathcal{E}$ and targets
$\mathcal{T}$, and minimise the average cross-entropy loss for the target
variables. The training loss can be expressed as
$\mathcal{L}^{T}=\frac{1}{|\mathcal{T}|}\sum_{X_{i}\in\mathcal{T}}-\log\hat{p}_{ij}$,
in which $\hat{p}_{ij}$ denotes the entry in $v_{\text{out}}$ corresponding to
the actual observed input value $x_{ij}$ of target variable $X_{i}$, i.e., the
predicted output probability for the target class of $X_{i}$.

Even though the model only receives discrete data samples and therefore never
observes the class probabilities it is meant to predict, this training
strategy steers the predicted probabilities towards the empirical
(conditional) class probabilities. It only requires that each sample’s
frequency of occurrence during training corresponds to its relative frequency
in the set of observations. Appendix A.1 proves this for the canonical case of
two binary variables, but the presented derivation can be extended to the
multi-variable case.

###  3.2 Training with causal structure

We now propose the following two training strategies to augment the neural
network with knowledge of the causal structure between the variables. They are
based on the assumption that the causal knowledge can be represented as a DAG,
reflecting the (conditional) independence relations that exist between the
variables.

#### 1\. Injecting independence relations through regularisation (REG):

In general, we can say a given DAG implies $M$ pairwise conditional
independence relations between two variables, conditioned on a set of observed
variables, shortly written as $X\perp Y\mid\mathcal{A}$, with
$X,Y\in\mathcal{V}$ and the conditioning set
$\mathcal{A}\subset\mathcal{V}\setminus\\{X,Y\\}$. The first proposed training
strategy is to inject these independence relations by constructing
regularisation loss terms $\mathcal{L}_{X\perp Y\mid\mathcal{A}}^{\text{REG}}$
that quantify how strongly the model violates them:

$\mathcal{L}_{X\perp
Y\mid\mathcal{A}}^{\text{REG}}=\frac{1}{n}\sum_{j=1}^{n}\big{(}\hat{p}(X=x_{j}\mid\mathcal{A}=a,Y=y)-\hat{p}(X=x_{j}\mid\mathcal{A}=a,Y=y^{\prime})\big{)}^{2}$
(1)

in which the summation runs over all $n$ possible values of $X$. The model’s
predicted probability for $X=x_{j}$ with as evidence the assignment
$\mathcal{A}=a$ and $Y=y$, is denoted as
$\hat{p}(X=x_{j}\mid\mathcal{A}=a,Y=y)$. The loss term expresses that if the
independence relation were satisfied by the model, its predicted probability
for any value of $X$, given any evidence $\mathcal{A}=a$, should be
independent of the value assumed by $Y$. The conditioning set $\mathcal{A}$ is
randomly instantiated every time the corresponding loss term applies during
training, and the values $y$ and $y^{\prime}$ of $Y$ are chosen randomly (with
$y\not=y^{\prime}$).

During training, the regular loss $\mathcal{L}^{T}$ per data sample is
augmented with the regularisation loss $\mathcal{L}^{\text{REG}}$ of one
sampled independence relation, weighted with a hyperparameter $\alpha$.

#### 2\. Injecting independence relations through evidence corruption (COR):

The second training strategy is based on the intuition that for particular
observed data samples, the value of some of the evidence variables may no
longer matter besides the other observed variables, when accounting for the
relevant independence relations. During training, we detect these cases, and
randomly corrupt (i.e., re-sample) those values in the evidence presented to
the network.

Consider training on a particular observed sample, the variables divided into
evidence $\mathcal{E}$ and targets $\mathcal{T}$. We then go through the known
independence relations $X\perp Y\mid\mathcal{A}$, to see which ones are
relevant to the training instance. This is the case if either of the following
two conditions hold: (1) $(\\{Y\\}\cup\mathcal{A})=\mathcal{E}$ and
$X\subset\mathcal{T}$, or (2) $(\\{X\\}\cup\mathcal{A})=\mathcal{E}$ and
$Y\subset\mathcal{T}$. When condition (1) holds, the predicted outcome of
target $X$ should not depend on the observed value of $Y$, as prescribed by
the independence relation, since the conditioning set that is needed for this
relation to hold is indeed part of the evidence. The observed input value for
$Y$ can hence be disregarded, and we randomly assign a new value from its
possible classes. Such a corruption of the input will condition the model to
ignore $Y$ when predicting $X$, given that all variables in $\mathcal{A}$ are
also provided as evidence. Note that predictions for other variables in
$\mathcal{T}$ need to be done based on the original (i.e., non-corrupted)
value for $Y$, as there is no guarantee that the same independence relation
holds for those targets as well. We apply a similar reasoning when condition
(2) holds, now corrupting the value of $X$ instead of $Y$. If no relevant
independence relation is found for a given selection of evidence and targets,
the input sample is not corrupted and passed as-is to the model. In Appendix
A.2, we consider the basic two-variable case and show that the corruption
strategy leads to the desired predicted probabilities.

##  4 Neural understudy of Bayesian network: proof-of-concept

The experiments presented in this section aim at answering the following
research questions:

  * -

RQ1: How does the basic neural network perform in comparison with a Bayesian
Network trained on the same set of samples, in terms of accuracy in predicted
probabilities and in terms of sample efficiency? (Section 4.2)

  * -

RQ2: What is the effect of injecting causal structure knowledge into our
neural network on its performance? (Section 4.2)

  * -

RQ3: How does each model’s performance change when faced with a partially
incorrect specification of the causal structure? (Section 4.3)

To this end, we compare the following models, trained on an artificial dataset
of samples generated from a given joint distribution, as specified in the next
section:

  * -

Bayesian network baseline (BN): Bayesian Network with the correct DAG, its
joint probability distribution estimated from the training samples using
Maximum Likelihood Estimation with a K2 prior (see Appendix A.3).

  * -

Neural network baseline (NN): Basic neural network approach as presented in
Section 3.1, implemented in a single-layer feed-forward model (see Appendix
A.4 for details).

  * -

NN with DAG-based regularisation (NN+REG): This model extends the basic NN
model with independence relation information, through the regularisation
technique REG (Section 3.2).

  * -

NN with DAG-based corruption (NN+COR): This model, also based on NN, injects
independence relations during training with the corrupted inputs strategy
introduced in Section 3.2.

###  4.1 Evaluation

Figure 2: The “Asia” BN model, describing a basic lung cancer detection system
(for details, see Appendix A.3).

For the experiments in this section, we make use of the “Asia” network based
on [Lauritzen and Spiegelhalter, 1988], as depicted in Fig. 2. Despite its
simplicity, with only 7 binary variables, it has sufficient connectivity
between the variables to make the problem non-trivial. Following the rules of
D-separation [Verma and Pearl, 1990], we can in total extract 191 unique
independence relations from its DAG, for the training strategies NN+REG and
NN+COR.

From the ground-truth “Asia” model (detailed in Appendix A.3), we can draw
data samples to use as training instances. Through variable elimination
[Russell and Norvig, 2009], we can build test queries with the ground-truth
conditional probabilities of target variables, for any particular assignment
of evidence variables. Just like during training, this means that one query
may contain multiple target variables. We can then evaluate models by
iterating through the test queries, and calculating the mean absolute error
(MAE) between the predicted target distribution and the desired one. For each
query, we compute the MAE per target variable, sum all their contributions and
divide by the number of targets in the query. We use two different test sets.
The first set contains all possible evidence assignments, in total 2,059
queries, and leads to the total MAE. The second set contains 1,000 queries,
each constructed by drawing a sample from the joint probability distribution,
then randomly selecting the set of evidence variables, and obtaining the
conditional probabilities for the target variables. It measures the sample
MAE, assigning more weight to the model’s performance for queries with
commonly observed evidence values, while the total MAE puts equal focus on
common as well as rare assignments of the evidence variables.

###  4.2 Performance of Bayesian network vs. neural understudy

Figure 3: Comparison of 4 models in terms of total MAE and sample MAE for
different sizes of the training set. The lines reflect the mean across 10
random seeds, used to sample the training set and initialise the neural
networks. The shading represents the 95% confidence interval.

In Fig. 3 we compare the sample efficiency of our 4 models in terms of total
and sample MAE. For the NN+REG model, we set the regularisation parameter
$\alpha$ to 10. The training hyperparameters (hidden layer dimensions, batch
size, learning rate, number of training epochs) are kept fixed for all neural
models to allow for a direct comparison. Details on the training process and
hyperparameters are given in Appendix A.4. The results allow us to answer the
first two research questions.

#### RQ1 (Sample efficiency)

The proposed training strategy leads to neural understudy models that can
approximately infer conditional probabilities. For smaller training sets, the
basic NN model shows similar performance as the BN baseline in terms of total
MAE. As the training set grows larger, the BN model’s knowledge of the DAG
structure allows it to significantly overtake the basic NN model in terms of
performance. The BN baseline outperforms the neural models across the board
when looking at sample MAE, though the gap is not large.

#### RQ2 (Effect of including DAG info)

Including information on the DAG structure in the neural model improves its
performance, with a more visible effect on the total MAE. In terms of that
metric, injecting the independence relations during training allows the neural
model to outperform the BN baseline for smaller sample sets. The NN+REG model
performs slightly better than the NN+COR model, albeit not significantly. We
hypothesise that this stems from the REG technique systematically iterating
over all conditional independence relations, whereas the COR technique only
applies a relation when the randomly sampled evidence forms a match. The
strength of the neural models appears to lie in their improved performance for
rare evidence combinations, which are more heavily disadvantaged in smaller
datasets.

###  4.3 Robustness against miss-specification of causal structure

Figure 4: Boxplots visualising the total MAE for models receiving the correct
DAG (“base”) vs. a partially miss-specified DAG (“miss”), i.e., either one
edge removed (left) or added (right), for different models (BN; NN+REG;
NN+COR) trained on 100 observed samples. All train runs are done for 5 random
seeds, and for “miss” runs, 5 random DAG corruptions are done (resulting in 25
runs per model for the “miss” setting).

We now explore the impact of injecting partially incorrect information on the
causal structure. To this end, we create 10 new DAGs, by randomly removing or
adding 222We were careful not to add edges which might introduce cycles in the
network. one edge at a time to the ground-truth DAG from Fig. 2. While still
being trained on samples from the correct “Asia” network, the BN, NN+REG and
NN+COR models now get their conditional independence relations from a
partially incorrect DAG. The results displayed in Fig. 4 allow answering the
third research question.

#### RQ3 (Partially incorrect DAG)

While the inclusion of DAG information in the neural understudy brings clear
additional benefits in terms of total MAE compared to its most basic form, we
must be careful to correctly specify this causal structure. When we assume two
variables to be independent when they are not (i.e. by removing an edge in the
DAG), both the Bayesian baseline and the neural models become less stable,
showing higher variation across sample sets. Adding an edge in the DAG seems
to do no harm, since this makes for a more conservative estimate of the DAG
structure (we assume two variables to have some dependency when they do not).

##  5 Conclusions and future work

We presented ideas for building neural networks that behave like Bayesian
Networks. In future research we aim to combine the benefits of neural networks
(i.e., encoding unstructured data) with some advantages of Bayesian networks,
such as their ability to combine uncertainty with causal structure. As a first
contribution in that direction, we presented a method to learn a neural
network model on observational sample data, and two strategies to encode known
causal relationships between the variables, by injecting independence
relations.

We tested our method on a single small-scale example, and saw that our
proposed training strategy generally works: our neural understudy models
(w.r.t. to their Bayesian Network counterpart trained on the same samples)
were able to make approximate predictions of conditional probabilities. The
inclusion of causal structure resulted in similar performance of the neural
understudy compared to the BN baseline, with the neural understudy slightly
outperforming the BN under low-sample regimes when testing on all possible
queries. We saw the performance of the NN models become less stable as soon as
incorrect independence relations were injected, though this behaviour was also
observed for the BN.

We see multiple avenues for future work. We will first extend our experiments
to larger, more realistic settings, to see whether the conclusions based on
our small example still hold. We then aim to investigate how our models can be
extended towards continuous variables as well as unstructured data, to be able
to answer probabilistic queries concerning a combination of any of these
inputs. For example, we envision the use of pre-trained language encoders,
although that will require considerable adaptation of the straightforward
model with corresponding discrete nodes at input and output.

## References

  * Ankur Ankan and Abinash Panda [2015] Ankur Ankan and Abinash Panda.  pgmpy: Probabilistic Graphical Models using Python.  In _Proceedings of the 14th Python in Science Conference_ , pages 6 – 11, 2015.
  * Arora et al. [2019] Paul Arora, Devon Boyne, Justin J. Slater, Alind Gupta, Darren R. Brenner, and Marek J. Druzdzel.  Bayesian networks for risk prediction using real-world data: A tool for precision medicine.  _Value in Health_ , 22(4), 2019.
  * Brouillard et al. [2020] Philippe Brouillard, Sébastien Lachapelle, Alexandre Lacoste, Simon Lacoste-Julien, and Alexandre Drouin.  Differentiable causal discovery from interventional data.  In _Advances in Neural Information Processing Systems_ , volume 33, pages 21865–21877. Curran Associates, Inc., 2020.
  * Garcia Satorras and Welling [2021] Víctor Garcia Satorras and Max Welling.  Neural enhanced belief propagation on factor graphs.  In _Proceedings of The 24th International Conference on Artificial Intelligence and Statistics_ , volume 130 of _Proceedings of Machine Learning Research_ , pages 685–693. PMLR, 2021.
  * Geiger et al. [2022] Atticus Geiger, Zhengxuan Wu, Hanson Lu, Josh Rozner, Elisa Kreiss, Thomas Icard, Noah Goodman, and Christopher Potts.  Inducing causal structure for interpretable neural networks.  In _Proceedings of the 39th International Conference on Machine Learning_ , volume 162 of _Proceedings of Machine Learning Research_ , pages 7324–7338. PMLR, 2022.
  * Harradon et al. [2018] Michael Harradon, Jeff Druce, and Brian E. Ruttenberg.  Causal learning and explanation of deep neural networks via autoencoded activations, 2018.  URL http://arxiv.org/abs/1802.00541.
  * Hinton [2018] Geoffrey Hinton.  Deep Learning—A Technology With the Potential to Transform Health Care.  _Journal of the American Medical Association (JAMA)_ , 320(11):1101–1102, 2018.
  * Ke et al. [2019] Nan Rosemary Ke, Olexa Bilaniuk, Anirudh Goyal, Stefan Bauer, Hugo Larochelle, Bernhard Schölkopf, Michael C. Mozer, Chris Pal, and Yoshua Bengio.  Learning neural causal models from unknown interventions, 2019.  URL https://arxiv.org/abs/1910.01075.
  * Kingma and Ba [2015] Diederik P. Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  In _3rd International Conference on Learning Representations, ICLR_ , 2015.
  * Kuck et al. [2020] Jonathan Kuck, Shuvam Chakraborty, Hao Tang, Rachel Luo, Jiaming Song, Ashish Sabharwal, and Stefano Ermon.  Belief propagation neural networks.  In _Advances in Neural Information Processing Systems_ , volume 33, pages 667–678. Curran Associates, Inc., 2020.
  * Kyono et al. [2020] Trent Kyono, Yao Zhang, and Mihaela van der Schaar.  Castle: Regularization via auxiliary causal graph discovery.  In _Advances in Neural Information Processing Systems_ , volume 33, pages 1501–1512. Curran Associates, Inc., 2020.
  * Kyrimi et al. [2020] Evangelia Kyrimi, Mariana Raniere Neves, Scott McLachlan, Martin Neil, William Marsh, and Norman Fenton.  Medical idioms for clinical bayesian network development.  _Journal of Biomedical Informatics_ , 108, 2020.
  * Kyrimi et al. [2021a] Evangelia Kyrimi, Kudakwashe Dube, Norman Fenton, Ali Fahmi, Mariana Raniere Neves, William Marsh, and Scott McLachlan.  Bayesian networks in healthcare: What is preventing their adoption?  _Artificial Intelligence in Medicine_ , 116, 2021a.
  * Kyrimi et al. [2021b] Evangelia Kyrimi, Scott McLachlan, Kudakwashe Dube, Mariana R. Neves, Ali Fahmi, and Norman Fenton.  A comprehensive scoping review of bayesian networks in healthcare: Past, present and future.  _Artificial Intelligence in Medicine_ , 117, 2021b.
  * Lauritzen and Spiegelhalter [1988] S. L. Lauritzen and D. J. Spiegelhalter.  Local Computation with Probabilities on Graphical Structures and their Application to Expert Systems (with discussion).  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 50(2):157–224, 1988.
  * MacKay [1992] David J. C. MacKay.  A practical bayesian framework for backpropagation networks.  _Neural Computation_ , 4(3):448–472, 1992.
  * Manhaeve et al. [2018] Robin Manhaeve, Sebastijan Dumančić, Angelika Kimmig, Thomas Demeester, and Luc De Raedt.  Deepproblog: Neural probabilistic logic programming.  In _Advances in Neural Information Processing Systems_ , volume 31, 2018.
  * Montani and Striani [2019] Stefania Montani and Manuel Striani.  Artificial intelligence in clinical decision support: a focused literature survey.  _Yearbook of Medical Informatics_ , 28(1):120–127, 2019.
  * Parafita and Vitrià [2020] Álvaro Parafita and Jordi Vitrià.  Causal inference with deep causal graphs, 2020.  URL https://arxiv.org/abs/2006.08380.
  * Paszke et al. [2019] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang, Zachary DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala.  Pytorch: An imperative style, high-performance deep learning library.  In _Advances in Neural Information Processing Systems_ , volume 32, pages 8024–8035. Curran Associates, Inc., 2019.
  * Pawlowski et al. [2020] Nick Pawlowski, Daniel Coelho de Castro, and Ben Glocker.  Deep structural causal models for tractable counterfactual inference.  In _Advances in Neural Information Processing Systems_ , volume 33, pages 857–869. Curran Associates, Inc., 2020.
  * Pearl [1988] Judea Pearl.  _Probabilistic Reasoning in Intelligent Systems: Networks of Plausible Inference_.  Morgan Kaufmann Publishers Inc., San Francisco, CA, USA, 1988.
  * Peiffer-Smadja et al. [2020] N. Peiffer-Smadja, T.M. Rawson, R. Ahmad, A. Buchard, P. Georgiou, F.-X. Lescure, G. Birgand, and A.H. Holmes.  Machine learning for clinical decision support in infectious diseases: a narrative review of current applications.  _Clinical Microbiology and Infection_ , 26(5):584–595, 2020.
  * Rohekar et al. [2018] Raanan Y. Rohekar, Shami Nisimov, Yaniv Gurwicz, Guy Koren, and Gal Novik.  Constructing deep neural networks by bayesian network structure learning.  In _Proceedings of the 32nd International Conference on Neural Information Processing Systems_ , NIPS’18, page 3051–3062, Red Hook, NY, USA, 2018. Curran Associates Inc.
  * Russell and Norvig [2009] Stuart J. Russell and Peter Norvig.  _Artificial Intelligence: a modern approach_.  Pearson, 3 edition, 2009.
  * Tayefi et al. [2021] Maryam Tayefi, Phuong Ngo, Taridzo Chomutare, Hercules Dalianis, Elisa Salvi, Andrius Budrionis, and Fred Godtliebsen.  Challenges and opportunities beyond structured data in analysis of electronic health records.  _Wiley Interdisciplinary Reviews: Computational Statistics_ , 13, 02 2021.
  * Verma and Pearl [1990] Thomas Verma and Judea Pearl.  Causal networks: Semantics and expressiveness.  In _Uncertainty in Artificial Intelligence_ , volume 9 of _Machine Intelligence and Pattern Recognition_ , pages 69–76. North-Holland, 1990.
  * Wachter et al. [2017] Sandra Wachter, Brent Mittelstadt, and Chris Russell.  Counterfactual explanations without opening the black box: Automated decisions and the gdpr.  _Harv. JL & Tech._, 31:841, 2017.

##  Appendix A Appendix

###  A.1 Sample-based training of neural network

Table 1: Absolute frequency table showing the occurrence of $X=x$ and $Y=y$ in
the hypothetical training set ($\nu=n_{00}+n_{01}+n_{10}+n_{11}$).
x | y | freq.
---|---|---
0 | 0 | $n_{00}$
0 | 1 | $n_{01}$
1 | 0 | $n_{10}$
1 | 1 | $n_{11}$

We provide a small-scale proof of how optimising the cross-entropy loss for
the observed samples leads to probabilistic outcomes. While the illustrative
setting for which we provide the proof only concerns two binary variables, the
proof can be extended to a more general setting.

Assume we have two binary variables, and observe $\nu$ samples $\\{X=x,Y=y\\}$
during every training epoch. The distribution of the training data is shown in
Table 1. Our model simultaneously optimises three partial training losses,
each corresponding to a possible evidence mask: $\mathcal{L}_{X\rightarrow Y}$
(target Y, evidence X), $\mathcal{L}_{Y\rightarrow X}$ (target X, evidence Y),
$\mathcal{L}_{XY}$ (target X and Y, no evidence). The random evidence/target
mask decides which loss is optimised during which iteration, and the sum of
these three losses forms the overall training objective. We will calculate the
optimum for each partial training loss and show that this leads to the desired
predictions for the targets at hand.

#### Evidence $X$, target $Y$:

Optimising the partial loss $\mathcal{L}_{X\rightarrow Y}$ should lead us to
an estimate for query $P(Y=y|X=x)$. The model returns a prediction
$\hat{y}_{x\rightarrow y}$ for target $Y=y$ taking evidence $X=x$ as an input.
Taking into account the frequencies of each sample as listed in Table 1, we
obtain the partial loss $\mathcal{L}_{X\rightarrow Y}$ as shown in eq. (3).
Here, we use the definition of binary cross-entropy loss as listed in eq. (2)
and simplify by filling in the possible values (0 or 1) for target y.

$\mathcal{L_{CE}}=-ylog(\hat{y})-(1-y)log(1-\hat{y})$ (2)

$\begin{split}\mathcal{L}_{X\rightarrow Y}&=-n_{01}log(\hat{y}_{0\rightarrow
1})-n_{00}log(\hat{y}_{0\rightarrow 0})-n_{11}log(\hat{y}_{1\rightarrow
1})-n_{10}log(\hat{y}_{1\rightarrow 0})\\\ &=-n_{01}log(\hat{y}_{0\rightarrow
1})-n_{00}log(1-\hat{y}_{0\rightarrow 1})-n_{11}log(\hat{y}_{1\rightarrow
1})-n_{10}log(1-\hat{y}_{1\rightarrow 1})\end{split}$ (3)

When given $X=0$ as evidence, we will optimise this loss for
$\hat{y}_{0\rightarrow 1}$. Equation (4) illustrates how calculating the
partial derivative and setting it to zero, leads to the optimum for
$\hat{y}_{0\rightarrow 1}$ (whereby $\hat{y}_{0\rightarrow 0}$ = 1 -
$\hat{y}_{0\rightarrow 1}$). A similar derivation for $X=1$ leads to the
optimum for $\hat{y}_{1\rightarrow 1}$ (and $\hat{y}_{1\rightarrow 0}$ = 1 -
$\hat{y}_{1\rightarrow 1}$). In other words, the predicted value for target Y
moves towards its relative frequency in the training set, conditioned on the
observed evidence values.

$\begin{split}\frac{\partial\mathcal{L}_{X\rightarrow
Y}}{\partial\hat{y}_{0\rightarrow 1}}=0\Rightarrow\hat{y}_{0\rightarrow
1}=\frac{n_{01}}{n_{00}+n_{01}};\quad\frac{\partial\mathcal{L}_{X\rightarrow
Y}}{\partial\hat{y}_{1\rightarrow 1}}=0\Rightarrow\hat{y}_{1\rightarrow
1}=\frac{n_{11}}{n_{10}+n_{11}}\end{split}$ (4)

#### Evidence $Y$, target $X$:

Optimising the partial loss $\mathcal{L}_{Y\rightarrow X}$ should lead us to
an estimate for query $P(X=x|Y=y)$. The derivation is symmetrical to the
previous case, with the roles of $X$ and $Y$ switched.

#### No evidence, targets $X$ and $Y$:

By optimising $\mathcal{L}_{XY}$, the model jointly optimises its predicted
outcome for queries $P(X=x)$ and $P(Y=y)$. A training case only contributes to
this partial loss when the evidence mask chooses both variables as targets.
Now, we can write the loss as in eq. (5), where $\hat{x}$ is the prediction
for $X=1$ given no evidence, and analogous for $\hat{y}$. We use the same
notation for the relative frequencies as before, and again use the definition
of binary CE-loss from eq. (2).

$\begin{split}\mathcal{L}_{XY}=-(n_{10}+n_{11})log(\hat{x})-(n_{00}+n_{01})log(1-\hat{x})\\\
-(n_{01}+n_{11})log(\hat{y})-(n_{00}+n_{10})log(1-\hat{y})\end{split}$ (5)

Again, we can optimise the loss above for $\hat{x}$ and $\hat{y}$. This leads
to the optima shown in eq. (6). The optimum for $\hat{x}$ simply corresponds
to the relative frequency of seeing $X=1$ in the training set, which is indeed
what we want as a prediction for $P(X=1)$. The same goes for $\hat{y}$.

$\begin{split}\frac{\partial\mathcal{L}_{XY}}{\partial\hat{x}}=0\Rightarrow\hat{x}=\frac{n_{10}+n_{11}}{n_{00}+n_{01}+n_{10}+n_{11}};\frac{\partial\mathcal{L}_{XY}}{\partial\hat{y}}=0\Rightarrow\hat{y}=\frac{n_{01}+n_{11}}{n_{00}+n_{01}+n_{10}+n_{11}}\\\
\end{split}$ (6)

We emphasise that the results above only hold if training happens uniformly
over the available training instances, so that the frequency of occurrence of
a target given a particular evidence actually corresponds with the empirical
probability in the training data. We ensure that this is the case by properly
shuffling our batches within each epoch.

###  A.2 Injecting independence relations through evidence corruption

We again consider the simple setting of two binary variables $X$ and $Y$, to
show that the NN+COR method works as expected. We now add the knowledge that
$X$ and $Y$ are independent. The total loss is made up of
$\mathcal{L}_{X\rightarrow Y}$, $\mathcal{L}_{Y\rightarrow X}$ and
$\mathcal{L}_{XY}$, as defined in Section A.1. Since the DAG-based corruption
is only executed when the evidence set is not empty, only the partial losses
$\mathcal{L}_{X\rightarrow Y}$ and $\mathcal{L}_{Y\rightarrow X}$ are
affected. We zoom in on how to adapt the first one according to this new
setting and how this affects the predicted outputs. The derivation for the
other partial loss is symmetrical.

Say we receive a sample $\\{X=x,Y=y\\}$ during training and the mask indicates
that $X$ is evidence, while $Y$ is the target. Since $X\perp Y$, the value of
the target $y$ should be independent of the observed evidence value $x$.
Therefore, we corrupt the value of $x$, setting it to 1 with probability
$\gamma$ and to 0 with probability $1-\gamma$. We can use the frequencies from
Table 1 to write out the contribution of each training sample to the partial
loss $\mathcal{L}_{X\rightarrow Y}$, taking into account that the evidence is
corrupted for a fraction of the training samples. This is depicted in eq. (7).
The predicted targets $\hat{y}_{0\rightarrow 1}$ and $\hat{y}_{1\rightarrow
1}$ are defined in the same way as described in Section A.1.

$\begin{split}\mathcal{L}_{X\rightarrow
Y}=&-(1-\gamma)(n_{01}+n_{11})log(\hat{y}_{0\rightarrow
1})-(1-\gamma)(n_{00}+n_{10})log(1-\hat{y}_{0\rightarrow 1})\\\
&-\gamma(n_{11}+n_{01})log(\hat{y}_{1\rightarrow
1})-\gamma(n_{10}+n_{00})log(1-\hat{y}_{1\rightarrow 1})\end{split}$ (7)

When taking the partial derivative of the loss as shown in eq. (8), we get the
optima for $Y=1$ with either value of $X$ as evidence. Due to the corruptions
we implemented in the training process, we now get $\hat{y}_{0\rightarrow 1}$
= $\hat{y}_{1\rightarrow 1}$. The prediction for $Y$ is indeed independent of
the value of $X$ (in accordance to $X\perp Y$) and simply equal to the
relative frequency of observing $Y=1$ in the training set.

$\begin{split}\frac{\partial\mathcal{L}_{X\rightarrow
Y}}{\partial\hat{y}_{0\rightarrow 1}}=0\Rightarrow\hat{y}_{0\rightarrow
1}=\frac{n_{01}+n_{11}}{n_{00}+n_{01}+n_{10}+n_{11}}\\\
\frac{\partial\mathcal{L}_{X\rightarrow Y}}{\partial\hat{y}_{1\rightarrow
1}}=0\Rightarrow\hat{y}_{1\rightarrow
1}=\frac{n_{01}+n_{11}}{n_{00}+n_{01}+n_{10}+n_{11}}\end{split}$ (8)

Note that the parameter $\gamma$ plays no role in the optimum for $\hat{y}$.
It does not matter according to which distribution we corrupt the evidence. In
our implementation, we sample uniformly over all possible classes for the
variable in question to corrupt its value. We could also opt to pull a random
sample from the training set and switch out the value of the evidence variable
to its value in this sample.

###  A.3 Bayesian network implementation

Figure 5: Conditional probability tables defining the ground-truth “Asia”
Bayesian network, based on [Lauritzen and Spiegelhalter, 1988].

We use the pgmpy Python library [Ankur Ankan and Abinash Panda, 2015] (version
0.1.17) for sampling, training and inference in our Bayesian networks.

To train our Bayesian networks from observed samples, we use the Maximum
Likelihood Estimator. This estimator studies the co-occurrence of particular
values of each variable and its parents in the training set, filling up the
CPTs as such. We use a K2 prior as a smoothing strategy, to counteract the
extremely skewed probability distributions that might appear in the CPTs when
particular combinations of variables are never observed in the training set.

Artificial samples are generated from the ground-truth Bayesian network using
the method of forward sampling with a particular seed. As a ground-truth
Bayesian model, from which our artificial training and test sets are created,
we use the “Asia” model (see Fig. 2 for its DAG structure). The CPTs that
define its joint distribution with 20 parameters are shown in Fig. 5.

We use the technique of variable elimination to perform exact inference on our
Bayesian networks. There are some cases where this method fails because the
query contains some evidence combination it has never seen before. When coming
across such a case during test time, we simply throw out the evidence for this
particular query and take $P(X)$ as an estimate. We believe this makes for a
fairer comparison than simply ignoring those queries, since our NN is in fact
able to provide an estimate for all queries, even if it has never seen a
particular evidence set before.

###  A.4 Neural network training details

Our neural network architecture and its training procedure are implemented in
Pytorch [Paszke et al., 2019].

As shown in Fig. 1, our neural network is made up of 3 layers. The network
receives an input vector of dimension $k$. In our “Asia” example, $k$ is equal
to 14 (all 7 variables have 2 classes). First, the values of the target
variables (as selected by the evidence mask) are substituted by their
respective value from $v_{0}$. The full input is then transformed to dimension
$h$, using an input-to-hidden linear layer with tanh activation. Then, a
hidden-to-hidden linear layer of size $h$, again with tanh activation is
applied. Finally, the hidden-to-output layer transforms the activations back
to dimension $k$. This layer applies $N$ softmax functions to transform the
activations belonging to each variable separately into normalised probability
values, as depicted in Fig. 1.

For initialisation of the input vector on the positions of the target
variables, we use the vector $v_{0}$ of size $k$, obtained by applying the
hidden-to-output layer to a trainable vector of size $h$, followed by the per-
variable softmax normalisation. When no evidence is present, $v_{0}$ serves as
the full input vector, and leads to the model predicting the empirical mean
probability for all variables.

We chose a hidden size $h$ of 50, since we noticed this allowed enough
flexibility while still constraining the parameter space sufficiently to avoid
overfitting. With these dimensions, our network has 4014 trainable parameters
in total. We use the Adam optimiser [Kingma and Ba, 2015] with a learning rate
of 0.001 and otherwise default parameters in Pytorch. The training samples are
fed to the neural model in batches of size 16, and the model is trained for
500 epochs. In the NN+REG model, we additionally use a regularisation batch
size of 16.

The hyperparameter $\alpha$ controls the trade-off between training and
regularisation loss. The MSE naturally has a smaller order of magnitude than
the cross-entropy loss, therefore rather large choices for $\alpha$ (10, 100,
1000) work best. We settled for $\alpha$ equal to 10 since this seemed to lead
to the best performance of the NN+REG method, though we did not observe a big
difference between any choice of $\alpha$ within the range of 1 to 1000.