# Learning Modulo Theories

Matt Fredrikson  Kaiji Lu  Saranya Vijayakumar  Somesh Jha  Vijay Ganesh
Zifan Wang

###### Abstract

Recent techniques that integrate _solver layers_ into Deep Neural Networks
(DNNs) have shown promise in bridging a long-standing gap between inductive
learning and symbolic reasoning techniques. In this paper we present a set of
techniques for integrating _Satisfiability Modulo Theories_ (SMT) solvers into
the forward and backward passes of a deep network layer, called
$\mathsf{SMTLayer}$. Using this approach, one can encode rich domain knowledge
into the network in the form of mathematical formulas. In the forward pass,
the solver uses symbols produced by prior layers, along with these formulas,
to construct inferences; in the backward pass, the solver informs updates to
the network, driving it towards representations that are compatible with the
solver’s theory. Notably, the solver need not be differentiable. We implement
$\mathsf{SMTLayer}$ as a Pytorch module, and our empirical results show that
it leads to models that _1)_ require fewer training samples than conventional
models, _2)_ that are robust to certain types of covariate shift, and _3)_
that ultimately learn representations that are consistent with symbolic
knowledge, and thus naturally interpretable.

Machine Learning, ICML

##  1 Introduction

A recent class of techniques aims at integrating solver layer(s) within deep
neural networks (DNNs) (Wang et al., 2019a; Pogancic et al., 2020; Huang et
al., 2021; Manhaeve et al., 2018), both during training and inference. A class
of problems which can benefit from such an integration is one that has both a
perceptual and a symbolic sub-problem, such as “visual” Sudoku (Wang et al.,
2019a), or the problem of determining the shortest path from a picture of a
map (Pogancic et al., 2020).

The most straightforward way to incorporate a solver layer into an ML model is
to learn models with representations that are compatible with symbols used by
the solver. For example, if one wanted to leverage symbolic domain knowledge
to classify images of birds, or diagnose ailments from CT scans, then one
could train a model in a fashion similar to “concept bottlenecking” (Koh et
al., 2020). This requires detailed labels for supervision, which may be
prohibitively expensive to obtain and keep consistent with a potentially
evolving domain theory.

We present a set of techniques for incorporating a _Satisfiability Modulo
Theories_ (SMT) solver into a DNN layer so that symbolic knowledge can be
leveraged to learn such a compatible representation, _without requiring label
supervision_. Our approach is general, and can handle a broad range of domain
knowledge encoded as SMT constraints, provided that they interface with the
surrounding neural network layers over propositional variables. Unlike the
most closely related prior work (Wang et al., 2019a), our approach does not
approximate the solver’s behavior by formulating a differentiable relaxation.
Rather, we extract information from the solver as it works on a set of
constraints, that is geared towards checking the correctness of the output of
the model that precedes the solver, and use that information to construct
updates to the model during training (Section 4.2).

We present two different approaches for this, one based on unsatisfiable
cores, and another based on weighted MaxSMT (Section 4.2). There are several
advantages to this approach. Aside from the mild interface constraints
mentioned earlier (i.e., solver and neural layers interface with each other
via boolean variables), our approach does not place any restrictions on the
theory solver embedded in the layer, such as linearity (Pogancic et al., 2020)
or even decidability—if the solver is capable of efficiently discharging the
relevant constraints, then the layer can operate as intended. Because there is
no need to provide a differentiable relaxation for each theory or solver
technique that one may want to incorporate, we can leverage the continuous and
unabated progress being made in solver technology.

We implement our approach as a PyTorch (Paszke et al., 2019) layer, using the
Z3 (De Moura & Bjørner, 2008) SMT solver as the solver layer to solve SMT and
MaxSMT constraints. On three applications involving vision and natural
language: visual arithmetic, algebraic equation solving, and a so-called
natural language “liar’s puzzle,” we demonstrate that our implementation can
be incorporated into DNN architectures to solve problems more effectively than
conventional DNNs (Section 5). In particular, our results show that the data
needed to train a DNN with symbolic knowledge may be much simpler than may be
necessary otherwise, and that while doing so is more expensive
computationally, often times the more efficient (i.e., not involving MaxSMT)
algorithms perform well in practice.

Our contributions are as follows:

  1. 1.

We present $\mathsf{SMTLayer}$, a framework for incorporating an SMT solver
into a DNN, as a layer that leverages symbolic knowledge during training and
inference.

  2. 2.

We prototype our approach in Pytorch111We plan to release our implementation
as an open source library upon publication of this paper, and show that it can
be applied to solve a range of problems that incorporate symbolic knowledge.

  3. 3.

Our empirical evaluation, over four diverse applications, shows that models
using $\mathsf{SMTLayer}$ require significantly less training data, can be
trained more efficiently, and are more robust than those based on closely-
related prior work (Wang et al., 2019a; Huang et al., 2021).

Section 3 provides background on ERM and the first-order theories used in our
framework. Section 4 describes $\mathsf{SMTLayer}$, Section 5 gives our
empirical evaluation, and Section 6 concludes the paper.

##  2 Related Work

Combining logical solvers and deep models can be difficult because logic has
discrete structure while the most successful way to construct neural networks
today requires differentiability (Riegel et al., 2020).

#### Combinatorial Solver Layers.

Vlastelica et al. (2019) integrate a blackbox, non-differentiable
combinatorial solver on top of a deep network. To propagate the gradient
through the solver on the backward pass, they linearly interpolate the loss
w.r.t the solver’s input and define the gradient of the solver as the slopes
of the line segments. CSL solves a set of problems where the solver’s
objective must be linear w.r.t its input, e.g. finding the shortest path and
travel salesman problem (TSP). Further, the authors assume that the only
labels available are the outputs of solvers, e.g. the minimum cost in TSP, and
hence their tool has to discover the label for the output of the network
itself. These requirements limit the choices one has for the solver layer.

#### Neural Logic Programming.

While SATNet integrates a logic-based solver on top of a network, DeepProbLog
takes the opposite approach, extending the capability of a probabilistic logic
programming language with neural predicates Manhaeve et al. (2018). In the
context of our work, the logic program can be viewed as a “solver layer” that
explicitly encodes symbolic knowledge. Scallop (Huang et al., 2021) extends
DeepProbLog to scale without sacrificing accuracy compared to DeepProbLog.
Similarly to DeepProbLog, each possible result of the sum of two digits in
MNIST is given a probability, in the form of a weighted Boolean formula. They
prune unlikely clauses of the formula, represented by proofs, only keeping the
top-$k$ most likely. Likelihood is computed using weighted model counting
(Huang et al., 2021; Chavira & Darwiche, 2008). These techniques are well-
suited to problems that benefit from probabilistic Datalog, but have inherent
limitations: they cannot handle quantifiers, general negation, and the range
of supported first-order theories is more restrictive.

#### SATNet.

Wang et al. (2019b) present SATNet, a network architecture with a
differentiable approximate MAXSAT solver layer. Their approximation is based
on a coordinate descent approach to solving the semidefinite program (SDP)
relaxation of the MAXSAT problem. SATNet does not assume that the logical
structure of the problem is given, and instead attempts to learn it. By
placing the MAXSAT solver layer on top of a convolution network to learn
representations from images, SATNet directly solve problems like Visual
Sudoku, for which neural networks alone are not well suited (Wang et al.,
2019b).

#### Differentiable Logic.

Another recent direction has explored differentiable logics (Fischer et al.,
2019; Varnai & Dimarogonas, 2020; van Krieken et al., 2022). These approaches
provide ways of integrating symbolic knowledge into training, by making
logical formulas differentiable, and therefore amenable to optimization when
included in a loss function. This line of work does not explicitly aim to make
use of symbolic information during inference. In contrast, the information
that our approach extracts from the solver during training is used to
condition the model towards a representation that will allow it to communicate
effectively with the solver during inference. Additionally, we do not require
the logical formulas, or the solver, to be differentiable.

##  3 Background

Let $\mathbf{X}$ denote a domain of features, $\mathbf{Y}$ a domain of labels,
and $\mathcal{D}$ a distribution over $\mathbf{X}\times\mathbf{Y}$. Formally,
$\mathcal{D}$ is a probability measure on a space given by a $\sigma$-algebra
over subsets of $\wp(\mathbf{X}\times\mathbf{Y})$. The goal of a learning
algorithm $A$ is to find a function $h:\mathbf{X}\to\mathbf{Y}$ that, for
$(x,y)\sim\mathcal{D}$, can be used to predict $y$ when given $x$. To do this,
$A$ is given a set of training examples $S=(x_{1},y_{1}),\ldots,(x_{m},y_{m})$
sampled i.i.d. from $\mathcal{D}$, and uses some criterion to select $h$ from
a hypothesis class $\mathcal{H}$ of functions. We refer to $h$ as the _model_
learned by $A$ on $S$. When the learning algorithm $A$ is clear from the
context, we will write $h_{S}$ to denote the model produced from the given
sample. Throughout this paper, we will generally assume that the loss is
either the 0-1 loss $\ell^{01}$ or binary cross-entropy $\ell^{\mathrm{bce}}$.

A theory $T$ consists of a signature $\Sigma$ of constant, predicate, and
function symbols, as well as a set of axioms over $\Sigma$. Formulas in a
theory are composed of elements of $\Sigma$, variables, and logical symbols
such as quantifiers and Boolean operations. We use the term _decision
procedure_ to refer to an algorithm that is given an open $T$-formula, and
returns _true_ if it is satisfiable, and _false_ otherwise. Additionally, it
may return an assignment to all of the variables that demonstrates
satisfiability, or if the formula is not satisfiable, then it may return an
_unsatisfiable core_ , which is a subset of clauses taken from the formula’s
representation in conjunctive normal form that remains unsatisfiable. Loosely,
we also refer to such an algorithm as a “solver”, but this term is more
general, and could also refer to an algorithm that identifies the maximal set
of clauses, possibly weighted by some user-defined values, that are
satisfiable when conjoined.

##  4 Constructing $\mathsf{SMTLayer}$

In this section, we present $\mathsf{SMTLayer}$, a set of algorithms for
computing the forward and backward passes of a layer whose behavior is defined
by a set of user-defined SMT constraints. $\mathsf{SMTLayer}$ does not have
trainable parameters, and its functionality is wholly defined by a set of SMT
constraints $\phi$ that are provided by the model designer.
$\mathsf{SMTLayer}$ can be used in modern deep-learning frameworks as a drop-
in replacement for more conventional neural network layers, e.g., dense,
convolutional, and LSTM (Hochreiter & Schmidhuber, 1997) are prominent
examples of widely-used layers.

Section 4.1 provides a high-level overview of our approach, Section 4.2
describes them in detail, and Section 4.3 begins an analysis of this setting
that we hope future work will continue developing.

Features $\mathbf{X}$ Symbolic Domain $\mathbf{Z}$ $\mathsf{0010000111}$
Neural Network $\begin{array}[]{ll}\phi(\quad z_{1}\|\ldots\|z_{10},\quad
y\quad)\ \equiv&\\\ a=\mathds{1}_{z_{1}>0}\|\ldots\|\mathds{1}_{z_{5}>0}\land
b=\mathds{1}_{z_{5}>0}\|\ldots\|\mathds{1}_{z_{10}>0}&\land\\\
a+b=y&\end{array}$ Prediction Logic $\phi$ Labels $\mathbf{Y}$
$\\{\mathsf{01011}\\}$ Satisfying Assignments Figure 1:  MNIST Addition
example.

