\section{#1}}
\renewcommand{\theequation}{\thesection.\arabic{equation}}

\def\langle{\langle}
\def\rangle{\rangle}
\def^{\dag}{^{\dag}}
\def{\cal K}{{\cal K}}
\def{\cal N}{{\cal N}}
\font\maj=cmcsc10 \font\itdix=cmti10
\def\mbox{I\hspace{-.15em}N}{\mbox{I\hspace{-.15em}N}}
\def\mbox{I\hspace{-.15em}H\hspace{-.15em}I}{\mbox{I\hspace{-.15em}H\hspace{-.15em}I}}
\def\mbox{I\hspace{-.15em}1}{\mbox{I\hspace{-.15em}1}}
\def\rm{I\hspace{-.1em}N}{\rm{I\hspace{-.1em}N}}
\def\mbox{Z\hspace{-.3em}Z}{\mbox{Z\hspace{-.3em}Z}}
\def{\rm Z\hspace{-.2em}Z}{{\rm Z\hspace{-.2em}Z}}
\def{\rm I\hspace{-.15em}R}{{\rm I\hspace{-.15em}R}}
\def{\rm I\hspace{-.10em}R}{{\rm I\hspace{-.10em}R}}
\def\hspace{3pt}{\rm l\hspace{-.47em}C}{\hspace{3pt}{\rm l\hspace{-.47em}C}}
\def\mbox{l\hspace{-.47em}Q}{\mbox{l\hspace{-.47em}Q}}
\def\begin{equation}{\begin{equation}}
\def\end{equation}{\end{equation}}
\def\begin{enumerate}{\begin{enumerate}}
\def\end{enumerate}{\end{enumerate}}
\renewcommand{\theequation}{\thesection.\arabic{equation}}
\newcommand{{\bf 1}}{{\bf 1}}
\newtheorem{lemma}{Lemma}
\makeatletter
\usepackage{latexsym}\usepackage{bm}\def1.0{1.0}

\makeatother

\begin{document}
\title{Conjugate spinor field equation for massless spin-$\frac{3}{2}$ field in de Sitter ambient space }

\author{S. Falahi
}
\email{sfalahi27@gmail.com} \affiliation{Department of
Physics,Kermanshah Branch, Islamic Azad University,\\
Kermanshah, Iran}
\author{S. Parsamehr}
\email{sajadparsamehr@iauksh.ac.ir} \affiliation{Department of
Physics,Kermanshah Branch, Islamic Azad University,\\
Kermanshah, Iran}

\date{\today}

\begin{abstract}

\noindent \hspace{0.35cm}
The quantum field theory is expanding in de Sitter space-time. The ambient space formalism let us to develop the quantum field theory in a careful mathematical framework. Using the gauge-covariant derivative in the de Sitter ambient space, the gauge invariant Lagrangian density has been found.
In this paper, the equations of the conjugate spinor for massless spin-$\frac{3}{2}$ gravitational field is obtained by Euler-Lagrange equation. Finally, these equations can be written in terms of the operator of the Casimir group de Sitter.
\end{abstract}

\maketitle
\section{Introduction}
The experimental data and observations show that the universe is expanding with at constant positive acceleration, i.e. the space-time is curved\cite{r1,r2,r3,r4,r5,r6}.
Since the simplest curved space-time that correspondence to these observations is de Sitter space-time and this space-time has maximal symmetry, the group of ten parametric SO (1,4) is the kinematic group of the de Sitter space. Therefore, quantum field theory and gauge theory are investigated in this space-time\cite{qft,qft1,qft2,qft3,qft4,qft5,qft6,qft7,qft8}.\\
 Supersymmetry has been introduced as one of the fundamental principles of all the efforts to achieve  the grand unified theory. It acts in such a way that the relationship between the bosons (integer-valued spin) and the fermions (half-integer spin) is established. In the supersymmetry, all the particles have a partner. For example, for a graviton ( gravity-carrier particle with a spin-2), a partner with a spin-$\frac{3}{2}$ called gravitino can be considered.
With the development of quantum field theory, gauge theories and their very successful results, everyone began to think about quantization and gauge everything. As far as gauge theories are concerned, it should be considered that the structure of a particular symmetry called gauge symmetry, should remain invariant.
Quantum gravity theories that use supersymmetry are called supergravity and seek to unify the gravitational interaction with other fundamental interactions. It should be noted that supergravity is a local supersymmetry theory\cite{s1}. Therefore, gauge theory can be extended to gravity\cite{s2}.\\
To understand physical systems, we must obtain the equation of motion of physical quantities or equally system's Lagrangian.
In Lagrangian mechanics, changes in a physical system are described through solving the Euler-Lagrange equation for that system's behavior. The spin-$\frac{3}{2}$ field equation was first introduced by William Rarita and Julian Schwinger in 1941 using the Euler Lagrange method\cite{rs}.\\
In the previous work\cite{qft6}, we introduced Lagrangian for massless spin-$\frac{3}{2}$ field in de Sitter space-time. Due to the complex and difficult calculations, we intend to present more details of the calculations in this article . The notation of ambient space is briefly reviewed in Section II. In Section III, through this notation, we calculate the conjugate spinor field equation also the equation of this field is invariant. A brief conclusion is presented in Section IV. Finally, details of mathematical calculations are given in two appendices.

\section{notations}

It has been discovered today that the universe is accelerating, with a small, but non-zero and positive cosmological constant, Therefore, it can be concluded that the shape and geometry of the universe is curved. So, in the first-order approximation, we can use  the de Sitter space-time for its. This space-time is a 4-dimensional hyperboloid that can be embedded in a Minkowski 5-dimensional space-time\cite{qft7,qft8}:
\begin{equation} \label{dSs} X_H=\{x \in {\rm I\hspace{-.15em}R}^5| \; \; x \cdot x=\eta_{\alpha\beta} x^\alpha
 x^\beta =-H^{-2}\},\;\; \alpha,\beta=0,1,2,3,4, \end{equation}
where  $\eta_{\alpha\beta}=$diag$(1,-1,-1,-1,-1)$ and H is the Hubble
parameter.
The metric is defined as follows:\begin{equation} \label
{dsmet}  ds^2=\eta_{\alpha\beta}dx^{\alpha}dx^{\beta}|_{x^2=-H^{-2}}=
g_{\mu\nu}^{dS}dX^{\mu}dX^{\nu},\;\; \mu=0,1,2,3,\end{equation}
$X^\mu$ is the four components of the space-time coordinates in a system of intrinsic coordinates on a hyperboloid and $x^\alpha$ is the five dimensional Minkowski space-time(de Sitter ambient space).
Two Casimir operators of the group include:
 \begin{equation} \label{casi1} Q^{(1)}=-\frac{1}{2}L_{\alpha\beta}L^{\alpha\beta},\;\;  \alpha, \beta=0,1,2,3,4, \end{equation}
      \begin{equation} Q^{(2)}=-W_\alpha W^\alpha\;\;,\;\;W_\alpha =\frac{1}{8}
      \epsilon_{\alpha\beta\gamma\delta\eta} L^{\beta\gamma}L^{\delta\eta},\end{equation}
			 $\epsilon_{\alpha\beta\gamma\delta\eta}$ is an anti-symmetric tensor, and $L_{\alpha\beta}$ are ten infinitismal generators in de Sitter space. They can be written as a linear combination:
			$L_{\alpha\beta}=M_{\alpha\beta}+S_{\alpha\beta}$.
			Where $M_{\alpha\beta}$ is the orbital part and $S_{\alpha\beta}$ is the spinoral part. In this formalism, the space $M_{\alpha\beta}$ is represented as follow:
			\begin{equation} \label{genm} M_{\alpha \beta}=-i(x_\alpha \partial_\beta-x_\beta
      \partial_\alpha)=-i(x_\alpha\partial^\top_\beta-x_\beta
        \partial^\top_\alpha),\end{equation}
				where $\partial^\top_\beta=\theta_\beta^{\;\;\alpha}\partial_\alpha$ is the transverse derivative ($x.\partial^\top=0$) and $\theta_{\alpha\beta}=\eta_{\alpha\beta}+H^2x_\alpha x_\beta$ is the projection tensor on de sitter hyperboloid. For half-integer spin fields  $s = l + 1/2$, the spinoral part is defined as:
\begin{equation}
S_{\alpha\beta}^{(s)}=S_{\alpha\beta}^{(l)}+S_{\alpha\beta}^{(\frac{1}{2})},
\end{equation}
where $S_{\alpha\beta}$for spin $\frac{1}{2}$field is:
\begin{equation}
S_{\alpha\beta}=-\frac{i}{4}[\gamma_{\alpha},\gamma_{\beta}],\end{equation}
the $\gamma$-matrices have to satisfy following relation:
\begin{equation}\{
{\gamma^{\alpha},\gamma^{\beta}}
\}=2\eta^{\alpha\beta}\mathbb{I}.\end{equation}
 A proper display for them is \cite{qft7,qft8}:
$$ \gamma^0=\left( \begin{array}{clcr} \mathbb{I} & \;\;0 \\ 0 &-\mathbb{I} \\ \end{array} \right)
      ,\gamma^4=\left( \begin{array}{clcr} 0 & \mathbb{I} \\ -\mathbb{I} &0 \\ \end{array} \right) ,  \label{gammam}
   \gamma^1=\left( \begin{array}{clcr} 0 & i\sigma^1 \\ i\sigma^1 &0 \\
    \end{array} \right)
   ,\gamma^2=\left( \begin{array}{clcr} 0 & -i\sigma^2 \\ -i\sigma^2 &0 \\
      \end{array} \right)
   , \gamma^3=\left( \begin{array}{clcr} 0 & i\sigma^3 \\ i\sigma^3 &0 \\
      \end{array} \right),$$ \begin{equation}\gamma^{\alpha\dagger}=\gamma^{0}\gamma^{\alpha}\gamma^{0} \qquad  (\gamma^4)^2=-1 \qquad (\gamma^0)^2=1 \end{equation}
			where $\mathbb{I}$ is  unit $2\times 2$ matrix and $\sigma^i$ are the pauli matrices.
It should be noted that for massless spin-$\frac{3}{2}$ field $Q^{(1)}_\frac{3}{2}$ is\cite{qft6}:
\begin{equation}
Q_{\frac{3}{2}}^{(1)}\Psi_\alpha=Q_0^{(1)}\Psi_\alpha+\not{x}\not{\partial}^\top\Psi_\alpha
+2x_\alpha\partial^\top\cdot\Psi-\frac{11}{2}\Psi_\alpha+\gamma_\alpha
\not{\Psi},
\end{equation}
$Q_0^{(1)}=-\partial_\alpha^\top\partial^{\alpha\top}$ is the "scalar" Casimir operator.

\section{conjugate spinor field equation for massless spin-$\frac{3}{2}$ field }

It is believed that the gauge theory is the basis of fundamental particle interactions. The Lagrangian for  massless spin-$\frac{3}{2}$ field is peresented, in the linear approximation, by using the gauge theory and defining the gauge covariant derivative\cite{qft6}. The vector-spinor field equation $\Psi_\alpha(x)$ is obtained from the usual Euler-Lagrange equations. This Lagrangian is invariant under the  gauge transformation (see more details [\cite{epjp}]):
\begin{equation}\label{A.1}
{\cal L}=\left(\tilde{\nabla}^\top_\alpha\tilde{
\Psi}_\beta-\tilde{\nabla}^\top_\beta\tilde{
\Psi}_\alpha\right)\left( \nabla^{\top\alpha}
\Psi^\beta-\nabla^{\top\beta} \Psi^\alpha\right),
\end{equation}
$\nabla^\top_\alpha$ is a transverse-covariant derivative which is defined to obtain an invariant Lagrangian according to the following equation:
\begin{equation} \label{cdsa} \nabla^\top_\beta \Psi_{\alpha_1....\alpha_l}\equiv
(\partial^\top_\beta+\gamma^\top_\beta\not
x)
\Psi_{\alpha_1....\alpha_l}-\sum_{n=1}^{l}x_{\alpha_n}\Psi_{\alpha_1..\alpha_{n-1}\beta\alpha_{n+1}..\alpha_l}.\end{equation}
also for the conjugate spinor$\tilde{\Psi}_\alpha$:
\begin{equation} \label{cdsa} \tilde{\nabla}^\top_\beta \tilde{\Psi}_{\alpha_1....\alpha_l}\equiv
\partial^\top_\beta
\tilde{\Psi}_{\alpha_1....\alpha_l}-\sum_{n=1}^{l}x_{\alpha_n}\tilde{\Psi}_{\alpha_1..\alpha_{n-1}\beta\alpha_{n+1}..\alpha_l}.\end{equation}
where $\not x=\gamma_\alpha x^\alpha$and $\gamma^\top_\alpha=\theta_{\alpha}^{\beta}\gamma_\beta$.
The above equations are specifically designed for our calculations:
\begin{equation}
\nabla^\top_\alpha \Psi_\beta= \partial^\top_\alpha \Psi_\beta+ \gamma^\top_\alpha \not x  \Psi_\beta -x_\beta \Psi_\alpha.
\end{equation}
\begin{equation}\tilde{\nabla}^\top_\beta \tilde{\Psi}_\alpha=\partial^\top_\beta
\tilde{\Psi}_\alpha-x_{\alpha}\tilde{\Psi}_\beta
\end{equation}

Now we want to obtain the field equation for the conjugate spinor($\tilde{\Psi}_\alpha={\Psi}^\dag_\alpha\gamma^0$) by using the Euler-Lagrange equation in the linear approximation. The Euler-Lagrange equation is:
\begin{equation}
\frac{\delta {\cal L}}{\delta\tilde{
\Psi}_m}-\partial_{l}^{\top}\frac{\delta {\cal
L}}{\delta(\partial^{\top}_{ l}\tilde{ \Psi}_m)}=0,
\end{equation}
First, we extend each terms of the equation (III.11):
\begin{equation}
A=\left( \nabla^{\top\alpha} \Psi^\beta-\nabla^{\top\beta}
\Psi^\alpha\right)=\left(\partial^{\top\alpha}
\Psi^\beta+\gamma^\alpha\not {x}\Psi^\beta -\partial^{\top\beta}
\Psi^\alpha-\gamma^\beta\not {x}\Psi^\alpha \right),
\end{equation}
\begin{equation}
B=\left(\tilde{\nabla}^\top_\alpha\tilde{
\Psi}_\beta-\tilde{\nabla}^\top_\beta\tilde{ \Psi}_\alpha\right)
=\left(\partial_{\alpha}^{\top}\tilde{ \Psi}_\beta-x_\beta\tilde{
\Psi}_\alpha-\partial_{\beta}^{\top}\tilde{
\Psi}_\alpha+x_\alpha\tilde{ \Psi}_\beta\right),
\end{equation}
By use the Euler-Lagrange equation, we consider the following terms:
\begin{equation}
\frac{\delta {\cal L}}{\delta \Psi_
{m}}=(\gamma^\alpha\not {x} \delta_{\beta}^{m}-\gamma^\beta\not {x} \delta_{\alpha}^{m})\left(\tilde{\nabla}^\top_\alpha\tilde{
\Psi}_\beta-\tilde{\nabla}^\top_\beta\tilde{ \Psi}_\alpha\right),
\end{equation}
and
\begin{equation}
\frac{\delta {\cal L}}{\delta(\partial^{\top}_{ l} \Psi_{m})}=
(\delta_{l}^{\alpha}\delta_{m}^{\beta}-\delta_{l}^{\beta}\delta_{m}^{\alpha})\left(\tilde{\nabla}^\top_\alpha\tilde{
\Psi}_\beta-\tilde{\nabla}^\top_\beta\tilde{ \Psi}_\alpha\right),
\end{equation}
if $\beta=m$, then, one obtains:
\begin{equation}
\frac{\delta {\cal L}}{\delta \Psi_
{m}}=\gamma^\alpha\not {x} \left(\tilde{\nabla}^\top_\alpha\tilde{
\Psi}_\beta-\tilde{\nabla}^\top_\beta\tilde{ \Psi}_\alpha\right),
\end{equation}
\begin{equation}
\frac{\delta {\cal L}}{\delta(\partial^{\top}_{ l} \Psi_{m})}=
\delta_{l}^{\alpha}\left(\tilde{\nabla}^\top_\alpha\tilde{
\Psi}_\beta-\tilde{\nabla}^\top_\beta\tilde{ \Psi}_\alpha\right).
\end{equation}
Now, by placing the above expressions in the Euler-Lagrange equation, the field equation is obtained as follows:
\begin{equation}
\gamma^\alpha\not {x} \left(\tilde{\nabla}^\top_\alpha\tilde{
\Psi}_\beta-\tilde{\nabla}^\top_\beta\tilde{ \Psi}_\alpha\right)-\partial^{\top\alpha} \left(\tilde{\nabla}^\top_\alpha\tilde{
\Psi}_\beta-\tilde{\nabla}^\top_\beta\tilde{ \Psi}_\alpha\right)=0,
\end{equation}
this equation can be written in the summarized form:
\begin{equation}
\left(\partial^{\top\alpha}-\gamma^\alpha\not x \right)
\left(\tilde{\nabla}^\top_\alpha\tilde{
\Psi}_\beta-\tilde{\nabla}^\top_\beta\tilde{ \Psi}_\alpha\right)=0.
\end{equation}

In the appendix A, the field equations for the conjugate spinor is obtained by using the second order Casimir operator:
\begin{equation}
\left(Q_{\frac{3}{2}}^{(1)}+\frac{5}{2}\right)\tilde
\Psi_\alpha+\partial_{\alpha}^\top\left(\not x\tilde{\not \Psi}
+\partial^\top\cdot\tilde\Psi\right)-2\left(\gamma_\alpha\tilde{\not \Psi}
+\not x\not \partial^\top\tilde\Psi_\alpha-\tilde\Psi_\alpha\right)=0.
\end{equation}

In this here we present the gauge invariant field equation for the conjugate spinor.
Given the definition of $Q_{\frac{3}{2}}^{(1)}$, we rewrite Equation (III.25) as follows:
\begin{equation}
-Q_{0}\tilde{ \Psi}_\beta+\tilde{ \Psi}_\beta-2x_\beta(\partial^{\top}\cdot\tilde{
\Psi})-\partial_{\beta}^{\top}(\partial^{\top}\cdot\tilde{
\Psi})+\not x\not \partial^{\top} \tilde{
\Psi}_\beta-x_\beta\not x \not \tilde{\Psi}-\not x\partial_{\beta}^{\top}\not\tilde{
\Psi}=0
\end{equation}
We show that the above equation is invariant under the following gauge transformation:
\begin{equation}
\tilde{\Psi}_\alpha \longrightarrow
\tilde{\Psi}^g_\alpha=\tilde{\Psi}_\alpha+\partial_{\alpha}^\top
\tilde{\psi}.
\end{equation}
Therefore, equation (III.26) comes as follows:
$$
-Q_{0}(\tilde{\Psi}_\beta+\partial_{\beta}^\top
\tilde{\psi})+\tilde{\Psi}_\beta+\partial_{\beta}^\top
\tilde{\psi}-2x_\beta\partial^{\alpha\top}(\tilde{\Psi}_\alpha+\partial_{\alpha}^\top
\tilde{\psi})$$
\begin{equation}
-\partial_{\beta}^{\top}\partial^{\alpha\top}(\tilde{\Psi}_\alpha+\partial_{\alpha}^\top
\tilde{\psi})+\not x\not \partial^{\top} (\tilde{\Psi}_\beta+\partial_{\beta}^\top
\tilde{\psi})-x_\beta\not x (\tilde{\slashed{\Psi}}+\slashed{\partial}^\top\tilde{\psi})
-\not x\partial_{\beta}^{\top}(\tilde{\slashed{\Psi}}+\slashed{\partial}^\top\tilde{\psi})=0
\end{equation}
After we calculate a little:

$$
\Longrightarrow=-Q_{0}\tilde{\Psi}_\beta+\tilde{\Psi}_\beta-2x_\beta(\partial^{\top}\cdot\tilde{
\Psi})-\partial_{\beta}^{\top}(\partial^{\top}\cdot\tilde{
\Psi})+\not x\not \partial^{\top} \tilde{
\Psi}_\beta-x_\beta\not x \not \tilde{\Psi}-\not x\partial_{\beta}^{\top}\not\tilde{
\Psi}$$
\begin{equation}
\underbrace{-Q_{0}\partial_{\beta}^\top\tilde{\psi}+\partial_{\beta}^\top
\tilde{\psi}-2x_\beta\partial^{\alpha\top}\partial_{\alpha}^{\top}\psi-\partial_{\beta}^{\top}\partial^{\alpha\top}\partial_{\alpha}^{\top}\psi
+\slashed{x}\slashed{\partial}^\top\partial_{\beta}^{\top}\tilde{\psi}-x_\beta\not x \not\slashed{\partial}^\top\tilde{\psi}-\slashed{x}
\partial_{\beta}^{\top}\slashed{\partial}^\top\tilde{\psi}}_{=0}=0
\end{equation}

For a gauge invariant, the latest terms must be zero, We have to prove this. To do this we use the following auxiliary relationships:
\begin{equation}
[\partial^\top_\alpha , \not\partial^\top]=\not x\partial^\top_\alpha-x_\alpha\not\partial^\top
\end{equation}
\begin{equation}
[\partial^\top_\alpha , Q_0]=-6\partial^\top_\alpha-2(Q_0+4)x_\alpha
\end{equation}
Thus, the latest terms of equation (III.29) are obtained as follows:
$$
-Q_{0}\partial_{\beta}^\top\tilde{\psi}+\partial_{\beta}^\top
\tilde{\psi}-2x_\beta\partial^{\alpha\top}\partial_{\alpha}^{\top}\psi-\partial_{\beta}^{\top}\partial^{\alpha\top}\partial_{\alpha}^{\top}\psi
+\slashed{x}\slashed{\partial}^\top\partial_{\beta}^{\top}\tilde{\psi}-x_\beta\not x \not\slashed{\partial}^\top\tilde{\psi}-\slashed{x}
\partial_{\beta}^{\top}\slashed{\partial}^\top\tilde{\psi}
$$
\begin{equation}
=2x_\beta Q_{0}\tilde{\psi}+\partial_{\beta}^\top
\tilde{\psi}-x_\beta\not x \not\slashed{\partial}^\top\tilde{\psi}-6
\partial_{\beta}^{\top}\tilde{\psi}-2Q_{0}x_\beta\tilde{\psi}-8x_\beta \tilde{\psi}+x_\beta\not x \not\slashed{\partial}^\top\tilde{\psi}
+\partial_{\beta}^{\top}\tilde{\psi}
\end{equation}

After a few of simplification:
\begin{equation}
=2x_\beta Q_{0}\tilde{\psi}-4\partial_{\beta}^\top
\tilde{\psi}-2Q_{0}x_\beta\tilde{\psi}-8x_\beta \tilde{\psi}=2(x_\beta Q_{0}-Q_{0} x_\beta )\tilde{\psi}-4\partial_{\beta}^\top
\tilde{\psi}-8x_\beta \tilde{\psi}
\end{equation}
And finally using an auxiliary relationship $[x_\alpha , Q_0]=2\partial^\top_\alpha+4x_\alpha$
 we have:
\begin{equation}
=2(2\partial^\top_\beta+4x_\beta)\tilde{\psi}-4\partial_{\beta}^\top
\tilde{\psi}-8x_\beta \tilde{\psi}=0
\end{equation}

Therefore, we prove that the equation is the invariant field.

\section{Conclusions}

In order to better understand the evolution of the universe, it is necessary to extend the theory of quantum fields, field interactions, or gauge theory, supersymmetry and supergravity in  the de Sitter space-time.
We studied the conjugate spinor field equation by the Euler-Lagrange equation in the de Sitter ambient space formalism. Furthermore, the field equation in terms of the Casimir operator is obtained also the gauge invariant field equation is presented. These studies could be appropriate for supergravity theory in curved space, which could be presented in the future.

\vspace{0.5cm}

\noindent {\bf{Acknowlegements}}:  I would like to express our heartfelt thank and sincere gratitude to Professor
M.V. Takook for his helpful discussions. This work has been supported by the Islamic Azad University, Kermanshah Branch, Kermanshah, Iran.


\section{Sensitivity studies}
\label{sec:appendix}

Here we provide results on a series of comprehensive systematic studies that aim to quantify the sensitivity of the resulting predictions on: (i) the neural network initialization, (ii) the total number of training and collocation points, (iii) the neural network architecture, and (iv) the adversarial training procedure. In all cases we have used the non-linear Burgers defined in section \ref{sec:burgers} as a prototype problem.

\subsection{Sensitivity with respect to the neural network initialization}

In order to quantify the sensitivity of the proposed methods with respect to the initialization of the neural networks, we have considered a noise-free data set comprising of $N_u = 150$ and $N_r = 10000$ training and collocation points, respectively, and fixed the architecture for generator neural networks to include 4 hidden layers with 50 neurons each and discriminator neural networks to include 3 hidden layers with 50 neurons each
, and a hyperbolic tangent activation function. Then we have trained an ensemble of 15 cases all starting from a normal Xavier initialization \citep{glorot2010understanding} for all network weights (with a randomized seed), and a zero initialization for all bias parameters. In table \ref{tab:sens_t1} we report the relative error between the predicted mean solution and the known exact solution for this problem for all 15 randomized trials using at set of $25600$ randomly selected test points. Evidently, our results are robust with respect to the the neural network initialization as in all cases the stochastic gradient descent training procedure converged roughly to the same solution. We can summarize this result by reporting the mean and the standard deviation of the relative $\mathbb{L}_2$ error as
\[
\hat{\mathcal{L}_2} \in [\mu_L-\sigma_L, \mu_L+\sigma_L] = [4.7\times 10^{-2} - 1.3\times 10^{-2}, 4.7\times 10^{-2} + 1.3\times 10^{-2}].
\]

\begin{table}
\centering
\begin{tabular}{|c|c|c|c|c|c|}
\hline
\multicolumn{5}{|c|}{Relative $\mathcal{L}_2$ error}\\
\hline
4.1e-02& 7.9e-02& 4.4e-02& 4.0e-02& 3.8e-02\\
\hline
3.2e-02& 5.7e-02& 4.7e-02& 6.5e-02& 4.0e-02\\
\hline
3.5e-02& 3.5e-02& 6.4e-02& 4.0e-02& 4.9e-02\\
\hline
\end{tabular}
\caption{Relative $\mathbb{L}_2$ prediction error for different neural network initializations using a randomized seed.}
\label{tab:sens_t1}
\end{table}

\subsection{Sensitivity with respect to the total number of training and collocation points}\label{B2}

In this study our goal is to quantify the sensitivity of our predictions with respect to the total number of training and collocation points $N_u$ and $N_r$, respectively. As before, we have considered noise-free data sets, and fixed the architecture for generator neural networks to include 4 hidden layers with 50 neurons each and discriminator neural networks to include 3 hidden layers with 50 neurons each, a hyperbolic tangent activation function, and a normal Xavier initialization \citep{glorot2010understanding} for all network weights and zero initialization for all network biases. The results of this study are summarized in table \ref{tab:sens_t2}, indicating that as the number of collocation points are increased, a more accurate prediction is obtained. This observation is in agreement with the original results of Raissi {\em et. al.} \citep{raissi2017physics1,raissi2017physics2} for deterministic physics-informed neural networks, indicating the role of the residual loss as an effective regularization mechanism for training deep generative models in small data regimes.

\begin{table}
\centering
\begin{tabular}{|c|c|c|c|c|c|c|c|c|c|}
\hline
\diagbox{$N_u$}{$N_r$} & 10& 100& 250& 500& 1000& 5000& 10000\\
\hline
60  & 9.3e-01& 5.6e-01& 4.8e-01& 5.0e-02& 1.9e-01& 5.0e-02& 5.1e-02\\
\hline
90  & 5.8e-01& 5.3e-01& 3.5e-01& 1.5e-01& 4.9e-02& 1.0e-01& 5.8e-02\\
\hline
150 & 6.7e-01& 1.4e-01& 3.0e-01& 3.6e-02& 4.9e-02& 1.2e-01& 4.7e-02\\
\hline
\end{tabular}
\caption{Relative $\mathbb{L}_2$ prediction error for different number of training and collocation points $N_u$ and $N_r$, respectively.}
\label{tab:sens_t2}
\end{table}

\subsection{Sensitivity with respect to the neural network architecture}

In this study we aim to quantify the sensitivity of our predictions with respect to the architecture of the neural networks that parametrize the generator, the discriminator, and the encoder. Here we have fixed the number of noise-free training data to $N_u = 150$ and $N_r =10000$, and we kept the number of layers for discriminator to always be one less than the number of layers for generator (e.g., if the number of layers for generator is two then the number of layers for discriminator is one, etc.).
In all cases, we have used a hyperbolic tangent non-linearity and a normal Xavier initialization \citep{glorot2010understanding}. In table \ref{tab:sens_t3} we report the relative $\mathbb{L}_2$ prediction error for different feed-forward architectures for the generator, discriminator, and encoder (i.e., different number of layers and number of nodes in each layer). The general trend suggests that as the neural network capacity is increased we obtain more accurate predictions, indicating that our physics-informed constraint on the PDE residual can effectively regularize the training process and safe-guard against over-fitting. We note number of neurons in each layer as $N_{n}$ and number of layers for generator (encoder) as $N_g$.

\begin{table}
\centering
\begin{tabular}{|c|c|c|c|}
\hline
\diagbox{$N_g$}{$N_n$} & 20& 50& 100\\
\hline
2 &    4.2e-01& 3.8e-01& 5.7e-01 \\
\hline
3 &    6.5e-02& 3.5e-02& 2.1e-02 \\
\hline
4 &    9.3e-02& 4.7e-02& 5.4e-02 \\
\hline
\end{tabular}
\caption{Relative $\mathbb{L}_2$ prediction error for different feed-forward architectures for the generator, encoder, and the discriminator. The total number of layers of the latter was always chosen to be one less than the number of layers for generator.}
\label{tab:sens_t3}
\end{table}

\subsection{Sensitivity with respect to the adversarial training procedure}

Finally, we test the sensitivity with respect to the adversarial training process. To this end, we have fixed the number of noise-free training data to $N_u = 150$ and $N_r =10000$, and the neural network architecture to be the same as \ref{B2}, and we vary the total number of training steps for the generator $K_g$ and the discriminator $K_d$ within each stochastic gradient descent iteration. The results of this study are presented in table \ref{tab:sens_t4} where we report the relative $\mathbb{L}_2$ prediction error. These results reveal the high sensitivity of the training dynamics on the interplay between the generator and discriminator networks, and pinpoint on the well known peculiarity of adversarial inference procedures which require a careful tuning of
$K_g$ and $K_d$ for achieving stable performance in practice.

\begin{table}
\centering
\begin{tabular}{|c|c|c|c|}
\hline
\diagbox{$K_g$}{$K_d$} & 1& 2& 5\\
\hline
1 &    3.5e-01& 5.0e-01& 1.5e+00 \\
\hline
2 &    4.3e-02& 3.2e-01& 5.4e-01 \\
\hline
5 &    4.7e-02& 2.3e-01& 7.0e-01 \\
\hline
\end{tabular}
\caption{Relative $\mathbb{L}_2$ error with different number of training for generator and discriminator in each epoch.}
\label{tab:sens_t4}
\end{table}

\section{Introduction}\label{sec:intro}
Recent advances in machine learning and data analytics have yielded transformative results across diverse scientific disciplines, including image recognition \cite{krizhevsky2012imagenet}, natural language processing \cite{lecun2015deep}, cognitive science \cite{lake2015human}, and genomics \cite{alipanahi2015predicting}.
In all aforementioned areas, the volume of data has increased substantially compared to even a decade ago, but analyzing big data is expensive and time-consuming. Data-driven methods, which have been enabled  by the availability of sensors, data storage, and computational resources, are taking center stage across many disciplines of science. We now have highly scalable solutions for problems in object detection and recognition, machine translation, text-to-speech conversion, recommender systems, and information retrieval \cite{lecun2015deep}. All of these solutions attain state-of-the-art performance when trained with large amounts of data.

However, more often than not, in laboratory experiments and large-scale simulations aiming to elucidate and predict complex phenomena, a large number of quality and error-free data is prohibitively costly to obtain. Under this setting, purely data-driven approaches for machine learning present difficulties when the data is scarce relative to the complexity of the system. The vast majority of state-of-the art machine learning techniques (e.g., deep neural nets, convolutional networks, recurrent networks, etc. \cite{goodfellow2016deep}) are lacking robustness and fail to provide any guarantees of convergence or quantify the error/uncertainty associated with their predictions. Hence, the ability to learn in a robust and sample-efficient manner is a necessity in these data-limited domains. Even less well understood is how one can constrain such algorithms to leverage domain-specific knowledge and return predictions that satisfy certain physical principles (e.g., conservation of mass, momentum, etc.).

These shortcomings often generate skepticism and disbelief among applied mathematicians and engineers regarding the solid grounding of purely data-driven machine learning approaches. In recent work, Raissi {\em et. al.} \cite{raissi2017physics1,raissi2017physics2,raissi2018multistep,raissi2018hidden,raissi2018deep} set foot exactly at this relatively unexplored interface between applied mathematics and contemporary machine learning by revisiting the idea of penalizing the loss function of deep neural networks using differential equation constraints, as first put forth by Psichogios and Ungar \cite{psichogios1992hybrid} and Lagaris {\em et. al.} \cite{lagaris1998artificial}.
This line of work has empirically demonstrated how such physics-informed constraints regularize learning in {\em small data} regimes, can lead to the discovery of governing equations and reduced-order models, as well as enable the prediction of complex dynamics from incomplete models and incomplete data. Despite a series of impressive results in canonical problems, Raissi {\em et. al.} \cite{raissi2017physics1,raissi2017physics2} have also pointed out cases in which the training phase of these algorithms faces severe difficulties for reasons that are currently poorly understood. In lack of supporting theory on convergence and a-posteriori error estimation, this naturally poses the need for   for scalable algorithms for uncertainty quantification.

A literature review of the current state-of-the-art in uncertainty quantification reveals a subtle dichotomy between different communities. On one hand, researchers in applied mathematics and scientific computing predominately rely on mathematical models that are rigorously derived from first physical principles. At the dawn of exascale computing, such models have enabled the accurate simulation of increasingly more complex phenomena (see for e.g., \cite{perdikaris2016multiscale,rossinelli2015silico}). They have also enabled {\em in-silico} systematic studies in which the behavior of a system can be probed in a controlled fashion for different conditions, parameter settings, external inputs, etc. \cite{vsukys2017optimal}. The latter aims to both elucidate the key mechanisms that govern the behavior of a system, but also characterize the robustness of the resulting predictions with respect to epistemic and aleatory uncertainty \cite{oden2010computer}. However, despite the fact that much progress has been made over the last two decades, the most popular  methods for scientific computing under uncertainty, such as polynomial chaos expansions \cite{ghanem1990polynomial, xiu2002wiener, najm2009uncertainty}, sparse grid quadratures \cite{gerstner1998numerical, eldred2009comparison}, multi-level/multi-fidelity Monte Carlo sampling \cite{barth2011multi, peherstorfer2016optimal}, proper orthogonal decomposition \cite{berkooz1993proper, le2010spectral}, and Gaussian process regression models \cite{bilionis2012multi, bilionis2013multi, perdikaris2016multifidelity}, all face severe limitations in view of the non-Gaussian likelihoods and high-dimensional posterior distributions commonly encountered in realistic applications.

On the other hand, the recent explosive growth of machine learning research has put forth new effective ways of learning and manipulating complex high-dimensional probability distributions.  Inference tools like variational auto-encoders \cite{kingma2013auto} and generative adversarial networks \cite{goodfellow2014generative}, formulated on top of flexible building blocks such as feed-forward/convolutional/recurrent neural networks \cite{goodfellow2016deep}
have introduced highly scalable solutions, albeit for problems where not much prior information is assumed, but instead, large amounts of data can be obtained at relatively low cost (e.g., image recognition \cite{krizhevsky2012imagenet}, natural language processing \cite{lecun2015deep}).  In this work, we aim to leverage recent developments in machine learning to put forth a scalable framework for uncertainty propagation in physical systems for which
the cost of data acquisition is high and training  data-sets  are  typically  small, but strong prior information exists by means of known governing laws expressed by partial differential equations. Specifically, we construct a class of probabilistic physics-informed neural networks that enables us to obtain a posterior characterization of the uncertainty associated with their predicted outputs. Moreover, we will develop a flexible variational inference framework that will allow us to train such models directly from noisy input/output data, and predict outcomes of non-linear dynamical systems that are partially observed with quantified uncertainty. This setting necessitates a departure from the classical deterministic realm of modeling and scientific computation, and, consequently, our main building blocks can no longer be crisp deterministic numbers and governing laws, but instead we must operate with {\em probabilistic models}.

This paper is structured as follows. In section \ref{sec:PINNs} we provide a brief overview of physics-informed neural networks in sync with the recent developments in \cite{raissi2017physics1,raissi2017physics2,raissi2018hidden,raissi2018deep}. In sections \ref{sec:PPINNs} and \ref{sec:ADVI} we provide an outline of the proposed probabilistic formulation and the proposed variational inference framework. Finally, in section \ref{sec:Results} we will demonstrate the effectiveness of our approach through a series of examples involving uncertainty propagation in non-linear conservation laws, and the discovery of constitutive laws for flow through porous media directly from noisy data.

\section{Methods}\label{sec:methods}

\subsection{Physics-informed neural networks}\label{sec:PINNs}
The recent works of Raissi {\em et. al.} \cite{raissi2017physics1,raissi2017physics2, raissi2018hidden,raissi2018deep} have  demonstrated how classical conservation laws and numerical discretization schemes can be used as structured prior information that can enhance the robustness and efficiency of modern machine learning algorithms, introducing a new class of data-driven solvers, as well as a {\em physics-informed machine learning} approach to model discovery. To this end, the authors have considered constructing deep neural networks that return predictions which are constrained by parametrized partial differential equations (PDE) of the form

\begin{equation}
\bm{u}_t + \mathcal{N}_{\bm{x}}\bm{u} = 0,
\end{equation}
where $\bm{u}(\bm{x},t)$ is represented by a deep neural network parametrized by a set of parameters $\theta$, i.e. $u(\bm{x},t) = f_{\theta}(\bm{x},t)$, $\bm{x}$ is a vector of space coordinates, $t$ is the time coordinate, and $\mathcal{N}_{\bm{x}}$ is a nonlinear differential operator. As neural networks are differentiable representations, this construction defines a so-called {\em physics informed neural network} that corresponds to the PDE residual, i.e. $r_{\theta}(\bm{x},t):= \frac{\partial}{\partial t} f_{\theta}(\bm{x},t)  + \mathcal{N}_{\bm{x}}f_{\theta}(\bm{x},t)$. This new network has the same parameters as the network representing $\bm{u}(\bm{x},t)$, albeit different activation functions due to the action of the differential operator \cite{raissi2017physics1, psichogios1992hybrid, lagaris1998artificial}. From an implementation perspective, this network can be readily obtained by leveraging recent progress in automatic differentiation \cite{baydin2015automatic, abadi2016tensorflow}.
The resulting training procedure allows us to recover the shared network parameters $\theta$ using a few scattered observations of $\bm{u}(\bm{x},t)$, namely $\{(\bm{x}_{i}, t_i), \bm{u}_i\}$, $i = 1,\dots,N_u$, along with a larger number of collocation points $\{(\bm{x}_{i}, t_i), \bm{r}_i = 0\}$, $i = 1,\dots,N_r$, that aim to penalize the PDE residual at a finite set of $N_r$ collocation nodes. This simple, yet remarkably effective regularization procedure allows us to introduce the PDE residual as a soft penalty constraint penalty in the likelihood function of the model \cite{raissi2017physics1,raissi2017physics2}, and the resulting optimization problem can be effectively solved using standard stochastic gradient descent without necessitating any elaborate constrained optimization techniques, simply by minimizing the composite loss function

\begin{equation}\label{eq:PINN_loss}
\mathcal{L}(\theta) = \frac{1}{N_u}\sum\limits_{i=1}^{N_u}\|f_{\theta}(\bm{x}_{i}, t_i) - \bm{u}_i\|^2  +   \frac{1}{N_r}\sum\limits_{i=1}^{N_r}\|r_{\theta}(\bm{x}_{i}, t_i) - \bm{r}_i\|^2,
\end{equation}
where the required gradients $\frac{\partial \mathcal{L}}{\partial \theta}$ can be readily obtained using automatic differentiation \cite{baydin2015automatic}.
Finally, as the resulting predictions are encouraged to inherit any physical properties imposed by the PDE constraint (e.g., conservation, invariance, symmetries, etc.), this approach showcases how one can approximately encode physical and domain-specific constraints in modern machine learning algorithms and introduce a new form of regularization for learning from {\em small} data-sets.

\subsection{Probabilistic physics-informed neural networks}\label{sec:PPINNs}

Here we put forth a probabilistic formulation for propagating uncertainty through physics-informed neural networks using latent variable models of the form

\begin{equation}
\label{eq:PSC}
p(\bm{u}|\bm{x},t,\bm{z}), \ \ \bm{z}\sim p(\bm{z}), \ \ \text{s.t} \ \  \bm{u}_t + \mathcal{N}_{\bm{x}}\bm{u} = 0
\end{equation}

This setting encapsulates a wide range of deterministic and stochastic problems, where $\bm{u}(\bm{x},t)$ is a potentially multi-variate field, and $\bm{z}$ is a collection of random latent variables. The ability to learn such a model from data is the cornerstone of probabilistic scientific computing and uncertainty quantification in physical systems. Knowledge of the conditional probability $p(\bm{u}|\bm{x},t,\bm{z})$ subject to domain knowledge constraints introduces a regularization mechanism that limits the space of admissible solutions to a manageable size (e.g., in fluid mechanics problems by discarding any non-realistic flow solutions that violate the conservation of mass principle), thus enables training of probabilistic deep learning algorithms in {\em small data} regimes.
Moreover, by providing a complete characterization of uncertainty, it enhances the robustness of our predictions, and provides a-posteriori error estimates for assessing model inadequacy. The latter, can also enable downstream tasks such as the formulation of adaptive data acquisition policies for active learning or Bayesian optimization \cite{shahriari2016taking} with domain knowledge constraints. Finally, thanks to the structure encoded by the PDE itself, the resulting latent variables $\bm{z}$  can potentially lead to the extraction of physically relevant and interpretable low-dimensional feature representations, which can  subsequently introduce new techniques for nonlinear model order reduction and coarse-graining of complex systems.

\subsection{Adversarial inference for joint distribution matching}\label{sec:ADVI}

Following the recent findings of \cite{li2018learning} we argue that matching the joint distribution of the generated data $p_{\theta}(\bm{x},t,\bm{u})$ with the joint distribution of the observed data $q(\bm{x},t,\bm{u})$ by
 minimizing the reverse Kullback-Leibler divergence $\mathbb{KL}[p_{\theta}(\bm{x},t,\bm{u})||q(\bm{x},t,\bm{u})]$ is a promising approach to train the generative model presented in equation \ref{eq:PSC}. This also implies that the respective marginal and conditional distributions are also encouraged to match. The use of the reverse Kullback-Leibler divergence (in contrast to the maximum likelihood setup) is  motivated by examining the following decomposition

\begin{align}\label{eq:KL}
\mathbb{KL}[p_{\theta}(\bm{x},t,\bm{u})||q(\bm{x},t,\bm{u})] & = -h(p_{\theta}(\bm{x},t,\bm{u}))) - \mathbb{E}_{p_{\theta}(\bm{x},t,\bm{u})}[\log(q(\bm{x},t,\bm{u}))],
\end{align}
where $h(p_{\theta}(\bm{x},t,\bm{u}))$ denotes the entropy of the generative model. The second term can be further decomposed as

\begin{align}
\mathbb{E}_{p_{\theta}(\bm{x},t,\bm{u})}[\log(q(\bm{x},t,\bm{u}))]  = & \int_{\mathcal{S}_{p_{\theta}\cap \mathcal{S}_q}} \log(q(\bm{x},t,\bm{u})) p_{\theta}(\bm{x},t,\bm{u}) d\bm{x} dt d\bm{u} \ + \label{eq:decomposition} \\
& \int_{\mathcal{S}_{p_{\theta}\cap \mathcal{S}_q^{o}}} \log(q(\bm{x},t,\bm{u})) p_{\theta}(\bm{x},t,\bm{u}) d\bm{x} dt d\bm{u} \nonumber,
\end{align}
where $\mathcal{S}_{p_{\theta}}$ and $\mathcal{S}_q$ denote the support of the distributions $p_{\theta}(\bm{x},t,\bm{u})$ and $q(\bm{x},t,\bm{u})$, respectively, while $\mathcal{S}_q^{o}$ denotes the complement of $\mathcal{S}_q$. Notice that by minimizing the Kullback-Leibler divergence in equation \ref{eq:KL} we introduce a mechanism that is trying to balance the effect of two competing objectives. Specifically, maximization of the entropy term $h(p_{\theta}(\bm{x},t,\bm{u})))$ encourages $p_{\theta}(\bm{x},t,\bm{u})$ to spread over its support set as wide, while the second integral term in equation \ref{eq:decomposition} introduces a strong (negative) penalty when the support of $p_{\theta}(\bm{x},t,\bm{u})$ and $q(\bm{x},t,\bm{u})$ do not overlap. Hence, the support of $p_{\theta}(\bm{x},t,\bm{u})$ is encouraged to spread only up to the point that $\mathcal{S}_{p_{\theta}}\cap \mathcal{S}_{q^{o}} = \emptyset$, implying that $\mathcal{S}_{p_{\theta}}\subseteq \mathcal{S}_{q^{o}}$. When $\mathcal{S}_{p_{\theta}}\subset \mathcal{S}_{q^{o}}$ the pathological issue of ``mode-collapse" (commonly encountered in the training of generative adversarial networks \cite{goodfellow2014generative}) is manifested \cite{salimans2016improved}.
This issue is present if one seeks to directly minimize the reverse Kullback-Leibler objective in equation \ref{eq:KL} as this provides no control on the relative importance of the two
terms. As discussed in \cite{li2018learning},  we may rather minimize $-\lambda h(p_{\theta}(\bm{x},t,\bm{u}))) - \mathbb{E}_{p_{\theta}(\bm{x},t,\bm{u})}[\log(q(\bm{x},t,\bm{u}))]$, with $\lambda \ge 1$ to allow for control of how much
emphasis is placed on mitigating mode collapse. It is then clear that the entropic regularization introduced by $h(p_{\theta}(\bm{x},t,\bm{u})))$ provides an effective mechanism for controlling and mitigating the effect of mode collapse, and, therefore, potentially enhancing the robustness adversarial inference procedures for learning $p_{\theta}(\bm{x},t,\bm{u})$.