###  4.1 Overview

We envision $\mathsf{SMTLayer}$ being used primarily at the top of a DNN
taking inputs from a stack of conventional DNN layers that convert raw input
features into ground terms for the constraints
$\phi(z_{0},\ldots,z_{p-1},y_{0},\ldots,y_{q-1})$ embedded in
$\mathsf{SMTLayer}$, and producing outputs that are consistent with $\phi$ and
the given ground terms. Figure 1 shows an illustrative example, with the
previously-studied problem of MNIST addition (Manhaeve et al., 2018; Huang et
al., 2021).

During the forward pass the outputs of the previous layer are mapped to
designated free variables $z_{0},\ldots,z_{p-1}$. The layer then checks the
satisfiability of $\phi$, a formula in an appropriate combination of first-
order theories, after substituting these ground terms for the $z_{i}$, and the
output of the layer consists of the solver’s model for $y_{0},\ldots,y_{q-1}$.
These outputs are converted from Boolean to floating-point values by mapping
_false_ to -1 and _true_ to 1. At the moment, the only restriction on $\phi$
that our layer requires is that $z$ and $y$ be vectors of Booleans, so that
they can be appropriately mapped to continuous values; any other symbols
appearing in $\phi$ can come from arbitrary domains (e.g. strings) supported
by the underlying SMT solver.

In the backward pass, the layer receives the gradient of its output with
respect to the function whose derivative is being computed, which we will
assume is the binary cross-entropy loss $\ell(y,y^{\star})$. Unless stated
otherwise, we will assume this loss for the remainder of the section. This
gradient is used, along with the inputs and outputs of the corresponding
forward pass, to first compute an amended output $\hat{y}$ which corresponds
to an output that would have yielded a smaller loss. Because the outputs are
Boolean, it is always possible to determine the ground truth output
$y^{\star}$ from this information. Using $\hat{y}$, the layer determines which
of components of its inputs are inconsistent with $\phi$ and $\hat{y}$, and
provides the corresponding gradients to the previous layer. Section 4.2
details the manner in which these gradients are computed.

###  4.2 $\mathsf{SMTLayer}$, forward and backward

Inputs: $z\in\mathbb{R}^{p}$ layer input
$\phi(z_{0},\ldots,z_{p-1},y_{0},\ldots,y_{q-1})$ $T$-formula

Output: $y\in\mathbb{R}^{q}$

begin

2 $z_{b}\leftarrow[z[i]>0~{}:~{}i=0\ldots p-1]$ $C\leftarrow\sum_{i\in
I}\mathsf{softmax}(|z|)[i]$
$y\leftarrow\operatorname*{arg\,max}_{y_{b}}\max_{I}C\cdot\mathds{1}(\phi\land\bigwedge_{i\in
I}y_{i}=y_{b}[i]\land z_{i}=z_{b}[i])$ return $y$

3 end

Algorithm 1 $\mathsf{F}_{\mathit{max}}^{\phi}(z)$
MaxSMT-based forward pass of $\mathsf{SMTLayer}$

Inputs: $z\in\mathbb{R}^{p}$ layer input
$\phi(z_{0},\ldots,z_{p-1},y_{0},\ldots,y_{q-1})$ $T$-formula

Output: $y\in\mathbb{R}^{q}$

begin

2 $z_{b}\leftarrow[z[i]>0~{}:~{}i=0\ldots p-1]$
$\hat{\phi}\leftarrow\phi(z_{b}[0],\ldots,z_{b}[p-1])$ if  _$\hat{\phi}$ is
satisfiable_ then

3
$y_{b}[0],\ldots,y_{b}[q-1]\leftarrow\mathsf{solve}(\hat{\phi},y_{0},\ldots,y_{q-1})$
$y\leftarrow[y_{b}[i]>0\quad~{}:~{}i=0\ldots q-1]$

4 else

5 $y\leftarrow\mathbf{0}$

6 end if

7 return $y$

8 end

Algorithm 2 $\mathsf{F}_{\mathit{smt}}^{\phi}(z)$
SMT-based forward pass of $\mathsf{SMTLayer}$

We now present the details of the forward and backward passes of
$\mathsf{SMTLayer}$. There are two algorithms for each pass,
$\mathsf{F}_{\mathit{max}}^{\phi}$ and $\mathsf{F}_{\mathit{smt}}^{\phi}$ are
forward passes, and $\mathsf{B}_{\mathit{max}}^{\phi}$,
$\mathsf{B}_{\mathit{core}}^{\phi}$ are backward passes.
$\mathsf{F}_{\mathit{max}}^{\phi}$ and $\mathsf{B}_{\mathit{max}}^{\phi}$ both
make use of MaxSMT solvers, whereas $\mathsf{F}_{\mathit{smt}}^{\phi}$ and
$\mathsf{B}_{\mathit{core}}^{\phi}$ rely on satisfiability solvers (SMT).
Despite the symmetry in which type of solver each algorithm uses, they are all
compatible with each other. That is, $\mathsf{F}_{\mathit{smt}}^{\phi}$ can be
used with either $\mathsf{B}_{\mathit{max}}^{\phi}$ or
$\mathsf{B}_{\mathit{core}}^{\phi}$, and the same for
$\mathsf{F}_{\mathit{max}}^{\phi}$.

#### Forward pass.

Algorithms 1 and 2 illustrate $\mathsf{F}_{\mathit{max}}^{\phi}$ and
$\mathsf{F}_{\mathit{smt}}^{\phi}$, the methods for computing the forward pass
based on weighted MaxSMT and SMT, respectively. Both of the algorithms are
parameterized by a user-provided first-order formula $\phi$, and take a single
vector-valued input consisting of unscaled floating-point values (_logits_).
These values are cast to Boolean constants by taking their sign on line 2 of
both algorithms, so that they can be equated with the corresponding free
variables $z_{0},\ldots,z_{p-1}$.

The key difference between $\mathsf{F}_{\mathit{max}}^{\phi}$ and
$\mathsf{F}_{\mathit{smt}}^{\phi}$ is the way in which they handle inputs that
are inconsistent with $\phi$ when interpreted as Booleans.
$\mathsf{F}_{\mathit{smt}}^{\phi}$ addresses this by providing an output that
is also inconsistent with $\phi$, i.e. a vector of zeroes, effectively
signaling that the network below it did not provide consistent inputs.
Alternatively, we can interpret the values provided by the network as Booleans
enriched with “confidence” values. Although we expect inputs to
$\mathsf{SMTLayer}$ to be unscaled floating-point values, Algorithm 1 scales
them to a formal probability distribution via the softmax function (line 3)
for use as weights to find the weighted MaxSMT solution of $\phi$. With this
approach, $\mathsf{SMTLayer}$ will always provide a valid, although not
necessarily correct, output that is consistent wrt $\phi$ with the inputs of
which the network below is most “confident.” (line 4).

Inputs: $z\in\mathbb{R}^{p}$ input of forward pass
$y\in\mathbb{R}^{q}$ output of forward pass
$\partial_{y}\ell(y,y^{\star})$ gradient with respect to output
$\phi(z_{0},\ldots,z_{p-1},y_{0},\ldots,y_{q-1})$ a $T$-formula

Output: $\partial_{z}\ell(y,y^{\star})\in\mathbb{R}^{p}$ approximate gradient
of $\ell$

begin

2 $G_{z}\leftarrow\partial_{z}\ell(z,\mathsf{sign}(z))$
$\hat{y}=\mathsf{sign}(y)-2\cdot\mathsf{sign}\left(\partial_{y}\ell(y,y^{\star})\right)$
if  _$\mathsf{sign}(y)\neq\mathsf{sign}(\hat{y})$_ then

3 $z_{b}\leftarrow[z[i]>0~{}:~{}i=0\ldots p-1]$
$\hat{y}_{b}\leftarrow[\hat{y}[i]>0~{}:~{}i=0\ldots q-1]$
$\phi_{y}\leftarrow\bigwedge_{0\leq i<q}y_{i}=\hat{y}_{b}[i]$
$C\leftarrow\sum_{i\in I}\mathsf{softmax}(|z|)[i]$
$I\leftarrow\operatorname*{arg\,max}_{I\subseteq[0,p)}\mathds{1}(\phi\land\phi_{y}\land\bigwedge_{i\in
I}z_{i}=z_{b}[i])\cdot C$ foreach  _$i\in\bar{I}$_ do

4 $G_{z}[i]\leftarrow\partial_{z[i]}\ell(z[i],1-\mathsf{sign}(z[i]))$

5 end foreach

6

7 end if

8 return $G_{z}$

9 end

Algorithm 3
$\mathsf{B}_{\mathit{max}}^{\phi}(z,y,\partial_{y}\ell(y,y^{\star}))$
MaxSMT-based backward pass of $\mathsf{SMTLayer}$

Inputs: $z\in\mathbb{R}^{p}$ input of forward pass
$y\in\mathbb{R}^{q}$ output of forward pass
$\partial_{y}\ell(y,y^{\star})$ gradient with respect to output
$\phi(z_{0},\ldots,z_{p-1},y_{0},\ldots,y_{q-1})$ a $T$-formula

Output: $\partial_{z}\ell(y,y^{\star})\in\mathbb{R}^{p}$ approximate gradient
of $\ell$

begin

2 $G_{z}\leftarrow\partial_{z}\ell(z,\mathsf{sign}(z))$
$\hat{y}=\mathsf{sign}(y)-2\cdot\mathsf{sign}\left(\partial_{y}\ell(y,y^{\star})\right)$
if  _$\mathsf{sign}(y)\neq\mathsf{sign}(\hat{y})$_ then

3 $z_{b}\leftarrow[z[i]>0~{}:~{}i=0\ldots p-1]$
$\hat{y}_{b}\leftarrow[\hat{y}[i]>0~{}:~{}i=0\ldots q-1]$
$\phi_{z},\phi_{y}\leftarrow\bigwedge_{0\leq
i<p}z_{i}=z_{b}[i],\,\,\,\bigwedge_{0\leq i<q}y_{i}=\hat{y}_{b}[i]$
$I\leftarrow\operatorname*{arg\,min}_{I\subseteq[0,p)}\mathds{1}(\lnot\phi\lor\lnot\phi_{y}\lor\bigvee_{i\in
I}z_{i}\neq z_{b}[i])\cdot|I|$ foreach  _$i\in I$_ do

4 $G_{z}[i]\leftarrow\partial_{z[i]}\ell(z[i],1-\mathsf{sign}(z[i]))$

5 end foreach

6 foreach  _$i\in\bar{I}$_ do

7 $G_{z}[i]\leftarrow 0$

8 end foreach

9

10 end if

11 return $G_{z}$

12 end

Algorithm 4
$\mathsf{B}_{\mathit{core}}^{\phi}(z,y,\partial_{y}\ell(y,y^{\star}))$
unsat core-based backward pass of $\mathsf{SMTLayer}$