Minimization of equation \ref{eq:KL} with respect to the generative model parameters $\theta$ presents two fundamental difficulties. First, the evaluation of both distributions $p_{\theta}(\bm{x},t,\bm{u})$ and $q(\bm{x},t,\bm{u})$ typically involves intractable integrals in high dimensions, and we may only have samples drawn from the two distributions, not their explicit analytical forms. Second, the differential entropy term $h(p_{\theta}(\bm{x},t,\bm{u})))$ is intractable as $p_{\theta}(\bm{x},t,\bm{u}))$ is not known a-priori. In the next sections we revisit the unsupervised formulation put forth in \cite{li2018learning} and derive a tractable inference procedure for learning $p_{\theta}(\bm{x},t,\bm{u}))$ from scattered observation pairs of $\bm{u}(\bm{x},t)$, namely $\{(\bm{x}_{i}, t_i), \bm{u}_i\}$, $i = 1,\dots,N_u$.

\subsubsection{Density ratio estimation by probabilistic classification}\label{sec:density_ratio}

By definition, the computation of the Kullback-Leibler divergence in equation \ref{eq:KL} involves computing an expectation over a log-density ratio, i.e.
$$\mathbb{KL}[p_{\theta}(\bm{x},t,\bm{u})||q(\bm{x},t,\bm{u})] := \mathbb{E}_{p_{\theta}(\bm{x},t,\bm{u})}\left[\log\left(\frac{p_{\theta}(\bm{x},t,\bm{u})}{q(\bm{x},t,\bm{u})}\right)\right].$$
In general, given samples from two distributions, we can approximate their density ratio by constructing a binary classifier that distinguishes between samples from the two distributions. To this end, we assume that  $N$ data points are drawn from $p_{\theta}(\bm{x},t,\bm{u})$ and are assigned a label $y=+1$. Similarly, we assume that $N$ samples are drawn from  $q(\bm{x},t,\bm{u})$ and assigned label $y=-1$. Consequently, we can write these probabilities in a conditional form, namely
$$p_{\theta}(\bm{x},t,\bm{u}) = \rho(\bm{x},t,\bm{u}|y=+1), \ \ q(\bm{x},t,\bm{u}) = \rho(\bm{x},t,\bm{u}|y=-1),$$
where $\rho(\bm{x},t,\bm{u}|y=+1)$ and $\rho(\bm{x},t,\bm{u}|y=-1)$ are the class probabilities predicted by a binary classifier $T(\bm{x},t,\bm{u})$.
Using Bayes rule, it is then straightforward to show that the density ratio of $p_{\theta}(\bm{x},t,\bm{u})$ and $q(\bm{x},t,\bm{u})$ can be computed as

\begin{align}
    \frac{p_{\theta}(\bm{x},t,\bm{u})}{ q(\bm{x},t,\bm{u})} & = \frac{\rho(\bm{x},t,\bm{u}|y=+1)}{\rho(\bm{x},t,\bm{u}|y=-1)} \nonumber\\
    & = \frac{\rho(y=+1|\bm{x},t,\bm{u})\rho(\bm{x},t,\bm{u})}{\rho(y=+1)} \bigg/ \frac{\rho(y=-1|\bm{x},t,\bm{u})\rho(\bm{x},t,\bm{u})}{\rho(y=-1)} \nonumber\\
    & = \frac{\rho(y=+1|\bm{x},t,\bm{u})}{\rho(y=-1|\bm{x},t,\bm{u})} = \frac{\rho(y=+1|\bm{x},t,\bm{u})}{1 - \rho(y=+1|\bm{x},t,\bm{u})} \nonumber\\
    & = \frac{T(\bm{x},t,\bm{u})}{1-T(\bm{x},t,\bm{u})}.
\end{align}
This simple procedure suggests that we can harness the power of deep neural network classifiers to obtain accurate estimates of the reverse Kullback-Leibler divergence in equation \ref{eq:KL} directly from data and without the need to assume any specific parametrization for the generative model distribution $p_{\theta}(\bm{x},t,\bm{u})$.

\subsubsection{Entropic regularization bound}\label{sec:entropy_bound}
Here we follow the derivation of Li {\em et. al} \cite{li2018learning} to construct a computable lower bound for the entropy $h(p_{\theta}(\bm{x},t,\bm{u}))$. To this end, we start by considering random variables $(\bm{x}, t, \bm{u}, \bm{z})$ under the joint distribution
$$p_{\theta}(\bm{x}, t, \bm{u}, \bm{z}) = p_{\theta}(\bm{u}, \bm{x}, t|\bm{z}) p(\bm{z}) = p_{\theta}(\bm{u}|\bm{x}, t, \bm{z}) p(\bm{x}, t, \bm{z}),$$
where $p_{\theta}(\bm{u}|\bm{x}, t, \bm{z}) = \delta(\bm{u}-f_{\theta}(\bm{x}, t, \bm{z}))$, and $\delta(\cdot)$ is the Dirac delta function. The mutual information between $(\bm{x}, t, \bm{u})$ and $\bm{z}$ satisfies the information theoretic identity
$$ I(\bm{x}, t, \bm{u}; \bm{z}) = h(\bm{x}, t, \bm{u})-h(\bm{x}, t, \bm{u}|\bm{z}) =
h(\bm{z}) - h(\bm{z}|\bm{x}, t, \bm{u}),$$
where $h(\bm{x}, t, \bm{u})$, $h(\bm{z})$ are the marginal entropies and $h(\bm{x}, t, \bm{u}|\bm{z})$, $h(\bm{z}|\bm{x}, t, \bm{u})$ are the conditional entropies \cite{akaike1998information}.
Since in our setup $\bm{x}$ and $t$ are deterministic variables independent of $\bm{z}$, and samples of $p_{\theta}(\bm{u}|\bm{x}, t, \bm{z})$ are generated by a deterministic function $f_{\theta}(\bm{x}, t, \bm{z})$, it follows that $h(\bm{x}, t, \bm{u}|\bm{z}) = 0$.
We therefore have
\begin{equation}\label{eq:info_identity}
h(\bm{x}, t, \bm{u}) = h(\bm{z}) - h(\bm{z}|\bm{x}, t, \bm{u}),
\end{equation}
where $h(\bm{z}) := -\int \log p(\bm{z}) p(\bm{z}) d\bm{z}$ does not depend on the generative model parameters $\theta$.