#### Backward pass.

The backward pass is responsible for computing the gradient of the loss with
respect to the layer inputs. It is given the gradient of the loss with respect
to the layer outputs, and is assumed to have memoized the inputs that it
received previously on the forward pass, as well as the outputs that they
produced. The gradients returned by this pass are then used by the backward
pass of previous layers, and ultimately to derive updates to trainable
parameters that will yield smaller loss.

The key issue in designing a backward pass for $\mathsf{SMTLayer}$ is the
geometry of the functions computed by either forward pass. For any vector
$v\in\\{-1,0,1\\}^{p}$ and $x,x^{\prime}$ with
$\mathsf{sign}(x)=\mathsf{sign}(x^{\prime})=v$, then
$\mathsf{B}_{\cdot}^{\phi}(x)=\mathsf{B}_{\cdot}^{\phi}(x^{\prime})$, so these
functions are piece-wise constant step functions ranging over the corners of
the $\mathbb{R}^{q}$ unit hypercube. Thus, while they are differentiable
almost everywhere, the gradient is not helpful for training because it is
always zero. Prior work on integrating such step functions into deep networks
primarily addresses this problem by relaxing the function computed by the
forward pass, so that its gradients are no longer constant and hopefully more
informative.

In contrast, $\mathsf{B}_{\mathit{max}}^{\phi}$ (Algorithm 3) and
$\mathsf{B}_{\mathit{core}}^{\phi}$ (Algorithm 4) do not attempt to provide
gradients for a relaxation of the forward pass. Instead, they use information
provided by the solver in its computation of the forward pass to identify
which components of the input may have contributed to higher loss. The
gradient is then computed by constructing a counterfactual variant of the
input provided to the forward pass, which differs on the identified
components, and returning the gradient of the binary cross-entropy loss of the
original input on this counterfactual variant. The two algorithms differ in
the information that they extract from the solver, i.e., either solutions to a
MaxSMT instance or an unsatisfiable core.

Both algorithms begin by initializing the gradient to be the loss between the
logit inputs, and their hard labels (line 2). Recall that we assume the loss
$\ell$ is binary cross-entropy, so the result will not be zero. The purpose of
this initialization is to emulate the dynamics of training with cross-entropy
loss with a conventional layer; when the rounded output matches the target,
the loss is not zero, and training will continue to move the parameters in a
direction that makes them agree “more” with the hard target. One line 3, they
then use the provided gradient from the subsequent layer together with the
memoized output from the forward pass to construct $\hat{y}$, a “corrected”
output that satisfies $\ell(\hat{y},y^{\star})\leq\ell(y,y^{\star})$. If the
sign of $\hat{y}$ is the same as that of $y$, then both algorithms return the
initialized gradient. Otherwise, they extract information from the solver
using the inputs to the forward pass $z$ and $\hat{y}$.

$\mathsf{B}_{\mathit{max}}^{\phi}$ constructs a set of clauses $\phi_{y}$ that
constrain the free $y_{0},\ldots,y_{q-1}$ to take the values of $\hat{y}_{b}$,
the Boolean conversion of $\hat{y}$. It then computes the softmax values of
the absolute incoming logits $|z|$, and uses them to find the maximally-
weighted set of clauses $(\mathsf{softmax}(|z|)[i],z_{i}=z_{b}[i])$ that are
consistent with $\phi\land\phi_{y}$. Intuitively, these are the inputs that
the previous layer is most confident in that can be made consistent with the
corrected label $\hat{y}$ by changing some of the less confident inputs.
$\mathsf{B}_{\mathit{max}}^{\phi}$ then updates the initialized gradient at
each index for which the solution to this instance does not match the sign of
the original input.

$\mathsf{B}_{\mathit{core}}^{\phi}$ also constructs $\phi_{y}$, but instead
identifies a set of constraints $z_{i}=z_{b}[i]$ that are inconsistent with
$\phi\land\phi_{y}$. Note that line 7 specifies a minimal unsatisfiable core,
but this is not necessary. All that is needed is that none of the clauses in
the core be superfluous, i.e., deleting any singleton clause from $I$ will
cause it to be satisfiable. If a superfluous clause remains in the core, then
the gradient returned for the corresponding input will have the incorrect
sign, which may lead to issues with training.
$\mathsf{B}_{\mathit{core}}^{\phi}$ then updates the gradient at each input
identified in the core using the loss of $z[i]$ with respect to
$1-\mathsf{sign}(z[i])$, which will lead to updates in a direction that would
have modified the input such that $i$ was not in the unsat core. The indices
not in the unsat core have their gradients set to zero, as their absence in
the core is not evidence that these inputs were correct or incorrect.

###  4.3 Analysis

To understand the settings where $\mathsf{SMTLayer}$ will provide optimal
results, we introduce a class of _decomposable_ learning problems (Definition
1).

######  Definition 1 (Decomposable problem).

Let $T$ be a first-order theory with constants in $\mathbf{Z}$. An ERM problem
$\mathcal{D}$, $\mathcal{H}$ is _decomposable_ by $T$ if there exists a
function $f:\mathbf{X}\to\mathbf{Z}$, companion hypothesis class
$\mathcal{H}_{f}\subseteq\mathbf{X}\to\mathbf{Z}$, and $T$-formula $\phi$ such
that:

  1. 1.

For any $h\in\mathcal{H}$, there exists $h_{f}\in\mathcal{H}_{f}$ and
$h^{\prime}$ such that $h=h^{\prime}\circ h_{f}$.

  2. 2.

There exists a random function $g:\wp(\mathbf{Y})\to\mathbf{Y}$ such that for
any $n>0$ and $\forall S$ in the support of $\mathcal{D}^{n}$,

$\Pr_{(x,y)\sim\mathcal{D}}[(x,y)\in
S]=\Pr_{(x,\cdot)\sim\mathcal{D}}[(x,g(\langle x\rangle_{f,\phi}))\in S]$

where $\langle x\rangle_{f,\phi}=\\{y:\phi(f(x),y)~{}\text{is satisfied}\\}$.

In (2), $f$ is called the _grounding function_ and $\phi$ is called the
_prediction logic_.

Intuitively, a learning problem defined in terms of a distribution
$\mathcal{D}$ and hypothesis class $\mathcal{H}$ is decomposable if members of
$\mathcal{H}$ can be decomposed into functions that are responsible for
grounding and prediction, and $\mathcal{D}$ can be expressed in terms of a
grounding function and a first-order formula $\phi$. There are a few important
things to note. First, there is no requirement that the grounding function $f$
be a member of $\mathcal{H}_{f}$. While this may be realized at times, we
should not assume that the data is actually generated, or otherwise described
with perfect fidelity, by a function in the class that one learns over. In
fact, we do not assume that $f$ is efficiently computable, as it may
correspond to a natural process, or an aspect of data generation that is not
understood well enough to make such computational claims.

Second, for a given $x$, there may be more than one satisfying assignment for
$y$ to $\phi(f(x),y)$. The function $g$ in (2) accounts for this, requiring
only that when solutions to $\phi(f(x),y)$ are sampled by $g$, the result is
distributed identically to $\mathcal{D}$. This paper will focus on the case
where satisfying assignments for $y$ are unique, as these are more in line
with ”classic” ERM classification problems. We leave exploration of the more
general setting to future work.

We note that if the grounding function is known, can be computed efficiently,
and $\phi$ is efficiently solvable, then the learning problem effectively has
a closed-form solution. Rather, we assume that only $\phi$ and perhaps $g$ are
known, and a sample of $\mathcal{D}$ is given. The remaining challenge is to
identify a grounding hypothesis $h_{f}\in\mathcal{H}_{f}$ for which the
construction in (2) is an effective solution to the end-to-end learning
problem posed by $\mathcal{D}$, $\mathcal{H}$. This stands in contrast to
traditional ERM, in which a good solution $h\in\mathcal{H}$ must either solve
both grounding and prediction, or find a “shortcut” that manages to predict
$\mathcal{D}$ as well as the decomposition.

#### Convergence.

Regarding the backward passes, Theorem 2 below demonstrates that when $\phi$
satisfies certain conditions, and the companion hypothesis class
$\mathcal{H}_{f}$ satisfies conditions that are sufficient to guarantee
convergence with SGD, then training with $\mathsf{F}_{\mathit{smt}}^{\phi}$
and $\mathsf{B}_{\mathit{max}}^{\phi}$ will converge to the optimal solution
in the number of iterations. The proof of this theorem is based on the
observation that when the conditions on $\phi$ are met, then training with
$\mathsf{B}_{\mathit{max}}^{\phi}$ obtains the same solution that would be
obtained if the labels of $\phi$ were available for supervised learning. Thus,
the conditions on $\mathcal{H}_{f}$ are sufficient to ensure the stated
convergence, as stated in a well-known result outlined in Chapter 14 of
(Shalev-Shwartz & Ben-David, 2014).

It is also worth noting that Theorem 2 does not necessarily hold if
$\mathsf{B}_{\mathit{core}}^{\phi}$ is used instead of
$\mathsf{B}_{\mathit{max}}^{\phi}$. The reason is that there may be many
unsatisfiable cores that are locally minimal in cardinality, and gradients are
set only for inputs that appear in the computed core. These gradients will not
match those of the loss on a grounding sample, so the training dynamics are
likely to be different. We believe that training with
$\mathsf{B}_{\mathit{core}}^{\phi}$ may have more in common with block
coordinate descent than gradient descent, and save a more detailed exploration
of the topic for future work.

######  Theorem 2.

Let $\mathcal{D}$, $\mathcal{H}$ be a $T$-decomposable problem with grounding
function $f$ and prediction logic $\phi$ where:

  1. 1.

$\mathbf{Z}$ and $\mathbf{Y}$ are Cartesian products of Booleans.

  2. 2.

For any $(x,y)\sim\mathcal{D}$ and $y^{\prime}\neq y$, $\phi(f(x),y^{\prime})$
is $T$-equivalent to false and there is exactly one $z$ such that $\phi(z,y)$
is $T$-equivalent to true.

  3. 3.

$\mathcal{H}_{f}$ is a convex set and for all $h_{f}\in\mathcal{H}_{f}$,
$\|h_{f}\|\leq B$, and the loss $\ell(h_{f}(\cdot),z)$ is $M$-Lipschitz and
convex in $x$ for any fixed $z$.

Then for any $\epsilon>0$, selecting $h_{f}$ by minimizing either
$L_{S}(\mathsf{F}_{\mathit{smt}}^{\phi}(h_{f}(\cdot)))$ with
$\tau\geq\nicefrac{{M^{2}B^{2}}}{{\epsilon^{2}}}$ iterations of stochastic
gradient descent, with gradients provided by
$\mathsf{B}_{\mathit{max}}^{\phi}$, and learning rate
$\eta=\sqrt{\nicefrac{{B^{2}}}{{M^{2}\tau}}}$ yields a grounding hypothesis
$\hat{h_{f}}\in\mathcal{H}_{f}$ that satisfies:
$\mathbb{E}[L_{\mathcal{D}}(\hat{h_{f}})]\leq\min_{h_{f}\in\mathcal{H}_{f}}L_{\mathcal{D}}(h_{f})+\epsilon$.
The randomness in this expectation is taken over the choices of the SGD
algorithm.

##  5 Experimental Evaluation

In this section we present an empirical evaluation of $\mathsf{SMTLayer}$ on
four learning problems that can be decomposed into perceptual and symbolic
subtasks. Our results demonstrate the following primary findings. _1)_
$\mathsf{SMTLayer}$ is effective: on every benchmark, it provides superior
results over “conventional” learning that takes place without encoded symbolic
knowledge. _2)_ $\mathsf{SMTLayer}$ has distinct advantages over prior
approaches. Compared with SATNet (Wang et al., 2019a), it requires
_significantly_ less training data to converge, and in all cases yields a more
accurate model; compared with Scallop (Huang et al., 2021), it is less
computationally expensive, requires less training data, and it is more
expressive in terms of the knowledge that it can encode. _3)_ Models trained
with $\mathsf{SMTLayer}$ may be more robust to certain types of covariate
shift that occur relative to the symbolic component of the problem; when
$\mathsf{SMTLayer}$ succeeds at learning a compatible representation, then it
will continue to produce correct inferences provided the perceptual component
remains stationary.

###  5.1 Datasets

Additional details on the datasets and corresponding architectures used in our
evaluation can be found in Appendix A.2, and specific hyperparameters used
when training on each dataset are in Appendix A.3.

#### MNIST Addition.

The MNIST addition problem is illustrated in Figure 1, and is similar to the
benchmark described by Huang et al. (2021). For training, we use “MNIST +$p$%”
to denote a training set of size 60,000 that contains $p\%$ of the possible
pairs of digits. So $p=100$ indicates all possible pairs of digits are used,
and for $p=10$, we only use pairs of the same digit. We use $p=10,25,50,75$
and $100$ in our experiments. In all cases, we use the same test set
consisting of instances from all possible pairs of digits.

#### Visual Algebra.

The task is to solve for the variable $x$ in a graphical depiction of the
equation $ax+b=c$, where $a,b$ and $c$ are randomly-chosen numbers, and each
symbol is depicted visually using EMNIST (Cohen et al., 2017a) and HASY
graphics (Thoma, 2017b). Similar to MNIST addition, the training sample
selects $a$ and $b$ uniformly from pairs of the same digit, and $x$ uniformly
from the odd numbers between 0 and 9. The test sample was generated by
sampling $a,b$ uniformly from all pairs of digits, and $x$ from all numbers 0
to 9.

#### Liar’s Puzzle.

The liar’s puzzle is comprised of three sentences spoken by three distinct
agents: Alice, Bob, and Charlie. One of the agents is “guilty” of an
unspecified offense, and in each sentence, the corresponding agent either
states that one of the other parties is either guilty or innocent. For
example, “Alice says that Bob is innocent.” It is assumed that two of the
agents are honest, and the guilty party is not. The solution to the problem is
an identification of the guilty party. A formal characterization of the
underlying logic is given in Appendix A.2. We note that the logic has non-
stratified occurrences of negation, so it cannot be encoded with Scallop. We
select a training sample that does not fully specify the logic, so
conventional training should be insufficient to identify a good model.

#### Visual Sudoku.

This task is to complete a $9\times 9$ Sudoku board where each entry is an
MNIST digit. We use the dataset from the SATNet evaluation (Wang et al.,
2019a), and examine three configurations obtained by sampling 10%, 50%, and
100% of the original training set. Although there are examples of Sudoku
solvers implemented as logic programs, we were not able to implement one in
Scallop without violating stratified negation. When calculating accuracy, we
check that the _entire_ Sudoku board is correct.

###  5.2 Setup

We implemented a prototype of our approach using Pytorch (Paszke et al., 2019)
and Z3 (De Moura & Bjørner, 2008), which will be made available in open-source
when this paper is published. When training models with $\mathsf{SMTLayer}$,
we use SGD with Nesterov momentum at rate $0.9$ and gradient clipping rate
0.1. Before training a model with $\mathsf{SMTLayer}$ (or a comparison
technique, unless stated otherwise), we first pre-train the neural network by
replacing $\mathsf{SMTLayer}$ with a dense network containing one hidden layer
of 512 neurons. This can potentially limit the number of training updates
needed at lower layers, but will not result in a model with a representation
that is compatible with symbolic knowledge, so further training is needed. The
models labeled “conventional” in our evaluation have the same architecture as
the one used for pre-training. Results were averaged over five runs of
training.

Our evaluation was performed on a machine with an Intel i9 1050K CPU, 64GB
memory, and a GeForce RTX 3080 accelerator running Ubuntu 20.04.4, with CUDA
11.1.0 and cuDNN 8.0.4. We developed and tested our prototype with Pytorch
version 1.7.0a0+7036e91 and Z3 4.8.14, and the results in our evaluation use
these versions as well.

| _Conventional_ |  _w/_ $\mathsf{SMTLayer}$ |  _w/_ SATNet |  _w/_ Scallop
---|---|---|---|---
_configuration_ | _test_ | _epoch_ | _test_ | _epoch_ | _test_ | _epoch_ | _test_ | _epoch_
| _acc. (%)_ | _time (sec.)_ | _acc.(%)_ | _time (sec.)_ | _acc.(%)_ | _time (sec.)_ | _acc. (%)_ | _time (sec.)_
MNIST+ 10% | 10.0 | 7.1 | 98.1 | 75.4 | 10.0 | 31.0 | 33.7 | 96.3
MNIST+ 25% | 32.5 | 7.1 | 98.3 | 74.8 | 34.2 | 30.9 | 65.8 | 96.4
MNIST+ 50% | 51.5 | 7.0 | 98.6 | 75.8 | 54.8 | 32.8 | 98.4 | 96.5
MNIST+ 75% | 76.1 | 7.0 | 98.5 | 75.0 | 78.4 | 31.9 | 93.5 | 96.4
MNIST+ 100% | 98.3 | 7.1 | 98.5 | 75.8 | 96.7 | 33.5 | 98.6 | 96.6
Vis. Alg. #1 | 24.1 | 13.2 | 98.2 | 168.2 | 19.6 | 80.1 | 18.7 | 602.8
Vis. Alg. #2 | 25.4 | 11.2 | 25.4 | 127.2 | 18.6 | 52.5 | 21.3 | 636.1
Liar’s Puzzle | 54.2 | 3.1 | 86.1 | 28.7 | 84.6 | 3.0 | — | —
Vis. Sudoku 10% | 0.0 | 6.3 | 66.0 | 135.7 | 0.0 | 9.9 | — | —
Vis. Sudoku 50% | 0.0 | 28.3 | 73.1 | 608.1 | 0.0 | 45.4 | — | —
Vis. Sudoku 100% | 0.0 | 26.7 | 79.1 | 1199.0 | 63.2 | 86.5 | — | —
Table 1:  Results after training and inference with $\mathsf{SMTLayer}$ versus
a conventional architecture. We use the publicly-available implementations of
SATNet (Wang et al., 2019a) and Scallop (Huang et al., 2021), with
hyperparameters matching those in their code. All $\mathsf{SMTLayer}$ test
accuracies were measured with the MaxSMT forward pass. Epoch times are
averaged over all epochs on which the model was trained. Cells marked — denote
that the problem is not compatible with the approach.

###  5.3 Results

#### Overall performance.

In terms of accuracy, Table 1 shows that $\mathsf{SMTLayer}$ outperforms both
the conventional network and prior work in terms of accuracy, training time,
or both, on all configurations. While training with $\mathsf{SMTLayer}$ (or
any of the above approaches) is more expensive than conventional,
$\mathsf{SMTLayer}$ is consistently faster than Scallop (nearly 4$\times$ in
the case of visual algebra). The per-epoch time to train the SATNet models is
less expensive than $\mathsf{SMTLayer}$, but this is not always conclusive. In
the case of visual sudoku, the 10% $\mathsf{SMTLayer}$ model achieved superior
error rates in 15 epochs, compared with 100 epochs for the 100% SATNet model;
this means that the $\mathsf{SMTLayer}$ model took less than one-tenth the
amount of time to train.

It is also worth noting that although Theorem 2 suggests that Algorithm 3
might have learning advantages over Algorithm 4, we found this not to be the
case on these datasets. All of the results in Table 1 were trained with
Algorithm 4, and test inference was done using Algorithm 1.

#### Training sample size.

Because $\mathsf{SMTLayer}$ encodes explicit knowledge that is essential to
correct inference on these datasets, our approach is able to perform well in
data-impoverished settings where the training sample is insufficient to fully
specify the symbolic component of the learning task. This is readily apparent
across the results in Table 1: in the MNIST addition and first visual algebra
configuration, $\mathsf{SMTLayer}$ yields a model that performs nearly
perfectly despite not being given a sufficient sample in most cases. Because
SATNet must learn the symbolic component, it is at a disadvantage, and in
these settings performs similarly to a conventional model. In theory, Scallop
should be able to perform as well as $\mathsf{SMTLayer}$, as it also encodes
explicit knowledge. However, it is unable to learn a useful model for either
visual algebra configuration, and does not learn the correct representation
for MNIST addition until it is exposed to half of the possible digit pairs
during training.

$\mathsf{SMTLayer}$ does particularly well on the visual Sudoku dataset
introduced by Wang et al.. When trained on just 10% of the original sample, it
learns a function that _exceeds_ the performance of the SATNet model by a
healthy margin, which continues to grow as it is exposed to more training
data. On the other hand, we found that SATNet failed to converge with less
than the full original training sample.

#### Robustness & interpretability.

The reason that $\mathsf{SMTLayer}$ is able to perform well, and often near
the optimum, in configurations that other approaches perform poorly on, is
that it learns a representation that is consistent with the symbolic knowledge
encoded in the $\mathsf{SMTLayer}$. For example, the constraints that we use
for MNIST addition, visual algebra, and visual sudoku all encode digits as
bitvectors. In order to make a correct inference, the neural network must
learn to encode MNIST digits in their correct bitvector representation. If
learning succeeds at this, then there are two positive outcomes that follow.
First, the model’s representation will be inherently interpretable, because it
will coincide with the provided symbolic domain knowledge, which is also
(presumably) interpretable. Second, the resulting model is naturally robust to
covariate shift that does not affect the distribution of perceptual data that
the network translates into theory symbols, but that does affect the
statistics of their composition.

This type of shift is on display in the MNIST 10% and visual algebra
experiments, where at training time, the model only sees pairs of same-
numbered digits, and at test time it is exposed to a substantially different
distribution of digit pairs or formulas. We verified this by examining the
representations learned by $\mathsf{SMTLayer}$ and Scallop on MNIST Addition
10%; it is unreasonable to expect SATNet to learn an interpretable
representation, as it is not provided with an interpretable theory during
training. As expected, $\mathsf{SMTLayer}$ produces the correct representation
at the rate of accuracy of a typical MNIST model ($\approx 99\%$), whereas
Scallop’s digit representation was correct roughly 50% of the time. However,
architecture plays a role in this robustness, as shown in the
$\mathsf{SMTLayer}$ results for the second visual algebra configuration.
Because the network is shown the full instance, and not the individual digits,
it learns the training bias. Despite having access to the symbolic formulas in
$\mathsf{SMTLayer}$, it cannot disentangle the perceptual symbols from their
covariance. Understanding this issue is an important direction for future
work.