Now consider a general variational distribution $q_{\phi}(\bm{z}|\bm{x}, t, \bm{u})$ parametrized by a set of parameters $\phi$. Then,

\begin{align}
   h(\bm{z}|\bm{x}, t, \bm{u})  = & -\mathbb{E}_{p_{\theta}(\bm{x}, t, \bm{u}, \bm{z})}[\log(p_{\theta}(\bm{z}|\bm{x}, t, \bm{u}))] \nonumber\\
    = & -\mathbb{E}_{p_{\theta}(\bm{x}, t, \bm{u}, \bm{z})}[\log(q_{\phi}(\bm{z}|\bm{x}, t, \bm{u}))] \nonumber  \\
    & -\mathbb{E}_{p_{\theta}(\bm{x}, t, \bm{u})}[\mathbb{KL}[p_{\theta}(\bm{z}|\bm{x}, t, \bm{u})||q_{\phi}(\bm{z}|\bm{x}, t, \bm{u})]] \nonumber \\
    \le & -\mathbb{E}_{p_{\theta}(\bm{x}, t, \bm{u}, \bm{z})}[\log(q_{\phi}(\bm{z}|\bm{x}, t, \bm{u}))]. \label{eq:entropy_bound}
\end{align}
Viewing $\bm{z}$ as a set of latent variables, then $q_{\phi}(\bm{z}|\bm{x}, t, \bm{u})$ is a variational approximation to the true intractable posterior over the latent variables $p_{\theta}(\bm{z}|\bm{x}, t, \bm{u})$. Therefore,
if $q_{\phi}(\bm{z}|\bm{x}, t, \bm{u})$  is introduced as an auxiliary inference model associated with the generative model $p_{\theta}(\bm{x}, t, \bm{u})$, for which $\bm{u} = f_{\theta}(\bm{x}, t,\bm{z})$ and $\bm{z}\sim p(\bm{z})$, then we can use equations \ref{eq:info_identity} and \ref{eq:entropy_bound} to bound the entropy term in equation \ref{eq:KL} as
\begin{equation}
    h(p_{\theta}(\bm{x},t,\bm{u})) \ge h(p(\bm{z})) + \mathbb{E}_{p_{\theta}(\bm{x}, t, \bm{u}, \bm{z})}[\log(q_{\phi}(\bm{z}|\bm{x}, t, \bm{u}))].
\end{equation}
Note that the inference model $q_{\phi}(\bm{z}|\bm{x}, t, \bm{u})$ plays the role of a variational approximation to the true posterior over the latent variables, and appears naturally using information theoretic arguments in the derivation of the lower bound.

\subsubsection{Adversarial training objective}

By leveraging the density ratio estimation procedure described in section \ref{sec:density_ratio} and the entropy bound derived in section \ref{sec:entropy_bound}, we can derive the following loss functions for minimizing the reverse Kullback-Leibler divergence with entropy regularization

\begin{align}
	\mathcal{L}_{\mathcal{D}}(\psi) = & \  \mathbb{E}_{q(\bm{x},t)p(\bm{z})}[\log\sigma(T_{\psi}(\bm{x},t,f_{\theta}(\bm{x},t,\bm{z})))] + \nonumber \\ & \ \mathbb{E}_{q(\bm{x},t,\bm{u})}[\log(1-\sigma(T_{\psi}(\bm{x},t,\bm{u})))] \label{eq:discriminator_loss}\\
	\mathcal{L}_{\mathcal{G}}(\theta, \phi) = & \ \mathbb{E}_{q(\bm{x},t)p(\bm{z})}[T_{\psi}(\bm{x}, t, f_{\theta}(\bm{x},t,\bm{z}))+ (1-\lambda)\log(q_{\phi}(\bm{z}|\bm{x},t,f_{\theta}(\bm{x},t,\bm{z})))] \label{eq:generator_loss},
\end{align}
where $\sigma(x)=1/(1+e^{-x})$ is the logistic sigmoid function. Notice how the binary cross-entropy objective of equation \ref{eq:discriminator_loss} aims to progressively improve the ability of the classifier $T_{\psi}(\bm{x},t,\bm{u})$ to discriminate between ``fake" samples $(\bm{x},t,f_{\theta}(\bm{x},t,\bm{z}))$ obtained from the generative model $p_{\theta}(\bm{x},t,\bm{u})$ and ``true" samples $(\bm{x},t,\bm{u})$ originating from the observed data distribution $q(\bm{x},t,\bm{u})$. Simultaneously, the objective of equation \ref{eq:generator_loss} aims at improving the ability of the generator $f_{\theta}(\bm{x},t,\bm{u})$ to generate increasingly more realistic samples that can ``fool" the discriminator $T_{\psi}(\bm{x},t,\bm{u})$. Moreover, the encoder $q_{\phi}(\bm{z}|\bm{x},t,f_{\theta}(\bm{x},t,\bm{z}))$ not only serves as an entropic regularization term than allows us to stabilize model training and mitigate the pathology of mode collapse, but also provides a variational approximation to true posterior over the latent variables. The way it naturally appears in the objective of equation \ref{eq:generator_loss} also encourages the cycle-consistency of the latent variables $\bm{z}$; a process that is known to result in disentangled and interpretable low-dimensional representations of the observed data \cite{friedman2001elements}, which could be subsequently used as good features for nonlinear model order reduction.

In theory, the optimal set of parameters $\{\theta^{\ast}, \phi^{\ast}, \psi^{\ast}\}$ correspond to the Nash equilibrium of the two player game defined by the loss functions in equations \ref{eq:discriminator_loss},\ref{eq:generator_loss}, for which one can show that the exact model distribution and the exact posterior over the latent variables can be recovered \cite{goodfellow2014generative, pu2017symmetric}. In practice, although there is no guarantee that this optimal solution can be attained, the generative model can be trained by alternating between optimizing the two objectives in equations \ref{eq:discriminator_loss},\ref{eq:generator_loss} using stochastic gradient descent as

\begin{align}
    & \mathop{\max}_{\psi} \ \mathcal{L}_{\mathcal{D}}(\psi) \label{eq:discriminator}\\
	& \mathop{\min}_{\theta, \phi} \ \mathcal{L}_{\mathcal{G}}(\theta, \phi) \label{eq:generator}.
\end{align}

\subsubsection{Adversarial training with physics-informed constrains}

In order to learn the physics-informed probabilistic model of equation \ref{eq:PSC} from data we can extend the adversarial  inference framework presented above by appropriately penalizing the loss function of the generator (see equation \ref{eq:generator_loss}). The available data correspond to scattered observation pairs $\{(\bm{x}_{i}, t_i), \bm{u}_i\}$, $i = 1,\dots,N_u$, originating from known initial or boundary conditions, or any other (potentially noisy) measurements of $\bm{u}(\bm{x},t)$. In analogy to the deterministic setting put for in \cite{raissi2017physics1} and summarized in section \ref{sec:PINNs}, by defining $\bm{r}_{\theta}(\bm{x},t,\bm{z}):= \frac{\partial}{\partial t} f_{\theta}(\bm{x},t,\bm{z})  + \mathcal{N}_{\bm{x}}f_{\theta}(\bm{x},t,\bm{z})$ we essentially introduce a new conditional probability model $p_{\theta}(\bm{r}|\bm{x},t,\bm{z})$ that shares the same parameters as $p_{\theta}(\bm{u}|\bm{x},t,\bm{z})$, albeit the underlying neural network that serves as its approximation has different activation functions. However, since we would like to encourage every sample $\bm{u} = f_{\theta}(\bm{x},t,\bm{z})$ produced by the generator to satisfy the PDE constraint, we can simply treat the residual as a deterministic variable, i.e,  $\bm{r}_{\theta}(\bm{x},t,\bm{z}) = \bm{r}_{\theta}(\bm{x},t)$, and enforce the constraint at a finite set of collocation points $N_r$ by simply minimizing the mean square loss
\begin{equation}
\label{eq:PI_loss}
	\mathcal{L}_{PDE}(\theta) = \frac{1}{N_r}\sum_{i=1}^{N_r}\|\bm{r}_{\theta}(\bm{x}_i, t_i) - \bm{r}_i\|^2.
\end{equation}
Then, the resulting adversarial game for training the physics-informed model of equation \ref{eq:PSC} takes the form

\begin{equation}
    \label{eq:objective}
    \begin{aligned}
    	& \mathop{\max}_{\psi} \ \mathcal{L}_{\mathcal{D}}(\psi)\\
    	& \mathop{\min}_{\theta, \phi} \ \mathcal{L}_{\mathcal{G}}(\theta, \phi) + \beta \mathcal{L}_{PDE}(\theta),
    \end{aligned}
\end{equation}
where positive values of $\beta$ can be selected to place more emphasis on penalizing the PDE residual. For $\beta>0$, the residual loss $\mathcal{L}_{PDE}(\theta)$ acts as a regularization term that approximately enforces the given physical constraint, and, therefore, encourages the generator $p_{\theta}(\bm{u}|\bm{x},t,\bm{z})$ to produce samples that satisfy the underlying partial differential equation. Also note that this structured approach also encourages the encoder $q_{\phi}(\bm{z}|\bm{x},t,f_{\theta}(\bm{x},t,\bm{z}))$ to learn a set of spatio-temporal latent variables $\bm{z}$ that are relevant to the underlying physics, possibly opening a new directions for probabilistic model order reduction of complex systems.

\subsubsection{Predictive distribution}

Once the model is trained we can construct a probabilistic ensemble for the solution $\bm{u}(\bm{u}|\bm{x}, t, \bm{z})$ by sampling latent variables from the prior $p(\bm{z})$ and passing them through the generator to yield samples $\bm{u} = f_{\theta}(\bm{x}, t, \bm{z})$ that are distributed according to the predictive model distribution $p_{\theta}(\bm{u}|\bm{x}, t, \bm{z})$. Note that although the explicit form of this  distribution is not known, we can efficiently compute any of its moments via Monte Carlo sampling. The cost of this prediction step is negligible compared to the cost of training the model, as it only involves a single forward pass through the generator function $f_{\theta}(\bm{x}, t, \bm{z})$. Typically, we compute the mean and variance of the predictive distribution at a new test point $(\bm{x}^{\ast}, t^{\ast})$ as

\begin{align}
    \mu_{\bm{u}}(\bm{x}^{\ast}, t^{\ast}) & = \mathbb{E}_{p_{\theta}}[\bm{u}|\bm{x}^{\ast}, t^{\ast}, \bm{z}] \approx \frac{1}{N_s}\sum\limits_{i=1}^{N_s} f_{\theta}(\bm{x}^{\ast}, t^{\ast}, \bm{z}_i), \label{eq:predictive_mean} \\
    \sigma^{2}_{\bm{u}}(\bm{x}^{\ast}, t^{\ast}) & = \mathbb{V}\text{ar}_{p_{\theta}}[\bm{u}|\bm{x}^{\ast}, t^{\ast}, \bm{z}] \approx \frac{1}{N_s}\sum\limits_{i=1}^{N_s} [f_{\theta}(\bm{x}^{\ast}, t^{\ast}, \bm{z}_i) - \mu_{\bm{u}}(\bm{x}^{\ast}, t^{\ast})]^2, \label{eq:predictive_variance}
\end{align}
where $\bm{z}_i \sim p(\bm{z})$, $i = 1,\dots,N_s$, and $N_s$ corresponds to the total number of Monte Carlo samples.

\subsubsection{Advantages and caveats of adversarial learning}

Since their recent introduction \cite{goodfellow2014generative, makhzani2015adversarial, dumoulin2016adversarially, mescheder2017adversarial}, adversarial learning  techniques have provided great flexibility for performing probabilistic computations with arbitrarily complex implicit distributions. Essentially, they have lifted  the over-simplified approximations typically used in variational inference (Gaussian approximations, exponential families, etc.) \cite{blei2017variational}, yielding very general and flexible schemes for statistical inference. However, this flexibility comes at a price, as such methods in practice require very careful tuning in order to achieve stable and accurate performance. To this end, recall the training objective defined in equation \ref{eq:objective} that introduces an adversarial game between the generator and discriminator networks \cite{goodfellow2014generative}. In practice, this mini-max optimization problem is solved by alternating stochastic gradient updates between the two competing objectives, and it is highly sensitive on the capacity of the neural networks modeling the generator and discriminator, as well as the relative frequency with which the parameters of each network are updated within each iteration of stochastic gradient descent. To this end, we provide a series of empirical observations and lessons we learned throughout this study that can enhance the robustness and stability of this training procedure:

\begin{itemize}
    \item Changing the relative number of stochastic gradient updates for the generator $K_g$ and the discriminator $K_d$ is equivalent to changing their neural network architecture. For example, we can reduce the capacity of discriminator by either performing more stochastic gradient updates for the generator, or remove one layer in the neural network architecture of the discriminator.
    \item Given enough collocation points $N_r$ for penalizing the PDE residual, we can obtain robust uncertainty estimates together with precise predictions simply by tuning the capacity of discriminator and generator networks.
    \item Typically, by fixing the generator, we expect the discriminator to have some capacity so that the model training dynamics remain stable. But, we do not want the discriminator to be very powerful as in that case there will by very little information from the discriminator that can help the generator to improve towards producing more realistic samples (this a common characteristic of adversarial inference procedures \cite{goodfellow2014generative}).
    \item  For cases with a small number of training data we should reduce the capacity of the discriminator.  This can be achieved by either changing the relative frequency of stochastic gradient updates for the generator and discriminator, or by reducing the capacity of the discriminator neural network architecture.
\end{itemize}

\section{Results}
\label{sec:Results}

In all examples we have trained the models for 30,000 stochastic gradient descent steps using the Adam optimizer \cite{kingma2014adam} with a learning rate of $10^{-4}$, while fixing a one-to-five ratio for the discriminator versus generator updates. Moreover, we have fixed the entropic regularization and the residual penalty parameters to $\lambda = 1.5$ and $\beta = 1.0$, respectively. The proposed algorithms were implemented in Tensorflow v1.10 \cite{abadi2016tensorflow}, and computations were performed in single precision arithmetic on a single NVIDIA Tesla P100 GPU card. All data and code accompanying this manuscript will be made available at \url{https://github.com/PredictiveIntelligenceLab/UQPINNs}.

\subsection{A pedagogical example}

Let us illustrate the basic capabilities of the proposed methods through a simple example corresponding to the following nonlinear second-order ordinary differential equation
\begin{equation}\label{eq:ODE}
\begin{aligned}
	&u_{xx} - u^2u_x = f(x), \quad\quad x\in[-1, 1], \\
	&f(x) = -\pi^2\sin(\pi x) - \pi\cos(\pi x)\sin^2(\pi x),\\
\end{aligned}
\end{equation}
subject to random boundary conditions $u(-1), u(1)\sim \mathcal{N}(\bm{0}, \sigma_n^2\bm{I})$.
For this simple example, the deterministic solution corresponding to $\sigma_n^2 = 0$ can be readily obtained as $u(x) = \sin(\pi x)$. Given $N_u$ observations of $u(x)$ corresponding to different realizations of the random boundary conditions our goal is to obtain a probabilistic representation of the solution $p_{\theta}(u|x,\bm{z})$ by training a physics-informed generative model of the form $u = f_{\theta}(x,\bm{z})$, $\bm{z}\sim p(\bm{z})$ that is constrained by equation \ref{eq:ODE}. To this end, we introduce three deterministic mappings parametrized by deep neural networks, namely $f_{\theta}(x, \bm{z})$, $q_{\phi}(x, u)$, and $T_{\psi}(x, u)$ corresponding to the generator, encoder, and discriminator functions introduced in section \ref{sec:ADVI}. By construction, we also obtain a physics-informed neural network $r_{\theta}(x)$ corresponding to the deterministic residual of equation \ref{eq:ODE} that will be used to approximately enforce the differential equation constraint at a set of $N_r = 100$ randomly distributed collocation points $x\in[-1, 1]$. All neural networks were chosen to have two hidden layers with 50 neurons in each layer, and a hyperbolic tangent activation function. Moreover, the dimensionality of the latent variables was set to one, i.e. $\bm{z}=z$, and we have assumed an isotropic standard normal prior, namely $p(z)\sim \mathcal{N}(0,1)$. As the training data for $u(x)$ reflects the uncertainty in the boundary conditions, the role of the latent variables $z$ is to enable the propagation of  this uncertainty into the predicted solution obtained through the generative model $p_{\theta}(u|x,z)$.

\begin{figure}
\centering
\includegraphics[width=\textwidth]{ode_solution.pdf}
\caption{{\em A pedagogical example:} (a) Mean and two standard deviations of $p_{\theta}(u|x,z)$ against the exact solution for deterministic boundary data. (b) Mean and two standard deviations of $p_{\theta}(u|x,z)$ against the reference Monte Carlo solution for random boundary data corresponding to $5\%$ Gaussian uncorrelated noise.}
\label{fig:ODE_solution}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=\textwidth]{ode_density.pdf}
\caption{{\em A pedagogical example:} Predicted marginal densities against the reference Monte Carlo solution. (a) $p_{\theta}(u|x=-0.5,z)$. (b) $p_{\theta}(u|x=+0.5,z)$.}
\label{fig:ODE_density}
\end{figure}

Here we have considered two cases corresponding to deterministic and random boundary conditions, namely (i) $\sigma_n^2=0$ (i.e., noise-free data), and (ii) $\sigma_n^2= 0.05$ (i.e., $\%5$ Gaussian uncorrelated noise). In all cases, the training data consists of $N_u=\textcolor{black}{20}$ realizations for each boundary point, $u(-1), u(1)$, and a total of $N_r=100$ collocation points for enforcing the residual of equation \ref{eq:ODE}. Our probabilistic predictions for this example are summarized in figures \ref{fig:ODE_solution} and \ref{fig:ODE_density}. Specifically, figure \ref{fig:ODE_solution}(a) shows the generative model predictive mean and two standard deviations, plotted against the exact solution of this problem. Note that this case corresponds to deterministic training data for the boundary conditions, hence the exact solution is deterministic, and the prediction error here is measured as $\mathcal{E}_{\mathbb{L}_2} = 1.36\cdot10^{-3}$ in the relative $\mathbb{L}_2$ norm

\begin{equation}
    \mathcal{E}_{\mathbb{L}_2} := \frac{\sqrt{\sum\limits_{i=1}^{N^{\ast}} [\mu_{\bm{u}}(x_i^{\ast}) - u(x_i^{\ast})]^2}}{\sqrt{\sum\limits_{i=1}^{N^{\ast}} [u(x_i^{\ast})]^2}},
\end{equation}
where $N^{\ast} = \textcolor{black}{200}$ denotes the total number of equidistant test points $x^{\ast}$ in the interval $[-1, 1]$. Moreover, the variance shown in the inset of figure \ref{fig:ODE_solution}(a) serves as an {\it a-posteriori} error estimate  that quantifies the uncertainty associated with the generative model predictions. Figure \ref{fig:ODE_solution}(b) shows the resulting prediction and uncertainty estimates corresponding to random boundary conditions, compared against a reference mean solution obtained numerically using a spectral method with \textcolor{black}{2,000} Monte Carlo samples. In this case, the predictive uncertainty of the generative model reflects the aggregate {\em total} uncertainty due to both randomness in the boundary conditions and the inherent epistemic uncertainty in the neural network approximation. As the generative model can return a complete statistical characterization of the solution by means of its conditional probability density $p_{\theta}(u|x,z)$, in figure \ref{fig:ODE_density} we provide a visual comparison of the one-dimensional marginals between our predictions and the reference Monte Carlo solution corresponding to the spatial locations $x = -0.5$ and $x = 0.5$.

Albeit simple, this example aims to demonstrate the basic capabilities of the proposed methodology in propagating uncertainty through non-linear partial differential equations. In contrast to previous approaches to inferring solutions of partial differential equations from data \cite{cockayne2016probabilistic,raissi2017inferring,raissi2018numerical,cockayne2017bayesian}, the proposed methodology does not rely on Gaussian assumptions, and it can directly tackle nonlinear problems without any need for linearization.

\subsection{Burgers equation}\label{sec:burgers}
In this example we aim to provide a comprehensive systematic study to quantify the robustness of the proposed methods with respect to different parameter choices. We will do so through the lens of a more challenging canonical problem involving the non-linear time-dependent Burgers equation in one spatial dimension:
\begin{equation}
\label{eq:Burgers}
\begin{aligned}
	&u_t + u u_x - \nu u_{xx} = 0, \quad\quad x\in[-1, 1], t\in[0,1],\\
	&u(0,x) = -\sin(\pi x),\\
	&u(t,-1) = u(t,1) = 0,\\
\end{aligned}
\end{equation}
where the viscosity parameter is chosen as $\nu = (0.01/\pi)$ in order to generate a strongly nonlinear response that leads to the development of shock discontinuities in finite time. This is one of the few nonlinear partial differential equations that admits an exact solution through the Cole-Hopf transformation \cite{hopf1950partial}; a solution that will be subsequently used to test the validity of our predictions.

Here we represent the unknown solution $u(x,t)$ using a physics-informed generative model of the form $u = f_{\theta}(x, \textcolor{black}{t}, \bm{z})$, and we will introduce parametric functions corresponding to a generator $f_{\theta}(x, t, \bm{z})$, an encoder $q_{\phi}(x, t, u)$, and a discriminator $T_{\psi}(x, \textcolor{black}{t}, u)$ all constructed using deep feed-forward neural networks. The baseline architectures for the generator and the encoder have 4 hidden layers with 50 neurons per layers, while the discriminator network has 3 hidden layers and 50 neurons per layer. The activation function in all cases is chosen to be a hyperbolic tangent non-linearity. The prior over the latent variables $p(\bm{z})$ is chosen again to be a one-dimensional isotropic Gaussian distribution, i.e. $\bm{z} = z$, $z \sim \mathcal{N}(0,1)$.

First we consider a baseline scenario, in which we train our probabilistic model using a data-set comprising of $N_u = \textcolor{black}{1}50$ noisy-free input/output pairs for $u(x,t)$ -- 50 points for the initial condition (see figure \ref{fig:Burgers_initial}(a)) and \textcolor{black}{50} points for each of the domain boundaries -- plus an additional  $N_r = 10,000$ collocation points for enforcing the residual of the Burgers equation using the loss of equation  \ref{eq:PI_loss}. All data points were randomly selected within the bounds given in equation \ref{eq:Burgers}. The result of this experiment is summarized in \ref{fig:Burgers_baseline} where we report the predicted mean solution, as well as the uncertainty associated with this prediction as quantified by two standard deviations of the generative model $p_{\theta}(u|x,t,z)$. As the training data for this case is noise-free, the solution to this problem is deterministic, and the resulting uncertainty captured in $p_{\theta}(u|x,t,z)$ can be viewed as an {\em a-posteriori} error estimate of  the neural network approximation error due to the finite number of training data, which is measured as $\mathcal{E}_{\mathbb{L}_2} = 4.1\cdot 10^{-2}$ in the relative $\mathbb{L}_{2}$ norm. As discussed in \cite{raissi2017physicsI}, a higher approximation accuracy can be achieved by training the generative model using a quasi-Newton optimizer (e.g. L-BFGS \cite{liu1989limited}), however here we chose to use stochastic gradient descent using Adam updates \cite{kingma2014adam} in order to highlight the ability of the proposed method to return uncertainty estimates when the model predictions are not perfectly accurate.

\begin{figure}
\centering
\includegraphics[width=\textwidth]{burgers_initial.pdf}
\caption{{\em Burgers equation:} (a) Exact initial condition and noise-free training data (50 points). (b) Training data corresponding to a single realization of the non-additive noise corruption process (100 points, generated by equation \ref{eq:noisy_data}).}
\label{fig:Burgers_initial}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=\textwidth]{NewData1.pdf}
\includegraphics[width=\textwidth]{NewError1.pdf}
\caption{{\em Burgers equation with noise-free data:} {\it Top:} Mean of $p_{\theta}(u|x,t,z)$, along with the location of the noisy training data $\{(x_{i}, t_i), u_i\}$, $i=1,\dots,N_u$. {\it Middle:} Prediction and predictive uncertainty at $t=0.25$, $t=0.5$ and $t=0.75$. {\it Bottom:} Variance of $p_{\theta}(u|x,t,z)$.}
\label{fig:Burgers_baseline}
\end{figure}

Second, we repeat the same test for a more complicated scenario in which the initial condition has been corrupted by non-additive, non-Gaussian noise as shown in figure \ref{fig:Burgers_initial}(b), where the noise variance is larger around $x=0$, therefore amplifying the effect of uncertainty on the shock formation. Here the neural network architecture as well as the number and location of training points have been kept fixed as described above, but the initial condition is now corrupted as
\begin{equation}
\label{eq:noisy_data}
u(x, 0) = -\sin(\pi(x+2\delta))+\delta, \qquad \delta = \frac{\epsilon}{\exp(3|x|)}, \qquad \epsilon\sim N(0, 0.1^2).
\end{equation}
The results of this experiment are summarized in figure \ref{fig:Burgers_noisy}. We observe that the resulting generative model $p_{\theta}(u|x,t,z)$ can effectively capture the uncertainty in the resulting spatio-temporal solution due to the propagation of the input noise process through the complex non-linear dynamics of the Burgers equation. As expected, the uncertainty concentrates around the shock. Although we only plot the first two moments of the solution, we must emphasize that the generative model $p_{\theta}(u|x,t,z)$ provides a complete probabilistic characterization of its non-Gaussian statistics.