##  6 Conclusion

Our approach for integrating logical theories into deep learning,
$\mathsf{SMTLayer}$, provides a pragmatic solution to the problem of
incorporating symbolic knowledge into learning for training and inference,
which we demonstrate on several problems involving both perceptual
tasks—vision and natural language—and logical reasoning. Notably, we show that
models which incorporate symbolic knowledge during training and inference can
outperform conventional models as well as prior work in this area, especially
in settings where training data is limited. Continued progress on automated
reasoning techniques has played a pivotal role in the development of several
fields over the past decades, and our hope is that the contributions in this
paper will aid in progress towards realizing their potential in challenges
that surpass the capabilities of existing learning techniques.

## References

  * Chang et al. (2020) Chang, O., Flokas, L., Lipson, H., and Spranger, M.  Assessing satnet’s ability to solve the symbol grounding problem.  _Advances in Neural Information Processing Systems_ , 33:1428–1439, 2020.
  * Chavira & Darwiche (2008) Chavira, M. and Darwiche, A.  On probabilistic inference by weighted model counting.  _Artificial Intelligence_ , 172(6-7):772–799, 2008.
  * Cohen et al. (2017a) Cohen, G., Afshar, S., Tapson, J., and Van Schaik, A.  Emnist: Extending mnist to handwritten letters.  In _2017 international joint conference on neural networks (IJCNN)_ , pp. 2921–2926. IEEE, 2017a.
  * Cohen et al. (2017b) Cohen, G., Afshar, S., Tapson, J., and van Schaik, A.  EMNIST: an extension of MNIST to handwritten letters.  _CoRR_ , abs/1702.05373, 2017b.  URL http://arxiv.org/abs/1702.05373.
  * De Moura & Bjørner (2008) De Moura, L. and Bjørner, N.  Z3: An efficient smt solver.  In _International conference on Tools and Algorithms for the Construction and Analysis of Systems_ , pp. 337–340. Springer, 2008.
  * Deng (2012) Deng, L.  The mnist database of handwritten digit images for machine learning research.  _IEEE Signal Processing Magazine_ , 29(6):141–142, 2012.
  * Fischer et al. (2019) Fischer, M., Balunovic, M., Drachsler-Cohen, D., Gehr, T., Zhang, C., and Vechev, M.  DL2: Training and querying neural networks with logic.  In _Proceedings of the 36th International Conference on Machine Learning_ , volume 97 of _Proceedings of Machine Learning Research_ , pp. 1931–1941. PMLR, 2019.
  * Hochreiter & Schmidhuber (1997) Hochreiter, S. and Schmidhuber, J.  Long short-term memory.  _Neural Comput._ , 9(8):1735–1780, nov 1997\.  ISSN 0899-7667.
  * Huang et al. (2021) Huang, J., Li, Z., Chen, B., Samel, K., Naik, M., Song, L., and Si, X.  Scallop: From probabilistic deductive databases to scalable differentiable reasoning.  _Advances in Neural Information Processing Systems_ , 34:25134–25145, 2021.
  * Koh et al. (2020) Koh, P. W., Nguyen, T., Tang, Y. S., Mussmann, S., Pierson, E., Kim, B., and Liang, P.  Concept bottleneck models.  In _International Conference on Machine Learning_ , pp. 5338–5348. PMLR, 2020.
  * Manhaeve et al. (2018) Manhaeve, R., Dumancic, S., Kimmig, A., Demeester, T., and De Raedt, L.  Deepproblog: Neural probabilistic logic programming.  _Advances in Neural Information Processing Systems_ , 31, 2018.
  * Paszke et al. (2019) Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga, L., Desmaison, A., Kopf, A., Yang, E., DeVito, Z., Raison, M., Tejani, A., Chilamkurthy, S., Steiner, B., Fang, L., Bai, J., and Chintala, S.  Pytorch: An imperative style, high-performance deep learning library.  In _Advances in Neural Information Processing Systems 32_. 2019.
  * Pennington et al. (2014) Pennington, J., Socher, R., and Manning, C. D.  Glove: Global vectors for word representation.  In _Proceedings of the 2014 conference on empirical methods in natural language processing (EMNLP)_ , pp. 1532–1543, 2014.
  * Pogancic et al. (2020) Pogancic, M. V., Paulus, A., Musil, V., Martius, G., and Rolínek, M.  Differentiation of blackbox combinatorial solvers.  In _8th International Conference on Learning Representations, ICLR 2020, Addis Ababa, Ethiopia, April 26-30, 2020_. OpenReview.net, 2020.  URL https://openreview.net/forum?id=BkevoJSYPB.
  * Riegel et al. (2020) Riegel, R., Gray, A., Luus, F., Khan, N., Makondo, N., Akhalwaya, I. Y., Qian, H., Fagin, R., Barahona, F., Sharma, U., et al.  Logical neural networks.  _arXiv preprint arXiv:2006.13155_ , 2020.
  * Shalev-Shwartz & Ben-David (2014) Shalev-Shwartz, S. and Ben-David, S.  _Understanding Machine Learning: From Theory to Algorithms_.  Cambridge University Press, 2014.  ISBN 1107057132.
  * Thoma (2017a) Thoma, M.  The hasyv2 dataset.  _CoRR_ , abs/1701.08380, 2017a.  URL http://arxiv.org/abs/1701.08380.
  * Thoma (2017b) Thoma, M.  Hasyv2 - handwritten symbol database, January 2017b.  URL https://doi.org/10.5281/zenodo.259444.
  * Topan et al. (2021) Topan, S., Rolnick, D., and Si, X.  Techniques for symbol grounding with satnet.  _Advances in Neural Information Processing Systems_ , 34:20733–20744, 2021.
  * van Krieken et al. (2022) van Krieken, E., Acar, E., and van Harmelen, F.  Analyzing differentiable fuzzy logic operators.  _Artificial Intelligence_ , 302:1–46, January 2022.
  * Varnai & Dimarogonas (2020) Varnai, P. and Dimarogonas, D. V.  On robustness metrics for learning stl tasks.  In _2020 American Control Conference (ACC)_ , pp. 5394–5399. IEEE, 2020.
  * Vlastelica et al. (2019) Vlastelica, M., Paulus, A., Musil, V., Martius, G., and Rolínek, M.  Differentiation of blackbox combinatorial solvers.  _arXiv preprint arXiv:1912.02175_ , 2019.
  * Wang et al. (2019a) Wang, P., Donti, P. L., Wilder, B., and Kolter, J. Z.  Satnet: Bridging deep learning and logical reasoning using a differentiable satisfiability solver.  In Chaudhuri, K. and Salakhutdinov, R. (eds.), _Proceedings of the 36th International Conference on Machine Learning, ICML 2019, 9-15 June 2019, Long Beach, California, USA_ , volume 97 of _Proceedings of Machine Learning Research_ , pp. 6545–6554. PMLR, 2019a.  URL http://proceedings.mlr.press/v97/wang19e.html.
  * Wang et al. (2019b) Wang, P.-W., Donti, P., Wilder, B., and Kolter, Z.  Satnet: Bridging deep learning and logical reasoning using a differentiable satisfiability solver.  In _International Conference on Machine Learning_ , pp. 6545–6554. PMLR, 2019b.
  * Yurichev (2020) Yurichev, D.  SAT/SMT by example, 2020.  Available at https://sat-smt.codes/ (January, 2023).

##  Appendix A Appendix

###  A.1 Proofs

######  Theorem 2.

Let $\mathcal{D}$, $\mathcal{H}$ be a $T$-decomposable problem with grounding
function $f$ and prediction logic $\phi$ where:

  1. 1.

$\mathbf{Z}$ and $\mathbf{Y}$ are Cartesian products of Booleans.

  2. 2.

For any $(x,y)\sim\mathcal{D}$ and $y^{\prime}\neq y$, $\phi(f(x),y^{\prime})$
is $T$-equivalent to false and there is exactly one $z$ such that $\phi(z,y)$
is $T$-equivalent to true.

  3. 3.

$\mathcal{H}_{f}$ is a convex set and for all $h_{f}\in\mathcal{H}_{f}$,
$\|h_{f}\|\leq B$, and the loss $\ell(h_{f}(\cdot),z)$ is $M$-Lipschitz and
convex in $x$ for any fixed $z$.

Then for any $\epsilon>0$, selecting $h_{f}$ by minimizing either
$L_{S}(\mathsf{F}_{\mathit{smt}}^{\phi}(h_{f}(\cdot)))$ with
$\tau\geq\nicefrac{{M^{2}B^{2}}}{{\epsilon^{2}}}$ iterations of stochastic
gradient descent, with gradients provided by
$\mathsf{B}_{\mathit{max}}^{\phi}$, and learning rate
$\eta=\sqrt{\nicefrac{{B^{2}}}{{M^{2}\tau}}}$ yields a grounding hypothesis
$\hat{h_{f}}\in\mathcal{H}_{f}$ that satisfies:
$\mathbb{E}[L_{\mathcal{D}}(\hat{h_{f}})]\leq\min_{h_{f}\in\mathcal{H}_{f}}L_{\mathcal{D}}(h_{f})+\epsilon$.
The randomness in this expectation is taken over the choices of the SGD
algorithm.

###### Proof.

To prove this result, we introduce the notion of a _grounding sample_.

######  Definition 4 (Grounding sample).

Let $\mathcal{D}$, $\mathcal{H}$ be a $T$-decomposable problem with grounding
function $f$. The grounding sample $S_{f}$ for $S\sim\mathcal{D}$ is given by
$[(x_{i},f(x_{i})):(x_{i},y_{i})\in S]$, i.e., tuples that consist of the
first element of each instance in $S$ and its image under $f$.

Now observe that the conditions stated in assumption (3) are sufficient to
yield the result if instead of optimizing
$L_{S}(\mathsf{F}_{\mathit{smt}}^{\phi}(h_{f}(\cdot)))$, we were given the
grounding sample $S_{f}$ and minimized $L_{S_{f}}(h_{f})$ (see (Shalev-Shwartz
& Ben-David, 2014), Theorem 14.8). The result follows as stated then because
of assumptions (1) and (2), which imply that the update vectors provided by
$\mathsf{B}_{\mathit{max}}^{\phi}$ are the gradients of $L_{S_{f}}(h_{f})$.

To understand why, observe that the sign of $\hat{y}$ computed on line 3 of
both algorithms must be equal to that of $y^{\star}$. This follows from two
facts:

  1. 1.

At any coordinate $i$ where $y[i]\neq y^{\star}[i]$,
$\mathsf{sign}(\partial_{y}\ell(y,y^{\star}))[i]=\mathsf{sign}(y)[i]$.

  2. 2.

At any coordinate $i$ where $y[i]=y^{\star}[i]$,
$\mathsf{sign}(\partial_{y}\ell(y,y^{\star}))[i]=-1\cdot\mathsf{sign}(y)[i]$.

Now there are two cases to consider.

#### Case 1: $\mathsf{sign}(y)=\mathsf{sign}(\hat{y})$.