\begin{figure}
\centering
\includegraphics[width=\textwidth]{NewData3.pdf}
\includegraphics[width=\textwidth]{NewError3.pdf}
\caption{{\em Burgers equation with noisy data:} {\it Top:} Mean of $p_{\theta}(u|x,t,z)$, along with the location of the training data $\{(x_{i}, t_i), u_i\}$, $i=1,\dots,N_u$. {\it Middle:} Prediction and predictive uncertainty at $t=0.25$, $t=0.5$ and $t=0.75$. {\it Bottom:} Variance of $p_{\theta}(u|x,t,z)$.}
\label{fig:Burgers_noisy}
\end{figure}

In order to further investigate the performance of the proposed methodology for different parameter settings, we have performed a series of comprehensive systematic studies that aim to quantify the sensitivity of the resulting predictions on: (i) the neural network initialization, (ii) the total number of training and collocation points, (iii) the neural network architecture, and (iv) the adversarial training procedure. The results of these systematic studies are provided in \ref{sec:appendix}.

\subsection{Discovery of constitutive laws for flow through porous media}
In our final example, we aim to demonstrate the ability of the proposed methods to discover unknown  constitutive relationships directly from data with quantified uncertainty. To this end, we revisit the Darcy flow example put forth in \cite{tartakovsky2018learning} corresponding to a two-dimensional nonlinear diffusion equation with an unknown state-dependent diffusion coefficient
\begin{equation}
\label{eq:diffusion}
\begin{aligned}
	\nabla_{\bm{x}}\cdot[K(u)\nabla_{\bm{x}} u(\bm{x})] &= 0, \quad\quad \bm{x}=(x_1, x_2)\in \Omega = (0, L_1)\times (0, L_2)\\
	u(\bm{x}) &= u_0, \quad\quad\quad\;\; x_1 = L_1\\
	-K(u)\frac{\partial u(\bm{x})}{\partial x_1} &= q, \quad\quad\quad\quad x_1 = 0\\
	\frac{\partial u(\bm{x})}{\partial x_2} &= 0, \quad\quad\quad\quad x_2 = \{0, L_2\},\\
\end{aligned}
\end{equation}
where $q = 8.25 \times 10^{-5}$m/s and $u_0 = - 10$m are known boundary conditions.
In order to benchmark and validate our model predictions we consider a realistic data-set generated using the Subsurface Transport
Over Multiple Phases (STOMP) code \cite{white1995modeling} with the van Genuchten model
\cite{van1980closed} for $K(u)$ which reads as
\begin{equation}
\label{eq:STOMP}
\begin{aligned}
	K(s(u)) &= K_s s^{\frac{1}{2}}[1-(1-s^{\frac{1}{m}})^m]^2\\
	s(u) &= \{1+ [\alpha (u_g - u)]^{\frac{1}{1-m}}\}^{-m},\\
\end{aligned}
\end{equation}
with the following parameter values: $K_s = 8.25 \times 10^{-4}$ m/s, $u_g = 0$, $m = 0.469$, $\alpha = 0.1$, $L_1 = 10$ m and $L_2 = 10$ m.

Our goal is twofold: we aim to construct a physics-informed probabilistic model for $p_{\theta}(u|\bm{x},\bm{z})$, and simultaneously learn the unknown state-dependent diffusion coefficient $K(u)$ directly from data on $u(\bm{x})$ (i.e., we assume no measurements of $K(u)$). To this end, in addition to the three deterministic mappings $f_{\theta}(\bm{x}, \bm{z})$, $q_{\phi}(\bm{x}, u)$, and $T_{\psi}(\bm{x}, u)$ corresponding to the generator, encoder, and discriminator described in section \ref{sec:ADVI}, here we also introduce
another neural network $f_{\gamma}(u)$ for approximating $K(u)$. The parameters of $f_{\gamma}(u)$ are essentially inherited by the physics-informed residual network $r_{\theta,\gamma}(\bm{x}) := \nabla_{\bm{x}}\cdot[f_{\gamma}(f_{\theta}(\bm{x}, \bm{z}))\nabla_{\bm{x}} f_{\theta}(\bm{x}, \bm{z}))] $ that aims to enforce the residual of equation \ref{eq:diffusion} at the $N_r$ collocation points for any set of latent variables $\bm{z}$. All neural networks are chosen to have 2 hidden layers with 50 neurons per each, and a hyperbolic tangent activation function, while the probabilistic model for $p_{\theta}(u|\bm{x},\bm{z})$ assumes a two dimensional latent space with an isotropic Gaussian prior, i.e., $\bm{z} \sim \mathcal{N}(\bm{0}, \bm{I})$.

By construction, our probabilistic model for $p_{\theta}(u|\bm{x},\bm{z})$ can return predictions of the unknown solution $u(\bm{x})$ with quantified uncertainty.
We can then use this model to propagate uncertainty in our predictions of $K(u)$ via Monte Carlo sampling. Specifically, once the model is trained end-to-end, we can easily generate samples of $u(\bm{x})$ from $p_{\theta}(u|\bm{x},\bm{z})$ and propagate them through $f_{\gamma}(u)$ to obtain a samples for $K(u)$. Essentially this results in an implicit generative model $p_{\theta,\gamma}(k|\bm{x},\bm{z})$ which can fully characterize uncertainty in our predictions of the unknown state-dependent diffusion coefficient.

Here we also have considered two cases corresponding to noise-free training data and noisy data corrupted by $5\%$ Gaussian uncorrelated noise. In the \textcolor{black}{noise-free} case we used $N_u = 600$ scattered measurements of the unknown solution $u(\bm{x})$ -- 200 inside the domain $\Omega$ and 100 on each one of the four boundaries -- and total number of $N_r = 10,000$ randomly selected collocation points inside the domain for penalizing the residual of equation \ref{eq:diffusion}. \textcolor{black}{For the noisy case we chose $N_u = 1,400$ scattered measurements of the unknown solution $u(\bm{x})$ -- 1,000 inside the domain $\Omega$ and 100 on each one of the four boundaries -- while still keeping $N_r = 10,000$ collocation points.} Figure \ref{fig:k_u} summarizes the results for both cases by showing the predictive mean and two standard deviations of the corresponding generative model $p_{\theta,\gamma}(k|\bm{x},\bm{z})$, against the reference (deterministic) solution obtained from the Subsurface Transport Over Multiple Phases (STOMP) code \cite{white1995modeling} with the van Genuchten model \cite{van1980closed} (see equation \ref{eq:STOMP}). Evidently, the generative model is able to recover a sensible prediction for the unknown state-dependent diffusion coefficient with quantitative uncertainty, even when the training data on $u(\bm{x})$ is heavily corrupted by noise. Moreover, notice that  $k(u)$ implicitly depends on the spatial coordinates $\bm{x}=(x_1, x_2)$. In figure \ref{fig:k_x} we present the resulting prediction for $K(u(x_1,x_2))$ corresponding to the noise-free case, against the reference solution, as well as their point-wise absolute error.

Again we must emphasize that these predictions are obtained without ever observing any data on $K(u)$, while they are accompanied by quantitative estimates that jointly characterize the uncertainty due to noise in the training data for $u(\bm{x})$, and the underlying approximation error of the neural networks. This theme of consistently inferring correlated continuous quantities of interest from a small set of measurements by leveraging the underlying laws of physics is a great example of the exciting capabilities that physics informed machine learning has to offer.

\begin{figure}
\centering
\includegraphics[width=\textwidth]{k_u.pdf}
\caption{Prediction with quantified uncertainty of unknown state-dependent diffusion coefficient $K(u)$ compared against the reference solution obtained from the Subsurface Transport
Over Multiple Phases (STOMP) code \cite{white1995modeling} with the van Genuchten model
\cite{van1980closed} (see equation \ref{eq:STOMP}). (a) Noise-free training data for $u(\bm{x})$. (b) Noisy training data for $u(\bm{x})$ with noise level of $5\%$.}
\label{fig:k_u}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=\textwidth]{k_x.pdf}
\caption{Prediction of unknown state-dependent diffusion coefficient $K(u(x_1,x_2))$. (a) Reference solution obtained from the Subsurface Transport
Over Multiple Phases (STOMP) code \cite{white1995modeling} with the van Genuchten model
\cite{van1980closed} (see equation \ref{eq:STOMP}). (b) Predictive mean of the generative model $p_{\theta,\gamma}(k|\bm{x},\bm{z})$ trained on noise-free data for $u(\bm{x})$. (c) Absolute point-wise prediction error.}
\label{fig:k_x}
\end{figure}

\section{Conclusions}

We presented a class of  probabilistic physics-informed neural networks that are capable of approximating arbitrary conditional probability densities, while being constrained to generate samples that approximately satisfy given partial differential equations. Moreover, we have derived a flexible regularized adversarial inference framework that enables the end-to-end training of such models directly from noisy and incomplete measurements. Uncertainty in the system inputs and/or outputs is captured through a set of latent variables that are relevant to the underlying physics, and could possibly open new directions for probabilistic model-order reduction of complex systems. These developments allow us to perform probabilistic computations for uncertain systems, train deep generative models in small data regimes, handle complex noise processes, and seamlessly carry out uncertainty propagation studies for physical systems without the need for repeated evaluation of experiments and numerical simulations.

Although the proposed adversarial inference framework provides great flexibility for performing probabilistic computations and approximating arbitrarily complex and high-dimensional probability distributions, it relies on carefully tuning the interplay between the generator and discriminator networks. This is a known limitation of adversarial algorithms, and, although several works have led to improvements \cite{tolstikhin2017wasserstein, arjovsky2017wasserstein}, it still largely remains an open research problem. An alternative path for enhancing the robustness of the inference procedure, while not compromising its ability to handle complex probability distributions, comes through the use of invertible transformations and flow-based generative models \cite{rezende2015variational, kingma2018glow}.
Future work will examine the applications of such models in the context of physics-informed neural networks with the goal of robustifying the proposed methods and scaling them to more realistic systems.

\section*{Acknowledgements}
This work received support from the US Department of Energy under the Advanced Scientific Computing Research program (grant DE-SC0019116) and the Defense Advanced Research Projects Agency under the Physics of Artificial Intelligence program. We would also like to thank Dr. Alexandre Tartakovsky from the Pacific Nortwest National Laboratory for providing the Darcy flow data-set.

\input{appendix.tex}


\section{Introduction} \label{sec:introduction}
Scene graph generation aims at localizing all pairs of interacting objects in an input image and identifying relationships between them in the form of graph representations, where the nodes represent objects, e.g.~\textit{a man}, \textit{a horse}, and the edges represent the pairwise interaction between objects, e.g.~man \textit{feeding} horse. In recent years, significant progress has been made in scene graph generation \cite{Xu2017, Zellers2017, Li2018, Newell2017, Li2017a}.
One approach is to integrate external knowledge via linguistic priors \cite{Li2017a, Lu2016, Yu2017, Zhang2017a, Zhuang2017} or knowledge graphs \cite{Chen2018}. Other works have focused on building robust visual features \cite{Goutsu2018, Newell2017, Yang2018b}. However, these methods do not fully explore the potential of feature-level message sharing in order to achieve visual relationship recognition.
Motivated by recent advances in modeling neural message passing \cite{Battaglia2018, Gilmer2017, Santoro2017, Raposo2017, Online2009, Kipf2018}, we aim to extract the relational component of scene graphs by inducing a relational inductive bias in the architecture. From isolated object detection as well as visual and spatial features, we refine the object detection and infer the predicates between them. Graph Networks~\cite{Battaglia2018} introduced the idea of a global attribute that encodes shared information across the graph. Inspired by this work, we propose to use as a global attribute, an image-based context agent, that encodes image-level properties. The context aggregates information from all the nodes and edges through an attentional mechanism to retain only the relevant information.

Previous works have shown interest in relational networks for scene graph generation \cite{Xu2017, Yang2018a, Dai2017, Li2017, Li2017a, Li2018}. For instance, Xu et al.~\cite{Xu2017} proposed an iterative algorithm to pass messages between the nodes and the edges through a primal-dual graph. However, in their approach, a node can pass messages only to its neighboring edges.
Messages between objects must pass through an edge connecting the two objects. Similarly, messages between edges must pass through objects. The information flow is therefore constrained by these node-to-edge and edge-to-node messages. In contrast, our image-level context agent allows fast information propagation across the graph.
Recently, Li et al.~\cite{Li2018} proposed an efficient method to tackle the problem of the quadratic combinations of possible relationships with a relational proposal network. Interactions between objects, predicates and subjects are then modeled with an attention mechanism.
The importance of global context was shown by Zellers et al.~\cite{Zellers2017}, who transformed the graph topology of the problem into a sequence by stacking objects and passing them through a bi-LSTM that encodes object-based contexts. First, the context is used to infer the objects and then the predicates. The context is therefore object-centric. Unlike the previous works that benefit from object-centric context, we opt for an efficient image-level attentional context that shares information across the entire graph. Our contributions are:
\begin{enumerate}

    \item  We demonstrate the importance of image-level context for object and predicate classification.

    \item We propose an iterative graph-propagation network to encode different semantic levels of scene understanding. The network jointly refines the object detection process by modeling its context while predicting relationships between object pairs.

    \item We propose an efficient approach that uses low-dimensional
    edge/node representations and an image-level context. Our method achieves comparable performance with the state-of-the-art, while using less parameters.
\end{enumerate}

\section{The impact of context} \label{sec:context}

Previous works have shown that objects are conditionally dependent on their context \cite{Chen2017, Marino2017, Chen2018}. For instance, knowing that a \textit{tree} is in the scene raises our expectation about seeing another \textit{tree}. Behind this intuition, we can draw a homogeneity principle stating that the objects appearing in a scene tend to be compatible with their environment. To validate this idea, we studied the information gain when we partially reveal the context of an object. We computed the empirical object distribution $\Tilde{p}(o)$, derived with frequency counts from the Visual Genome training data (Xu et al. \cite{Xu2017} cleansed version). We then computed the object distribution conditioned by one randomly chosen object from the same image $I$, denoted by $\Tilde{p}(o|c_i)$, where the context $c=\{c_1, c_i,.., c_n\}$ represents the objects in $I$. Due to the limited dataset size, obtaining a meaningful empirical distribution for $\Tilde{p}(o|c)$ is not feasible. We assume that the lessons learned from studying $\Tilde{p}(o|c=c_i)$ can be generalized to $\Tilde{p}(o|c)$. To determine whether this reasoning applies to relationships, we computed the predicate distribution with and without conditioning. From the empirical distributions, we computed the entropy, conditional entropy, and the maximum entropy for the objects and the predicates (see Table~\ref{tab:entropy}).

\begin{table}[h]
  \caption{Empirical entropy study of the Visual Genome dataset.}
  \label{tab:entropy}
  \centering
  \begin{tabular}{lccc}
  \toprule
    Category     &   $H_{max}$ & $H(x)$ & $H(x|c_i)$ \\
    \midrule
    Objects &   $5.01$  & $4.61$ & $4.08$     \\
    Predicates & $3.91$ & $2.44$ & $2.37$     \\
    \bottomrule
  \end{tabular}
\end{table}

By exposing an object to its context, we observe a significant decrease of its entropy, implying that conditioning can help object detection. To a lesser extent, the same behavior is observed for predicates.
Predicate classification is also a strongly imbalanced class problem with a long tail distribution as shown by the large difference between the maximum entropy and the actual entropy.
Previous works have already shown that objects are highly predictive of the predicates \cite{Zellers2017, Xu2017}. These observations highlight the relational nature of scene graphs and enforce the idea of building contextualized models.

\section{Methodology} \label{sec:approach}

ILAC takes as input an image with the bounding boxes of the objects, and outputs a scene graph. We build a graph $G=\{V,E,C\}$ represented by nodes, edges and a context. The nodes represent the objects, the edges represent the predicates and the context encodes image-level properties (e.g.~whether the scene captures a street or a room).
We propose a composable directed graph-to-graph approach, where we iteratively refine the edge/node/context representations. The pipeline overview is shown in Figure~\ref{fig:pipeline}.

\begin{figure}[h]
\centering
\def\columnwidth{\columnwidth}
\input{scene_graph_pipeline_v3.pdf_tex}
\caption{Pipeline overview.
In (a), we start with bounding boxes (nodes) and union boxes (edges) associated with appearance features. In (b), we compute the image-level attentional context model (ILAC) from the node/edge aggregation (denoted by Agg). The context is used to refine each node and edge through the node and edge LSTM. In (c), we recursively apply this process until convergence. In (d), after $N$ steps, the object/predicate labels are derived with a classification layer. A confidence score is associated to each predicted phrase.
}
\label{fig:pipeline}
\end{figure}
\textbf{Initialization} The input features of the graph are a set of $N_o$ objects.
Each object $o_i$ is represented by its bounding box $b_i$, visual features $f_i \in \mathbb{R}^{4096}$ extracted from a pre-trained net and soft labels predicted by an object detector in isolation \footnote{We provide more details on the object detector in isolation in the supplemental.}, denoted by $\Tilde{p}(o_i) \in \mathbb{R}^{|Cls|}$.
We build a fully connected graph by connecting each object, resulting in $N_o(N_o-1)$ edges. The initial node $v_i \in \mathbb{R}^{d_v}$ and edge $e_{ij} \in \mathbb{R}^{d_e}$ representations are constructed as:
\begin{align}
    v_i^{(0)} = W_o [\Tilde{p}(o_i); b_i; f_i] \quad \mbox{and} \quad  e_{ij}^{(0)} = W_r f_{ij}
\end{align}
where $[.;.]$ denotes the concatenation operation. $f_i$ and $f_{ij}$ are visual features from object $i$ and from the union box of objects $i$ and $j$ respectively. $W_o \in \mathbb{R}^{(d_v \times (|Cls|+|b_i|+|f_i|)}$ and $W_r \in \mathbb{R}^{(d_e \times (|f_{ij}|)}$  are learned embeddings.
The image context $c^{(n)} \in \mathbb{R}^{d_c}$ is initialized with zeros.

\textbf{Recursion}
We first infer the image context from the aggregation of the nodes and the edges using an attention mechanism \cite{Bahdanau2015, Velickovic2017}. The intuition behind is that some objects from the scene are more explanatory of the overall context than others. By computing a weighted sum of the objects, we allow the network to learn discriminative objects. For instance, to infer that we are in a \textit{street}, knowing that a \textit{traffic light} is in the scene is more explanatory than the presence of a \textit{man}. A similar reasoning applies to relations. The context gathers information from the objects $c_o^{(n)}$ and from the relations $c_{r}^{(n)}$. We model the object context $c_o^{(n)}$ as:
\begin{align} \label{eq:context_o}
    c_o^{(n+1)} = \sum_i^{N_o} \alpha_i^{(n)} v_i^{(n)} \quad \mbox{with} \quad \alpha_i^{(n)} = \frac{\exp(w_o^T \sigma(\Phi_o [v_i^{(n)} ; c^{(n)}]))}{\sum_j^{N_o} \exp(w_o^T \sigma(\Phi_o [v_j^{(n)} ; c^{(n)}]))}
\end{align}
where $\sigma$ is the $\tanh$ activation function, $\Phi_o \in \mathbb{R}^{d_{\Phi} \times (d_v + d_c)}$ and $w_o \in \mathbb{R}^{d_\Phi}$ are learned parameters.
A similar attentional approach is used to model the relation context $c_r^{(n)}$.
The node and edge context are then used to update the image context $c^{(n)}$ as:
\begin{align} \label{eq:context}
    c^{(n+1)} = f_c([c_o^{(n+1)} ; c_{r}^{(n+1)}], c^{(n)})
\end{align}
The edges are then updated as follows:
\begin{align}
    e_{ij}^{(n+1)} = f_e([v_i^{(n)} ; v_j^{(n)} ; c^{(n+1)}] , e_{ij}^{(n)})
\end{align}
where $c^{(n+1)}$ is the previously computed context update, $v_i^{(n)}$, $ v_j^{(n)}$ are the node of object $i$ and $j$, respectively.
We then update the nodes as:
\begin{align}
    v_i^{(n+1)} = f_v([e_{i\bullet}^{(n)} ; c^{(n)}], v_i^{(n)})
\end{align}
where $e_{i\bullet}$ is the aggregation of all the incoming edges at node $i$ computed as in \eqref{eq:context_o}.
Following previous works \cite{Palm2017, Sanchez-Gonzalez2018, Xu2017}, $f_c$, $f_e$ and $f_v$ are RNNs, namely LSTMs, where the cell hidden state represents the current edge/node/context. The hidden state is therefore iteratively refined as information from the graph diffuses. After $n$ iterations, the node and edge hidden states are passed through a classification layer to derive the labels.
The hyperparameters and network dimensions are listed in Table \ref{tab:hyperparameters}.
\begin{table}[h]
  \caption{Network dimensions and hyperparameters used.}
  \label{tab:hyperparameters}
  \centering
  \begin{tabular}{cccccccc}
  \toprule
    Name & $d_c$ & $d_v$ & $d_e$ & $d_{\Phi}$ & learning rate & batch size & iterations \\
    \midrule
    Value &   $512$  & $512$ & $512$ & $256$ & $10^{-4}$ & $8$  & $2$   \\
    \bottomrule
  \end{tabular}
\end{table}
\section{Experimental results} \label{sec:results}
\vspace{-1mm}

We used the Visual Genome dataset \cite{Krishna2017} cleansed by \cite{Xu2017} for training. It contains $150$ object categories and $50$ different predicates. During training, we use as loss function the sum of the cross-entropy loss of the objects and predicates. We evaluate the top-K recall (R@K) for the \textbf{predicate classification} (PredCls) and \textbf{scene graph classification} (SGCls). Note that we focus on the relational part of scene graphs and therefore we assume that the object location is known.
We also report the contextualized object-detection accuracy
\footnote{We provide more details on the evaluation metrics in the supplemental section.}.

Experimental results are presented in Table \ref{table:result}. We compare our approach, ILAC, with IMP \cite{Xu2017}, neural motifs (NM) \cite{Zellers2017}, and a frequency baseline derived from statistics and object detection in isolation (FREQ) \cite{Zellers2017}.
ILAC is found to be superior to IMP and FREQ, and have comparable performance with NM. However, ILAC uses $\approx 3 \times$ fewer parameters than NM,
as can be seen in the last column of the table. This is achieved by keeping a unique global context and embedding the visual features into a lower-dimensional space.
More importantly, we improve the object detection in isolation (FREQ) by $1.5\%$, thus exhibiting the power of ILAC to extract the relational nature of scenes.
We also performed an ablation study to understand the contribution of each component. We first bypassed the image-level context module, bounding the model to node-to-edge and edge-to-node messages only (No-context in Table \ref{table:result}). The importance of the image-level context module is proven by experimental results as we observe a $3\%$ absolute performance gain for SGCls and $\approx 1\%$ for PredCls. This effect is more pronounced for SGCls than PredCls, as SGCls requires more relational reasoning and can therefore fully benefit from ILAC. Note that this observation matches the result of the entropy study. Beside having two iterations for ILAC, we also trained it with a single-pass graph network (1-iter). Having two iterations leads to an improvement of $\approx 0.6\%$,
which shows the gain by iteratively refining the node/edge/context representations. It should be noted that we observed no further improvements upon adding additional iterations.
\vspace{-2mm}
\begin{table}[h]
  \caption{Results expressed in $\%$. Number of parameters expressed in millions. All results were generated using the same object detector (from~\cite{Zellers2017}), evaluation service and test/train splitting to ensure fair comparisons. ILAC results have settings as listed in Table \ref{tab:hyperparameters}.}
  \label{table:result}
  \centering

  \begin{tabular}{lcccccccccc}
    \toprule
    \multicolumn{1}{c}{} &
    \multicolumn{2}{c}{PredCls} &
    \multicolumn{2}{c}{SGCls}  &
    \multicolumn{1}{c}{Object accuracy} &
    \multicolumn{1}{c}{Parameters}                   \\
    \cmidrule(r){2-3}
    \cmidrule(r){4-5}
    Model & R$@50$ & R$@100$ & R$@50$ & R$@100$     \\
    \midrule
    FREQ \cite{Zellers2017} & $75.2$ & $83.6$ & $39.0$ & $43.4$ & $66.1$ &  - \\
    IMP \cite{Xu2017} & $75.7$ & $82.9$  & $43.4$ & $47.2$ & $-$& $10.0$ \\
    NM \cite{Zellers2017} & $81.1$ & $88.3$ & $44.5$ & $47.7$ & $67.3$ & $40.0$ \\
    \textbf{ILAC} &  $\mathbf{80.4}$ & $\mathbf{88.4}$ & $\mathbf{43.7}$ & $\mathbf{47.6}$ & $\mathbf{67.6}$ & $\mathbf{14.8}$ \\
    \midrule
    No-context &  $78.8$ & $87.6$ & $40.3$ & $44.6$ & $66.2$ & $14.8$  \\
    ILAC (1-iter) &  $79.8$ & $87.9$ & $43.0$ & $47.0$  & $67.4$ & $14.8$\\
    \bottomrule
  \end{tabular
\end{table}

\vspace{-4mm}
\section{Conclusion} \label{sec:conclusion}
\vspace{-1mm}
This paper introduces an efficient method for achieving scene graph generation, called the image-level attentional context (ILAC) modeling. The ILAC model learns to refine graph nodes and edges iteratively, then propagates contextual information across the graph. We demonstrate that sharing the information on a graph through a context agent allows to create competitive models that have less parameters than the state-of-the-art.
More importantly, we show that image-level contexts can refine object detection by providing relational information.
Moreover, ILAC does not require domain-specific or external knowledge, which makes it easily pluggable on top of existing object detectors (with or without predicate information) to account for the relational nature of objects in images.

\section*{Acknowledgements}

The authors would like to thank their IBM colleagues Anca-Nicoleta Ciubotaru, Andrea Giovannini, Antonio Foncubierta-Rodriguez and Apostolos Krystallidis for their feedback on this manuscript.

\medskip
\small
\bibliographystyle{abbrv}