In this case the algorithm returns $\partial_{z}\ell(z,\mathsf{sign}(z))$.
Because solutions to $\phi(\cdot,y)$ are unique, $\mathsf{sign}(z)$ is the
correct grounding, i.e., $z=f(x)$ for the original features $x$.

#### Case 2: $\mathsf{sign}(y)\neq\mathsf{sign}(\hat{y})$.

Because of assumption (2), the set of indices computed on line 8 will contain
all of the coordinates at which $z$ matches the correct value
$z^{\star}=f(x)$. Note that at these coordinates, the vector returned by the
algorithm matches the gradient of $\ell(z,z^{\star})$, which is
$\partial_{z[i]}\ell(z[i],\mathsf{sign}(z[i]))$. In the remaining coordinates,
the vector will contain $\partial_{z[i]}\ell(z[i],1-\mathsf{sign}(z[i]))$,
which also matches the gradient of $\ell(z,z^{\star})$. The result follows. ∎

###  A.2 Dataset details

Two of the three problems that we examine are based on the MNIST handwritten
digit dataset (Deng, 2012), which consists of 60,000 28x28 gray-scale images
of handwritten numerals for training and 10,000 instances for testing. The
digits on the left of Figure 1 are examples of instances from this data. To
generate data for the visual algebra problem, we additionally drew from EMNIST
(Cohen et al., 2017b), which extends MNIST with handwritten letters, and HASY
(Thoma, 2017a), which contains handwritten symbols with similar
characteristics to MNIST. For the liar’s puzzle, inspired by examples
(Yurichev, 2020) which formulate similar examples as SMT constraints, we
constructed examples using a set of common phrases that we devised ourselves,
and did not otherwise draw from publicly-available data.

Below, we describe the way in which we used these data sources to construct
training and test samples, and the neural network architectures that we used
with $\mathsf{SMTLayer}$ to solve them.

#### MNIST Addition.

The MNIST addition problem is described in Example 1. In each instance, two
MNIST digits are presented as features, and the task of the model is to
provide their sum represented as a bitvector. The architecture that we use
consists of four convolutional layers with kernels of size 3, depths in the
order 64, 64, 128, 128, and a stride of width 2 on the first layer, and two
dense layers of width 256 and 4. This network is applied to each digit, and
the results are concatenated to obtain a vector of size 8 that is passed to an
instance of $\mathsf{SMTLayer}$ with SMT constraints from Example 1, which
ultimately produces a vector of width 5 that represents the bitvector sum of
the digits.

We generated five training samples starting with one containing only pairs of
the same digit, i.e. $(0,0),(1,1),\ldots$. We then added progressively more
from the full set of possible pairs, using 25%, 50%, and 100%, and trained on
batches of 128 across all datasets. Note that although we change the number of
digit pairs that appear between samples, we always map these pairs to random
MNIST images to obtain 60,000 training instances. This is to ensure that the
training sample contains a sufficient sample of MNIST images to be able to
perform well on the test data. In all cases, we use the same test set
consisting of instances from all possible pairs of digits. The purpose of this
is to demonstrate that the conventional network will not generalize until it
has seen the full distribution, whereas the model with $\mathsf{SMTLayer}$
should be able to generalize after seeing many fewer examples.

#### Visual Algebra.

The visual algebra problem is described in Example 2 and Example 3. Recall
that features depict handwritten depictions of linear equations of the form
$ax+b=c$. Values for $a,x$ and $b$ are randomly drawn from the range 1-9 to
ensure that solutions are unique. Then the corresponding value of $c$ is
decomposed into $c=10\cdot c_{1}+c_{2}$, and MNIST digits are selected at
random to represent $a,b,c_{1},c_{2}$. A random EMNIST alphabet character is
drawn for the variable, and random multiplication, addition, and equality
symbols are drawn from HASY. A minor note is that HASY does not contain the
standard equality symbol “$=$”, so we instead use “$\doteq$”. These images are
then concatenated horizontally in the appropriate order.

We evaluate two architectures for this problem. The first uses the same neural
network that was used for MNIST addition, and an instance of
$\mathsf{SMTLayer}$ with the constraints given in Example 2. It assumes that
the four numeric digits in the problem have already been extracted, e.g. by a
separate vision routine that can recognize digits from letters and arithmetic
symbols, and are provided directly to the model. The four inputs are given
separately to the neural network, which produces four 4-bit bitvectors that
are concatenated and passed to $\mathsf{SMTLayer}$, which produces a 4-bit
bitvector result. We refer to this as configuration #1 in our results.

The second uses an architecture which takes the entire image containing the
problem as a whole, and produces a 16-bit bitvector that is passed directly to
$\mathsf{SMTLayer}$. This architecture uses a similar stack of convolutional
layers, but has a larger initial dense layer containing 26,112 neurons, as it
is given a larger image. The difference in the convolutional stack is at the
second layer, which also has a stride of width 2, to reduce the size of the
feature map and mitigate the need for an even larger dense connection. The
difference between these architectures relates to one of the challenges of
this problem. Much of the information contained in the features is irrelevant
to the solution, e.g., it is irrelevant which letter is chosen for the
variable, or what the arithmetic operators look like, so this architecture
must also learn to disregard these parts of the instance. We refer to this as
configuration #2 in our results and in Example 3.

We generated a training sample by selecting $a$ and $b$ uniformly from pairs
of the same digit, i.e. $(0,0),(1,1),\ldots$, and sampling $x$ uniformly from
the odd numbers between 0 and 9. The test sample was generated by sampling
$a,b$ uniformly from all pairs of digits, and $x$ from all numbers 0 to 9. We
then map these values to random MNIST, EMNIST, and HASY images to obtain
60,000 samples. The intention is to study a problem wherein the model is not
shown all possible problems (modulo representation as digits), or all of the
solutions. This is more challenging than MNIST addition for two reasons: for a
given solution, there are many more compatible ground terms, and the model
does not see examples of some of the solutions it must provide for the test
set. Thus, in order for $\mathsf{SMTLayer}$ to succeed, it must use the
provided symbolic knowledge to approximate the correct grounding function,
despite these deficiencies in the data.

#### Liar’s Puzzle.

The liar’s puzzle is comprised of three sentences spoken by three distinct
agents: Alice, Bob, and Charlie. One of the agents is “guilty” of an
unspecified offense, and in each sentence, the corresponding agent either
states that one of the other parties is either guilty or innocent. It is
assumed that two of the agents are honest, and the guilty party is not. The
solution to the problem is an identification of the guilty party. An example
is described in Example 4

We synthesized a dataset for the liar’s puzzle based on a limited set of
utterances about who speaks in each sentence, which agent is the subject, and
whether the subject is guilty or innocent. There are five ways of denoting the
speaker: _“* says”_ , _“* says that”_ , _“* said”_ , _“* said that”_ , and a
colon _“* :”_ separating the speaker’s name from the rest of the sentence.
There are five ways of uttering either innocence or guilt: _“* did it/did not
do it”_ , _“* is guilty/innocent”_ , _“* is definitely guilty/innocent”_ , and
_“* definitely did it/did not do it”_ , _“* is the criminal/is a good
person”_. We generated all of the combinations of subject, speaker, and
proclaimed innocence or guilt, and took the product with all possible
combinations of these utterances. The result is a dataset of 375,000
instances, each containing three natural language sentences.

The prediction logic for this problem assumes a set of ground predicates
$\mathsf{speaker}(\mathit{agent},\mathit{sentence})$,
$\mathsf{subject}(\mathit{agent},\mathit{sentence})$,
$\mathsf{accuses}(\mathit{sentence})$, and $\mathsf{guilty}(\emph{agent})$.
For example, if the first sentence was _“Alice says that Bob is innocent”_ ,
the ground predicates would be $\mathsf{speaker}(\mathsf{alice},1)$,
$\mathsf{subject}(\mathsf{bob},1)$, and $\lnot\mathsf{guilty}(bob)$. Then the
prediction logic is shown in Equation 1.

$\begin{array}[]{ll}&|\\{a:\mathsf{guilty}(a)\\}|=1\\\ \land&\forall
a.|\\{s:\mathsf{speaker}(a,s)\\}|=1\\\ \land&\forall
s.|\\{a:\mathsf{subject}(a,s)\\}|=1\\\ \land&\forall
s,a_{1},a_{2}.\mathsf{speaker}(a_{1},s)\land\mathsf{subject}(a_{2},s)\rightarrow\\\
&\ \ \ \
\textbf{}(\lnot\mathsf{guilty}(a_{1})\leftrightarrow\mathsf{accuses}(s)\leftrightarrow\mathsf{guilty}(a_{2}))\par\end{array}$
(1)

In our implementation, agents and sentence identifiers are encoded unary as
3-bit bitvectors. The quantifiers are removed by substituting for each ground
term or sentence identifier, and the cardinality constraints are expanded into
propositional formulas. The architecture that we adopt is a two-layer
bidirectional long short-term model (LSTM) with 512 dimensions at each layer,
and a 300-dimension trainable embedding layer initialized from GloVe-6B
(Pennington et al., 2014). The hidden units of the last LSTM layer were
connected to 2-layer dense network containing 128 followed by 7 neurons. This
network is applied to each sentence in the input, and the concatenated results
are passed to the $\mathsf{SMTLayer}$, which solves the formula in Equation 1
to produce a unary encoding of the guilty party.

To evaluate solver layers on this problem, we selected training and test
samples by first subsampling half of the full 375,000 available instances. We
then selected half of the $\mathsf{speaker}$, $\mathsf{subject}$,
$\mathsf{accuses}$ predicate configurations for all three sentences appearing
in this subsample to appear in the training sample, and the other half to
appear in the test sample. To further limit the amount of information in the
training sample, we randomly selected one ordering for each predicate
configuration to remain for training. There were 9,400 resulting training
instances, and 28,062 test instances. Restricting the training set as
described ensures that the model is trained on a limited subset of possible
sentence configurations, and one that is logically disjoint from those that
appear in the test sample. Because there is not enough information in the
training sample to learn Equation 1, we expect only the model with
$\mathsf{SMTLayer}$ to succeed, but to do so it must approximate the grounding
function well from a limited sample.

Features Symbolic Domain $\mathbf{Z}$ $\mathsf{0010000111}$ Grounding Function
$f$ $\begin{array}[]{ll}\phi(\quad z_{1}\|\ldots\|z_{10},\quad y\quad)\
\equiv&\\\ a=\mathds{1}_{z_{1}>0}\|\ldots\|\mathds{1}_{z_{5}>0}\land
b=\mathds{1}_{z_{5}>0}\|\ldots\|\mathds{1}_{z_{10}>0}&\land\\\
c=\mathds{1}_{z_{5}>0}\|\ldots\|\mathds{1}_{z_{10}>0}&\land\\\
ax+b=c&\end{array}$ Prediction Logic $\phi$ Labels $\mathbf{Y}$
$\\{\mathsf{00101}\\}$ Satisfying Assignments Figure 2:  Visual Algebra
configuration 1 example. Features Symbolic Domain $\mathbf{Z}$
$\mathsf{0100010000100100}$ Grounding Function $f$
$\begin{array}[]{ll}\phi(\quad z_{1}\|\ldots\|z_{10},\quad y\quad)\ \equiv&\\\
a=\mathds{1}_{z_{1}>0}\|\ldots\|\mathds{1}_{z_{5}>0}\land
b=\mathds{1}_{z_{5}>0}\|\ldots\|\mathds{1}_{z_{10}>0}&\land\\\
ax+b=c&\end{array}$ Prediction Logic $\phi$ Labels $\mathbf{Y}$
$\\{\mathsf{00101}\\}$ Satisfying Assignments Figure 3:  Visual Algebra
configuration 2 example. InputAlice said that Bob did not do it.Bob: Alice is
definitely innocent Charlie: Alice did it  Symbolic Domain $\mathbf{Z}$
$\mathsf{100010001010100100110}$ Grounding Function $f$
$\begin{array}[]{ll}&|\\{a:\mathsf{guilty}(a)\\}|=1\\\ \land&\forall
a.|\\{s:\mathsf{speaker}(a,s)\\}|=1\\\ \land&\forall
s.|\\{a:\mathsf{subject}(a,s)\\}|=1\\\ \land&\forall
s,a_{1},a_{2}.\mathsf{speaker}(a_{1},s)\land\mathsf{subject}(a_{2},s)\rightarrow\\\
&\ \ \ \
\textbf{}(\lnot\mathsf{guilty}(a_{1})\leftrightarrow\mathsf{accuses}(s)\leftrightarrow\mathsf{guilty}(a_{2}))\end{array}$
Prediction Logic $\phi$ Labels (Liar) $\mathbf{Y}$ $\\{\mathsf{2}\\}$
Satisfying Assignments Figure 4:  Liar’s Puzzle example.

###  A.3 Hyperparameters

| _Conventional_ |  _w/_ $\mathsf{SMTLayer}$ |  _w/_ SATNet |  _w/_ Scallop
---|---|---|---|---
| _optimizer_ | _epochs_ | _optimizer_ | _epochs_ | _optimizer_ | _epochs_ | _optimizer_ | _epochs_
MNIST+ 10% | SGD(1.0) | 0/5 | SGD(1.0) | 3/5 | Adam(2e-3, 1e-5) | 3/5 | Adam(1e-3) | 3/5
MNIST+ 25% | SGD(1.0) | 0/5 | SGD(1.0) | 3/5 | Adam(2e-3, 1e-5) | 3/5 | Adam(1e-3) | 3/5
MNIST+ 50% | SGD(1.0) | 0/5 | SGD(1.0) | 3/5 | Adam(2e-3, 1e-5) | 3/5 | Adam(1e-3) | 3/5
MNIST+ 75% | SGD(1.0) | 0/5 | SGD(1.0) | 3/5 | Adam(2e-3, 1e-5) | 3/5 | Adam(1e-3) | 3/5
MNIST+ 100% | SGD(1.0) | 0/5 | SGD(1.0) | 3/5 | Adam(2e-3, 1e-5) | 3/5 | Adam(1e-3) | 3/5
Vis. Alg. #1 | SGD(1.0) | 0/5 | SGD(1.0) | 3/5 | Adam(2e-3, 1e-5) | 3/5 | Adam(1e-3) | 3/5
Vis. Alg. #2 | SGD(1.0) | 0/5 | SGD(1.0) | 3/5 | Adam(2e-3, 1e-5) | 3/5 | Adam(1e-3) | 3/5
Liar’s Puzzle | Adam(2e-3) | 0/15 | Adam(1e-3) | 15/5 | Adam(2e-3, 1e-5) | 15/5 | — | —
Vis. Sudoku 10% | SGD(1.0) | 0/100 | SGD(1.0) | 30/15 | Adam(2e-3, 1e-5) | 30/100 | — | —
Vis. Sudoku 50% | SGD(1.0) | 0/100 | SGD(1.0) | 30/5 | Adam(2e-3, 1e-5) | 30/100 | — | —
Vis. Sudoku 100% | SGD(1.0) | 0/100 | SGD(1.0) | 30/5 | Adam(2e-3, 1e-5) | 0/100 | — | —
Table 2:  Training hyperparameters. Numbers in parentheses after the optimizer
denote the learning rate; when there are multiple numbers, different learning
rates were applied to different parameter groups, as detailed in the text.
Each epoch pair corresponds to the pre-training and training phases, i.e., 3/5
denotes three epochs of pre-training and five subsequent training epochs with
the solver layer.

The four problems that our evaluation studies vary considerably in size and
complexity, and the models used to train them require different
considerations. This section details these differences. Table 2 relates the
optimizers and epochs for each dataset and solver layer configuration. For
SATNet and Scallop, we use the optimization settings described in their
respective papers, and found in their public implementations. For SATNet
models, the MaxSAT clause parameters were trained at a rate of 2e-3, and the
convnet at rate 1e-5. When SGD(1.0) is stated, we used a warmup period
spanning the first epoch, and cosine annealing for the remainder of training.

#### MNIST Addition.

The MNIST addition problem is the easiest of the problems that we study, at
least in its full (100%) configuration. We find that for the 100%
configuration, all of the solver layers converge to a nearly optimal solution
with three epochs of supervised pre-training, and five epochs of subsequent
training with the solver layer attached. For the subsample configurations, all
of the solver layers converge to a stable, although in many cases suboptimal,
solution within these parameters as well. For $\mathsf{SMTLayer}$, we clipped
gradients for all parameters at 0.1, and did not clip gradients for the other
solver layer models. For all configurations, we used batches of size 128.

#### Visual Algebra.

Although visual algebra is a more difficult learning problem than MNIST
addition, as evidenced by the results in Table 1, we find that the same
parameters allow all of the configurations studied in our evaluation to
converge. After 3/5 epochs of training, the models stabilize, and in some
cases, further training yields an overfit model. For $\mathsf{SMTLayer}$, we
clipped gradients for all parameters at 0.1, and did not clip gradients for
the other solver layer models. For all configurations, we used batches of size
64.

#### Liar’s Puzzle.

The Liar’s Puzzle is the only problem to use a recurrent model, and we found
that it required more epochs of pre-training to reduce the variance of the
final model with the solver layer. Additionally, the SGD optimizer used with
$\mathsf{SMTLayer}$ on other datasets caused the model to converge at local
minima. We found that pre-training at a higher learning rate, and using Adam
with a default learning rate, let to the best results. Additionally, we did
not clip gradients for the $\mathsf{SMTLayer}$ model. We used the same
parameters, but the normal SATNet optimizer, for the SATNet model. For all
configurations, we used batches of size 32.

#### Visual Sudoku.

Visual Sudoku is the most challenging problem that we studied, for all solver
layers as well as the conventional model. We did not use supervised pre-
training, as the supervision in this problem leaks the correct labels directly
to the model, bypassing the solver layer and the need for its updates (Chang
et al., 2020). We instead used the unsupervised pre-training method described
in (Topan et al., 2021), and found that 30 epochs of unsupervised pre-training
was sufficient to yield consistent and quick convergence with the solver
layer. For the most data-scarce configuration (10%), 15 epochs of training
with the solver layer were needed to converge, and for the others five epochs
were sufficient. For the $\mathsf{SMTLayer}$ model, we used batches of size 1
after pretraining (batch size 64 during pre-training), primarily due to the
fact that $\mathsf{SMTLayer}$ returns only the indices masked as non-hint
elements on the Sudoku board. Because each instance has a different number of
hint elements, this would lead to ragged tensors during training, which
Pytorch does not support.

When assessing SATNet on Visual Sudoku, we were unable to converge to a useful
model on any except the full (100%) configuration, as discussed in Section
5.2. Using the authors’ public implementation and training script, we
attempted the 10% and 50% configurations with and without unsupervised pre-
training, with and without the measures taken in (Chang et al., 2020) to
prevent label leakage, and with batches of size 40 (as used in the original
paper) as well as 1, to no avail. For the 100% configuration, we were able to
reproduce a useful model; we use the accuracy reported in the original paper
in Table 1 for consistency, as the average that we obtained did not differ
significantly from this.


# Virial clouds to explain rotational asymmetry in galactic halos

Asghar Qadir1,∗; Noraiz Tahir1,†; Francesco De Paolis2,3,‡;
Achille A. Nucita2,3,§
1 School of Natural Sciences, National University of Science and Technology,
H-12,
Islamabad, 46000, Pakistan
2 Department of Mathematics and Physics “E. De Giorgi”, University of Salento,
Via per Arnesano, I-73100 Lecce, Italy.
3 INFN, Sezione di Lecce, Via per Arnesano, 73100 Lecce, Italy.
∗<EMAIL_ADDRESS>
†<EMAIL_ADDRESS>
‡<EMAIL_ADDRESS>
§<EMAIL_ADDRESS>

###### Abstract

The rotation of the galactic objects has been seen by asymmetric Doppler shift
in the CMB data. Molecular hydrogen clouds at virial temperature may
contribute to the galactic halo dark matter and they might be the reason for
the observed rotational asymmetry in the galactic halos. We present a method
to constrain the parameters of these virial clouds given that they are
composed of a single fluid. The method is such that it should be possible to
extend it to more than one fluid.

##  1 Introduction

Half of the baryonic matter in the Universe is dark [1, 2, 3, 4], some of it
may consist of molecular hydrogen clouds left over from the process of star
formation, which would continue to collapse till their virial temperature
reaches the CMB. To observe them one could look for $\gamma$-rays
scintillation due to strikes by cosmic rays on ${\rm H_{2}}$ molecules. A
$\gamma$-ray corona around the Milky Way was discovered [5] with the expected
angular distribution [6], but it was not not clear how much, if any, of the
carpet might be due to the proposed clouds. An alternative suggestion was to
look for an asymmetric Doppler shift of these clouds in M31 [7].

The predicted effect was seen at a sufficient but low level in the WMAP data
in 2011 [8] and confirmed in the Planck data in 2014 [9]. Objections that the
blue shift may be due to “hot spots” were countered by the argument that the
red shift could not be “cold spots” and were soon dispelled. The precision was
sufficient to be able to use for studying the rotation of the M31. This was
followed by the same effect being seen in NGC 5128, M33, M82, M81 [10, 11, 12,
13]. However, the verification of the effect that it is due to the molecular
clouds which led to its prediction, does not prove that it is due to the cause
that led to the prediction. It could be partially or totally due to other
interstellar matter. To be able to determine the nature of the cause, we need
to model the molecular hydrogen clouds or dust clouds produced by the
mechanism mentioned, or some mixture of the two, and compare with the data.
One could then try to look at whether contamination of the spectrum by
radiation from higher temperature sources could be distinguished.

##  2 The Virial Cloud Model

As these clouds are immersed in a heat bath of the CMB, they are isothermal.
Since they have formed due to the Jeans instability, they must have a
precisely defined Jeans radius and consequent Jeans mass. The density profile
would then be given by the solution of the Lane-Emden equation subject to the
boundary conditions that the density be literally zero (not approximately
zero) at the boundary and the density profile be flat at the core. The latter
condition says that the derivative of the density be zero at the core. We call
them “virial clouds”. Since James Jeans had applied the analysis for the
formation of a normal star, he had taken the boundary condition that the
density of the cloud merge with that of the interstellar medium (ISM) (which
has been called “the Jeans fiddle”) [14]. In that case the central density had
to be put in on an ad hoc basis. We do not need to do this. Further, we will
want to extend the analysis (in later work) to a two-fluid model. As such, we
obtain the equations from first principles in any case.

For an isothermal gas sphere, the temperature, $T$, is directly related to the
mass of the molecule, $m_{H}$ in our case, and the average speed. The speed of
the perturbation would have to be taken to be the speed of sound in the cloud,
$c_{s}$. This gives

$c_{s}=\sqrt{\frac{\gamma kT_{CMB}}{m_{H}}}\leavevmode\nobreak\ ,$ (1)

where $\gamma$ is the adiabatic factor given by the ratio of the heat capacity
at constant pressure to that at constant volume. At the densities and
temperatures that we are considering no higher degrees of freedom would be
excited and we would only have the translational degree, making the ideal gas
approximation extremely good and so $\gamma=5/3$. For our purpose,
$m_{H}\approx 2.016\leavevmode\nobreak\ {\rm g/mol}$ and $T\approx
2.726^{o}K$, so

$c_{s}\simeq 1.110\times 10^{4}\leavevmode\nobreak\ {\rm
km\leavevmode\nobreak\ s^{-1}}\leavevmode\nobreak\ .$ (2)

Now the virial theorem $2K+U=0$ gives the Jeans mass as [15]

$M_{J}^{2}\simeq\left(\frac{81}{32\pi\rho_{c}}\right)\left(\frac{3c_{s}^{2}}{5G}\right)^{3}\leavevmode\nobreak\
,$ (3)

and the corresponding Jeans radius as

$R_{J}^{2}=\frac{27c_{s}^{2}}{20\pi\rho_{c}G}\leavevmode\nobreak\ ,$ (4)

where, $\rho_{c}$ is the central density of the virial cloud, and $G$ is
Newton’s gravitational constant. We need to use the canonical ensemble
distribution in order to calculate the density profile of the clouds which is
given by the relation [16]

$f(r,p)=\frac{1}{h^{3N}N!}\frac{1}{Z}exp\left(-\frac{H(r,p)}{kT_{CMB}}\right)\leavevmode\nobreak\
,$ (5)

where, $H(r,p)$ is the Hamiltonian of the system. The partition function $Z$
is the normalization factor, thus being the integral over phase space of the
rest of this distribution function. The Hamiltonian is position dependent

$H(r,p)=\frac{p^{2}}{2m_{H}}+\frac{GM(r)m_{H}}{r}\leavevmode\nobreak\ ,$ (6)

and so the density distribution, is given by

$\displaystyle\rho(r)=8m_{H}^{5/2}\left(\frac{G\leavevmode\nobreak\
\rho_{c}}{3k\leavevmode\nobreak\
T_{CMB}}\right)^{3/2}exp\left(-\frac{GM(r)m_{H}}{rkT_{CMB}}\right),$ (7)

where $M(r)$ is the mass of the cloud interior to $r$ given by

$GM(r)=4\pi\int_{0}^{r}\rho(q)q^{2}dq\leavevmode\nobreak\ .$ (8)

Our problem, then, is that the density is defined in terms of the integral of
the density, i.e. we have an integral equation. Taking the natural logarithm
of eq.(7) we obtain

$\int_{0}^{r}q^{2}\rho(q)dq=-\frac{rkT_{CMB}}{4\pi
Gm_{H}}\ln\left(\frac{\rho_{(}r)}{\zeta}\right)\leavevmode\nobreak\ ,$ (9)

where $\zeta=8m_{H}^{5/2}(G\rho_{c}/3kT_{CMB})^{3/2}$. Taking the derivative
of the above equation with respect to $r$, and substituting eq.(7) we obtain
the ordinary differential equation

$r\frac{d\rho(r)}{dr}-r^{2}\left(\frac{4\pi
Gm_{H}}{kT_{CMB}}\right)\rho^{2}(r)-\rho(r)\ln\left(\frac{\rho(r)}{\zeta}\right)=0\leavevmode\nobreak\
.$ (10)

It will be noticed that this is a nonlinear differential equation that cannot
be solved by separation of variables. However, we do have the initial
condition $\rho^{\prime}(r)|_{r=0}=0$ and the boundary condition that
$\rho(R_{J})=0$. As such, we can solve it numerically by assuming a central
density and the initial condition. If the density at the boundary is too high,
we lower the central value and try again, till it becomes too low. At that
stage we raise the central density a bit till it again becomes too high.
Iteratively we bring the boundary density close enough to zero to be
acceptable. The resulting density profile is shown in Fig.1. One sees that the
central density of the cloud is ${\displaystyle\rho_{c}\simeq 1.60\times
10^{-18}\leavevmode\nobreak\ kgm^{-3}}$, the radius is $R_{J}\simeq 0.030{\rm
pc}$, and the mass is $M_{J}\simeq 0.798M_{\odot}$.

Figure 1: The density of the virial molecular hydrogen cloud. It decreases
monotonically from the central density, ${\displaystyle\rho_{c}\simeq
1.60\times 10^{-18}\leavevmode\nobreak\ kgm^{-3}}$, to zero at $r=R_{J}\simeq
0.030\leavevmode\nobreak\ {\rm pc}$.

##  3 Luminosity of the Clouds

Since the virial clouds are at the CMB temperature, their luminosity is the
CMB luminosity, except for their Doppler shift. Thus, we only need to look at
the differential frequency shift, $\Delta\nu/\nu=-(v_{rot}/c)\cos\theta$,
where $v_{rot}$ is the rotational velocity of the galactic halo, $c$ the speed
of light and $\theta$ the angle between the direction of the velocity and the
line of sight. Now $\Delta\nu/\nu=\Delta E/E=\Delta T/T$, so the differential
energy is given by $\Delta E=(\Delta T/T)E$, where $E$ is the energy in the
radiation for a patch of the sky corresponding to the size of the cloud. Now,
as the luminosity, $L$, is proportional to the energy, we get the luminosity
of the CMB for that patch. Then the luminosity of the cloud will be [17]

$L_{c}=(\Delta\nu/\nu)L\simeq 1.486\times
10^{-4}L\cos\theta\leavevmode\nobreak\ .$ (11)

Thus, as we scan across the halo, the luminosity will not change linearly but
with $\cos\theta$. We can use this check for the Doppler shift explanation and
obtain the energy associated with each cloud. This would directly give an
estimate of the total number of clouds being seen.

##  4 Results and Discussions

The virial clouds are completely constrained by physical requirements, and so
precisely determinable. They provide a satisfactory explanation of the
observed asymmetry of the CMB Doppler shift in galactic halos. We have only
considered pure molecular hydrogen clouds, but we also need to model these
molecular hydrogen clouds contaminated with more or less interstellar dust
which will be dealt with elsewhere. We have not discussed the possible
contamination of radiation by matter at high temperature. That needs to be
done. The effect is probably negligible, but cannot be neglected without
checking. One also needs to investigate how the clouds would evolve over
billions of years, from the time when the CMB temperature was higher than it
is now to the present value. Precisely how did the clouds lose energy over
this period?

## Acknowledgments

We are grateful to the Salento University of Lecce for hospitality at Lecce,
where part of this work was done, and AQ is grateful to IUPAP, ESA, ICRANet
and AS-ICTP for support to participate in MG15.

## References

  * [1] G. Bertone, Particle dark matter: observations, models and searches (Cambridge University Press, 2010).
  * [2] R. Cen and J. P. Ostriker, Where are the baryons?, The Astrophysical Journal 514, p. 1 (1999).
  * [3] R. Cen and J. P. Ostriker, Where are the baryons? ii. Feedback effects, The Astrophysical Journal 650, p. 560 (2006).
  * [4] A. Fraser-McKelvie, K. A. Pimbblet and J. S. Lazendic, An estimate of the electron density in filaments of galaxies at $z\sim 0.1$, Monthly Notices of the Royal Astronomical Society 415, 1961 (2011).
  * [5] D. Dixon, D. Hartmann, E. Kolaczyk, J. Samimi, R. Diehl, G. Kanbach, H. Mayer-Hasselwander and A. Strong, Evidence for a Galactic gamma-ray halo, New Astronomy 3, 539 (1998).
  * [6] F. De Paolis, G. Ingrosso, P. Jetzer and M. Roncadelli, Gamma-ray astronomy and baryonic dark matter, The Astrophysical Journal Letters 510, p. L103 (1998).
  * [7] F. De Paolis, G. Ingrosso, P. Jetzer, A. Quadir and M. Roncadelli, Observing molecular hydrogen clouds and dark massive objects in galactic halos, Astronomy & Astrophysics 299, p. 647 (1995).
  * [8] F. De Paolis, V. Gurzadyan, G. Ingrosso, P. Jetzer, A. Nucita, A. Qadir, D. Vetrugno, A. Kashin, H. Khachatryan and S. Mirzoyan, Possible detection of the M31 rotation in WMAP data, Astronomy & Astrophysics 534, p. L8 (2011).
  * [9] F. De Paolis, V. Gurzadyan, A. Nucita, G. Ingrosso, A. Kashin, H. Khachatryan, S. Mirzoyan, E. Poghosian, P. Jetzer, A. Qadir et al., Planck confirmation of the disk and halo rotation of M31, Astronomy & Astrophysics 565, p. L3 (2014).
  * [10] F. De Paolis, V. Gurzadyan, A. Nucita, G. Ingrosso, A. Kashin, H. Khachatryan, S. Mirzoyan, G. Yegorian, P. Jetzer, A. Qadir et al., Planck revealed bulk motion of Centaurus a lobes, Astronomy & Astrophysics 580, p. L8 (2015).
  * [11] V. Gurzadyan, F. De Paolis, A. Nucita, G. Ingrosso, A. Kashin, H. Khachatryan, S. Sargsyan, G. Yegorian, P. Jetzer, A. Qadir et al., Planck view of the M82 galaxy, Astronomy & Astrophysics 582, p. A77 (2015).
  * [12] F. De Paolis, V. Gurzadyan, A. Nucita, L. Chemin, A. Qadir, A. Kashin, H. Khachatryan, S. Sargsyan, G. Yegorian, G. Ingrosso et al., Triangulum galaxy viewed by Planck, Astronomy & Astrophysics 593, p. A57 (2016).
  * [13] V. Gurzadyan, F. De Paolis, A. Nucita, A. Kashin, A. Amekhyan, S. Sargsyan, G. Yegorian, A. Qadir, G. Ingrosso, P. Jetzer et al., Messier 81’s Planck view versus its halo mapping, Astronomy & Astrophysics 609, p. A131 (2018).
  * [14] J. H. Jeans, The stability of a spherical nebula, Proceedings of the Royal Society of London 68, 454 (1901).
  * [15] S. Chandrasekhar, An introduction to the study of stellar structure (Courier Corporation, 1957).
  * [16] P. M. Morse and B.W. Shore, Thermal physics, Physics Today 17, p. 60 (1964).
  * [17] N. Tahir, F. De Paolis, A. Qadir and A. Nucita, Constraining baryons in the M31 halo using Planck data, International Journal of Modern Physics D, accepted in (2019).