Neural-symbolic AI (NeSy) allows neural networks to exploit symbolic background knowledge in the form of logic. It has been shown to aid learning in the limited data regime and to facilitate inference on out-of-distribution data.
Probabilistic NeSy focuses on integrating neural networks with both logic and probability theory, which additionally allows learning under uncertainty.
A major limitation of current probabilistic NeSy systems, such as , is their restriction to finite probability distributions, i.e., discrete random variables.
In contrast, deep probabilistic programming (DPP) excels in modelling and optimising continuous probability distributions.
Hence, we introduce , a neural probabilistic logic programming language that incorporates DPP techniques into NeSy.
Doing so results in the support of inference and learning of both discrete and continuous probability distributions under logical constraints.
Our main contributions are
1) the semantics of and its corresponding inference algorithm,
2) a proven asymptotically unbiased learning algorithm, and
3) a series of experiments that illustrate the versatility of our approach.

§ INTRODUCTION

Neural-symbolic AI (NeSy) [Garcez et al., 2002, De Raedt et al., 2021] focuses on the integration of symbolic and neural methods.
The advantage of NeSy is that it combines the reasoning power of logical representations with the learning capabilities of neural networks.
Additionally, it has been shown to converge faster during learning and to be more robust [Rocktäschel and Riedel, 2017, Xu et al., 2018, Evans and Grefenstette, 2018].

The challenge of NeSy lies in combining discrete symbols with continuous and differentiable neural representations.
So far, such a combination has been realised for Boolean variables by interpreting the outputs of neural networks as the weights of these variables.

These weights can then be given either a fuzzy semantics [Badreddine et al., 2022, Diligenti et al., 2017]
or a probabilistic semantics [Manhaeve et al., 2021, Yang et al., 2020].
The latter is also used in neural probabilistic logic programming (NPLP), where neural networks parametrise probabilistic logic programs.

A shortcoming of traditional probabilistic NeSy approaches is that they fail to capture models that integrate continuous random variables and neural networks –
a feature already achieved with mixture density networks [Bishop, 1994] and more generally deep probabilistic programming (DPP) [Tran et al., 2017, Bingham et al., 2019].

However, it is unclear whether DPP can be generalised to enable logical and relational reasoning.
Hence, a gap exists between DPP and NeSy as reasoning is, after all, a fundamental component of the latter.
We contribute towards closing this DPP-NeSy gap by introducing [Sea stands for the letter C, as in continuous random variable.],
an NPLP language
with support for discrete-continuous random variables that retains logical and relational reasoning capabilities.
We achieve this integration by allowing

arbitrary and differentiable probability distributions expressed in a modern DPP language
while combining knowledge compilation [Darwiche and Marquis, 2002] with the reparametrisation trick [Ruiz et al., 2016] and continuous relaxations [Petersen et al., 2021].

Our main contributions are [label=(*)]

* the well-defined probabilistic semantics of (Section <ref>) with an inference algorithm based on

weighted model integration (WMI) [Belle et al., 2015]

(Section <ref>),

* a proven asymptotically unbiased gradient estimate for WMI that turns into a differentiable, discrete-continuous NPLP language (Section <ref>), and

* an experimental evaluation showing the versatility of discrete-continuous reasoning and the efficacy of our approach (Section <ref>)

§ LOGIC PROGRAMMING CONCEPTS

A term t is either a constant c, a variable V or a structured term of the form
where f is a functor and each t@$_i$@ is a term.
Atoms are expressions of the form
Here, q/@$K$@ is a predicate of arity $K$
and each
t@$_i$@ is a term.
A literal is an atom or the negation of an atom
$\neg \text{\probloginline{q(t@$_1$@,...,t@$_K$@)}}$.
A definite clause (also called a rule) is an expression of the form
h:- b@$_1$@,...,b@$_K$@
where h is an atom and each b@$_i$@ is a literal.
Within the context of a rule, h is called the head and the conjunction of b@$_i$@'s is referred to as the body of the rule. Rules with an empty body are called facts.
A logic program is a finite set of definite clauses.
If an expression does not contain any variables, it is called ground.
Ground expressions are obtained from non-ground ones by means of substitution.
A substitution
$\theta = \{
\text{\probloginline{V@$_1$@}}
\text{\probloginline{t@$_1$@}}
\dots,
\text{\probloginline{V@$_K$@}}
\text{\probloginline{t@$_K$@}}
\}$
is a mapping from variables V@$_i$@ to terms t<EMAIL_ADDRESS>Applying a substitution $\theta$ to an expression e (denoted e@$\theta$@)
replaces each occurrence of V@$_i$@ in e with the corresponding t<EMAIL_ADDRESS>
While pure Prolog (or definite clause logic) is defined using the concepts above, practical implementations of Prolog extend definite clause logic with  an external arithmetic engine <cit.>. Such engines enable the use of system specific routines in order to handle numeric data efficiently.
Analogous to standard terms in definite clause logic, as defined above, we introduce numeric terms. A numeric term  n@$_i$@ is either a numeric constant (a real, an integer, a float, etc.), a numeric variable N@$_i$@,
or a numerical functional term, which is an expression of the form
@$\varphi$@(n@$_1$@,...,n@$_K$@) where @$\varphi$@ is an externally defined numerical function.
The difference between a standard logical term and a numerical term is that ground numerical terms are evaluated and yield a numeric constant. For instance,  if add is a function, then
add(3, add(5, 0)) evaluates to the numerical constant 8.

Lastly, numeric constants can be compared to each other using a built-in binary comparison operator
$\bowtie\ \in
\left\{
\text{\probloginline{<}},
\text{\probloginline{=<}},
\text{\probloginline{>}},
\text{\probloginline{>=}},
\text{\probloginline{=:=}},
\text{\probloginline{=\=}}
\right\}$.
Here we use Prolog syntax to write comparison operators, which correspond to $\{ <, \leq, >, \geq, =, \neq \}$ in standard mathematical notation.
Comparison operators appear in the body of a rule, have two arguments, and are generally written as
\text{\probloginline{@$\varphi_l$@(n@$_{l,1}$@,...,n@$_{n,K}$@)}}
\bowtie
\text{\probloginline{@$\varphi_r$@(n@$_{r,1}$@,...,n@$_{r,K}$@)}}.
They evaluate
their left and right side and subsequently compare the results, assuming everything is ground.
If the stated comparison holds, the comparison is interpreted by the logic program as true, else as false.

§.§ Syntax

While facts in pure Prolog are deterministically true, in probabilistic logic programs they are annotated with the probability with which they are true.
These are the so-called probabilistic facts [De Raedt et al., 2007]. When working in discrete-continuous domains, we need to use the more general concept of distributional facts [Zuidberg Dos Martires, 2020],
inspired by the distributional clauses of Gutmann et al., 2011.

Distributional facts are expressions of the form
x   distribution(n@$_1$@,...,n@$_K$@), where x denotes a term,  the n@$_i$@'s are numerical terms and
distribution expresses the probability distribution according to which x is distributed.

[Distributional fact]
To declare a Poisson distributed variable x with rate parameter $\lambda$, one would write x   poisson(@$\lambda$@).

The meaning of a distributional fact is that all ground instances $\text{\probloginline{x}}\theta$ serve as random variables that are distributed according to $\text{\probloginline{distribution(n@$_1\theta$@,...,n@$_K\theta$@)}}$.
To obtain a neural-symbolic interface, we will allow neural networks to parametrise these distributions.

A neural distributional fact (NDF) is a distributional fact x   distribution(n@$_1$@,...,n@$_K$@) in which
a subset

of the set of numerical terms $\{ \text{\probloginline{n@$_i$@}}\}_{i=1}^{K}$ is implemented by neural networks that depend on a set of parameters $\Neuralparams$.

Random variables defined by NDFs can then be used in the logic in the form of comparisons, e.g., $\probloginline{x} > \probloginline{y}$, to reason about desired ranges of the variables.

(Probabilistic comparison formula)

A probabilistic comparison formula (PCF) is an expression of the form $(g(\boldsymbol{x}) \bowtie 0)$,
where $g$ is a function applied to the set of random variables $\variables{x}$ and $\bowtie\ \in
\left\{
\text{\probloginline{<}},
\text{\probloginline{=<}},
\text{\probloginline{>}},
\text{\probloginline{>=}},
\text{\probloginline{=:=}},
\text{\probloginline{=\=}}
\right\}$ is a binary comparison operator.
A PCF is called valid if $\left\{\boldsymbol{x}\ |\ g(\boldsymbol{x}) \bowtie 0\right\}$ is a measurable set.

[Probabilistic comparison formula]
If x is Poisson distributed and represents the number of chocolate pieces put in a chocolate biscuit, then we can use a simple PCF to define when such a biscuit passes a quality test through the rule passes_test :- (x > 11).

Note that the general form of a PCF in Definition <ref> has a $0$ on the right hand side, which can always be obtained by subtracting the right hand-side from both sides of the relation.
With the definitions of NDFs and PCFs, a program can now be formally defined.

A program consists of
a finite set of  NDFs $\dfacts$ (defining random variables),
a finite set $\compset$ of  valid PCFs
and a set of logical rules $\lrules$ that can use any of those valid PCFs in their bodies.

humid denotes a Bernoulli
random variable that takes the value $1$ with a probability $p$ given by the output of a neural network humid_detector. temp denotes a normally distributed variable whose parameters are predicted by a network temperature_predictor.
The program further contains two rules that deduce whether we have good weather or not. The first one expresses the case of snowy weather, while the second holds for a rather temperate and dry situation.
The atom query(good_weather($\inlineimage{Imagery/world.png}$)) declares that we want to compute the probability of good_weather when evaluated on the data $\inlineimage{Imagery/world.png}$.

It illustrates the neural-symbolic nature of , as its ground argument is a sub-symbolic representation ($\inlineimage{Imagery/world.png}$) of the world.

temp(Data, T)

humid(Data) =:= 1, temp(Data) < 0.

humid(Data) =:= 0, temp(Data) > 15.

Notice how the random variables humid and temp appear in the body of a logical rule with comparison operators.
In our probabilistic setting, the truth value of a comparison depends on the value of its random variables and is thus random itself.

generalises a range of existing PLP languages.

For instance, if we were to remove the distributional fact on temp and all the PCFs using them, we would obtain a program [Manhaeve et al., 2021]. If we additionally replace the neural network in humid with a fixed probability p, we end up with a probabilistic logic program [De Raedt et al., 2007]. Replacing that constant probability p by a constant 1 yields a non-probabilistic Prolog program.
Alternatively, considering all rules and facts in Example <ref> but
replacing the neural parameters of the normal distribution with numeric constants results in a Distributional Clause program [Gutmann et al., 2011].
We further discuss these connections

in Appendix <ref>, where we also formally prove that strictly generalises .

§.§ Semantics

programs are used to

compute the probability that a ground atom q is entailed.
That probability follows from the semantics of a program.
As is custom in (probabilistic) logic programming, we will define the semantics of with respect to ground programs.
We will assume that each ground distributional fact $f \in \dfacts$ defines a different random variable, as each random variable can only have one unique distribution.
Also notice that any ground neural distributional facts will contain the inputs to their neural functions.
In a sense, a program is conditioned on these neural network inputs.

To define the semantics of ground programs, we first introduce the possible worlds over the PCFs.
Every subset $\compsubset$ of a set of PCFs $\compset$ defines a possible world $\world_{\compsubset} = \{ \compsubset \cup h \theta \mid \lrules \cup \compsubset  \models h\theta \text{ and } h\theta \text{ is ground} \}$.
Intuitively speaking, the comparisons in such a subset are considered to be true and all others false. A rule with a comparison in its body that is not in this subset can hence not be used to determine the truth value of atoms. The deterministic rules $\lrules$ and the subset $\compsubset$ together define a set of all ground atoms $h\theta$ that are derivable, i.e., entailed by the program, and thus considered true. Such a set is called a possible world. We refer the reader to the paper of De Raedt and Kimmig, 2015 for a detailed account of possible worlds in a PLP context.
Following the distribution semantics of Sato, 1995 and by taking inspiration from Gutmann et al., 2011, we define the probability of a possible world.

Let be a ground program
and $\compsubset=\{ c_1,\dots, c_H \} \subseteq \compset$ a set of PCFs that depend on the random variables declared in the set of distributional facts $\dfacts$. The probability $P(\world_{\compsubset})$ of a world $\world_{\compsubset}$
is then defined as

\begin{align}
% P(\world_{\compsubset}) =
\int
\left[
\Big(\prod_{c_i\in \compsubset} \indicator(c_i) \Big)
\Big( \prod_{c_i \in \compset\setminus \compsubset}  \indicator(\bar{c}_i) \Big)
\right]
\ \differential P_{\dfacts}.
\label{eq:world_probability}
\end{align}

Here the symbol $\indicator$ denotes the indicator function, $\bar{c}_i$ is the complement of the comparison $c_i$ and $\differential P_{\dfacts}$ represents the joint probability measure of the random variables defined in the set of distributional facts $\dfacts$.

[Probability of a possible world]
Given as in Example <ref>,
where humid_detector(data1) predicts $p(\text{\probloginline{data1}})$ and temperature_predictor(data1) predicts the tuple $(\mu(\text{\probloginline{data1}}),\sigma(\text{\probloginline{data1}}))$, the probability of the possible world $\omega_{\left\{\text{\probloginline{temp(data1)>15}},\ \text{\probloginline{humid(data1)=:=1}}\right\}}$ is given by

\begin{align}
\int\indicator(x {>} 15)
\frac{
\exp
\left(
\left(x-\mu(\text{\probloginline{data1}}) \right)^2
2 \sigma^2(\text{\probloginline{data1}})
\right)
}{\sqrt{2\pi} \sigma(\text{\probloginline{data1}})}
\ \differential x.
\end{align}

Indeed, the measure $\differential P_{\dfacts}$ decomposes into a counting measure and the product of a Gaussian density function with a differential. The counting measure leads to the factor $p(\text{\probloginline{data1}})$, since that is the probability that humid(data1)=:=1. Hence, the products in Equation <ref> reduce to a single indicator of the PCF $(x > 15)$.

The probability of a ground atom $q$ is given by

\begin{align}
P(q) = \sum_{\compsubset \subseteq \compset : q \in \world_{\compset} } P(\world_{\compsubset})
\label{eq:query_probability}.
\end{align}

[Measureability of query atom]propositionquerymeasurability

Let be a program, then
defines, for an arbitrary query atom $q$, the probability that $q$ is true.

See Appendix <ref>.

§ INFERENCE AND LEARNING

§.§ Inference via weighted logic

A popular technique to perform inference in probabilistic logic programming uses a reduction to so-called weighted model counting (WMC); instead of computing the probability of a query, one computes the weight of a propositional logical formula [Chavira and Darwiche, 2008, Fierens et al., 2015]. For , the equivalent approach is to map a ground program onto a satisfiability modulo theory (SMT) formula [Barrett and Tinelli, 2018]. The analogous concept to WMC for these formulas is weighted model integration (WMI) [Belle et al., 2015], which can handle infinite sample spaces.
In all that follows, for ease of exposition, we assume that all joint probability distributions are continuous.

[Inference as WMI]propositioninferenceaswmi

Assume that the measure $\differential P_{\dfacts}$ decomposes into a joint probability density function $\weight(\variables{x})$ and a differential $\differential \variables{x}$, then the probability $P(q)$ of a query atom $q$ can be expressed as the weighted model integration problem

\begin{align}
% P(q) =
\int \left[
\sum_{\compsubset \subseteq \compset : q \in \world_{\compset}} \prod_{c_i\in \compsubset {\cup} \negcompsubset}  \indicator(c_i(\variables{x}))
\right]
\weight(\variables{x})
\ \mathrm{d}\variables{x},
\label{eq:probability_query_as_wmi}
\end{align}

\negcompsubset \coloneqq \left\{\bar{c}_i\ |\ c_i \in \compset {\setminus} \compsubset\right\}

See Appendix <ref>.

Being able to express the probability of a queried atom in as a weighted model integral allows us to adapt and deploy inference techniques developed in the weighted model integration literature for .

We opt for the approximate inference algorithm Sampo presented in Zuidberg Dos Martires et al., 2019 because of its more scalable nature.

Sampo uses knowledge compilation [Darwiche and Marquis, 2002], a state-of-the-art technique for probabilistic logic inference [Chavira and Darwiche, 2008, Fierens et al., 2015].

Intuitively, knowledge compilation is a two-step procedure applied to a logical formula with PCFs, i.e., an SMT formula.
First, it infers the exact probability of all PCFs containing discrete variables through symbolic inference.
Then, it converts the remainder of the SMT formula into a polynomial in terms of those exact probabilities and the PCFs containing continuous random variables (Figure <ref>).
This polynomial is the integrand of Equation <ref>. All that remains is to approximate the integration of this polynomial by sampling from the joint probability distribution $\weight(\variables{x})$ of the continuous random variables.

In other words, Sampo computes the expression

\begin{align}
P(q) =
\int \amc(\variables{x}) \cdot  \weight(\variables{x}) \, \differential\variables{x}
% \ \approx\
\approx
\frac{1}{|\mathcal{X}|}
\sum_{\variables{x} \in \mathcal{X}} \amc(\variables{x}),
\label{eq:probability_query_approx}
\end{align}

where $\mathcal{X}$ denotes a set of samples drawn from $\weight(\variables{x})$ and $\amc(\variables{x})$ is the result of knowledge compilation, i.e., the sum of products of indicator functions in Equation <ref>.

Diagrammatic representation of the result of knowledge compilation for the query in Example <ref>.
The blue boxes originate from PCFs over discrete variables, while the orange ones are PCFs over continuous variables. Note how the discrete variable PCFs are reduced to their exact probabilities while the continuous PCFs still need to be inferred.

We stress that the Sampo algorithm only samples random variables whose expected value with respect to the function $\amc(\variables{x})$ can not be computed exactly.
Hence, in the absence of continuous random variables, our implementation of using Sampo coincides with on both a semantics level (Proposition <ref>) and inference level.

§.§ Learning via differentiation

A program depend on a set of (neural) parameters
$\Neuralparams$ (Definition <ref>).
In order to optimise these parameters, we need to take their gradients of a loss function that compares the probability $P(q)$ to a training signal. More precisely, we need to compute the derivative

\begin{align}
\deriv \mathcal{L}(P_{\Neuralparams}(q)) =
%    \ \approx\
\partial_{P_{\Neuralparams}(q)} \mathcal{L}(P_{\Neuralparams}(q))\cdot \deriv P_{\Neuralparams}(q),
\end{align}

where we explicitly indicate the dependency of the probability on $\Neuralparams$ and $\neuralparam\in \Neuralparams$.
Differentiating $P_{\Neuralparams}(q)$ with respect to $\lambda$ presents two obstacles.
the question of differentiating through the sampling  process of Equation <ref>
and second,
the non-differentiability of the indicator functions in $\amc(\variables{x})$.

The non-differentiability of sampling is tackled using the reparametrisation trick [Ruiz et al., 2016].
Reparametrisation offers better estimates than other approaches, such as REINFORCE [Williams, 1992] and is readily utilised in modern probabilistic programming languages such as Tensorflow Probability [Tran et al., 2017] and Pyro [Bingham et al., 2019].
Conversely, the non-differentiability of the indicator functions prevents swapping the order of differentiation and integration [Flanders, 1973], which we resolve by applying continuous relaxations following the work of Petersen et al., 2021.
Together, we obtain the gradient estimate

\begin{align}
\deriv P_{\Neuralparams}(q)
\deriv
\int \amc(\variables{x}) \cdot  \weight_{\Neuralparams}(\variables{x}) \, \differential\variables{x} \\
% &\approx\
% \deriv \int \amc_\softened(\reparam) \cdot
% p(\variables{u})
% \ \differential \variables{u} \\
\int
\left[
\deriv \amc_\softened(\reparam)
\right] \cdot
\ \differential \variables{u}
\label{eq:probability_derivative}
\end{align}

where the subscript $s$ in $\amc_\softened(\variables{x}) $ denotes the continuously relaxed or softened version of $\amc(\variables{x})$ and $r(\variables{u}, \Neuralparams)$ is the reparametrisation function.

Our gradient estimate using relaxations is asymptotically unbiased.
As an example of these relaxations, consider the indicator of a PCF $(g(\boldsymbol{x}) > 0)$, which is relaxed into the sigmoid $\sigma(\coolness\cdot g(\boldsymbol{x}))$.

Appendix <ref> provides more details on relaxations of general PCFs.
The coolness parameter $\coolness\in \mathbb{R}_+^0$

determines the strictness of the relaxation. Hence, we recover the hard indicator function when $\coolness{\rightarrow} +\infty$.
Note that relaxing indicator functions
introduces bias.

Petersen et al., 2021 already stated in their work that, in the infinite coolness limit, a relaxed function coincides with the non-relaxed one. Proposition <ref> extends this result

to the derivatives of relaxed and non-relaxed functions, proving that our gradient estimate is asymptotically unbiased.

[Unbiased in the infinite coolness limit]propositionunbiasedapprox

Let $\mathbb{P}$ be a program with PCFs $(g_i(\boldsymbol{x}) \bowtie 0)$ and corresponding coolness parameters $\coolness_i$.

If all $\deriv (g_i \circ r)$ are locally integrable over $\mathbb{R}^k$ and
every $\coolness_i \rightarrow +\infty$,
then we have, for any query atom $q$, that

\begin{align}
\deriv P(q)
\int
\deriv\amc_\softened(\reparam)
\cdot p(\variables{u})
\ \differential \variables{u}
\end{align}

The proof makes use of the mathematical theory of distributions [Schwartz, 1957], which generalise the concept of functions, and is given in Appendix <ref>.

Finally, we obtain a practical and unbiased estimate of $\deriv P_{\Neuralparams}(q)$ using a set of samples $\mathcal{U}$ drawn from

\begin{align}
\deriv P(q)
\int
\left[
\deriv      \amc_\softened(\reparam)
\right]
\cdot p(\variables{u})
\ \differential \variables{u} \\
\frac{1}{|\mathcal{U}|} \sum_{\variables{u} \in \mathcal{U}} \deriv \amc_\softened(\reparam).
\label{eq:WMIDerFinalApprox}
\end{align}

Computing this gradient estimate does not require drawing new samples. Implementing the relaxations of PCFs in a straight-through manner allows us to directly apply automatic differentiation on the inferred probability.

§.§ Probabilistic programming connections

Since knowledge compilation symbolically infers discrete random variables, we only have to sample from a continuous joint probability distribution.
To sample such distributions, we can fully exploit the advanced inference and learning techniques [Hoffman et al., 2014] of modern probabilistic programming languages [Tran et al., 2017, Bingham et al., 2019]. Our implementation of utilises Tensorflow Probability for this task, effectively using knowledge compilation as a differentiable bridge between logical and probabilistic reasoning.
While this bridge is limited to sampling techniques for now, it presents an interesting direction for future work to completely unify NeSy with DPP.

§.§ Limitations

While the use of relaxations is well-known and used in recent gradient estimators [Tucker et al., 2017, Grathwohl et al., 2018], the bias they introduce is often hard to deal with in practice.
In our case, this bias only reduces to zero in the infinite coolness limit (Proposition <ref>), meaning the use of annealing can be necessary. Finding a good annealing scheme for any problem is non-trivial and effectively introduces another component in need of optimisation. However, as relaxations allow the use of the reparametrisation trick, the resulting lower variance estimates together with our theoretical guarantees support our choice. A more detailed discussion of the current limitations of can be found in Appendix <ref>.

§ RELATED WORK

From a NeSy perspective the formalism most closely related to is that of Logic Tensor Networks (LTNs) [Badreddine et al., 2022]. The main difference between LTNs and is the fuzzy logic semantics of the former and the probabilistic semantics of the latter.
Interestingly, LTNs and other NeSy approaches based on fuzzy logic also require relaxations to incorporate continuous values.

However, fuzzy-based approaches require these relaxations at the semantics level, in contrast to .
Even more, they can only compare continuous point values instead of more general continuous random variables.
LTNs' fuzzy semantics also exhibit drawbacks on a more practical level.
Unlike with its probabilistic semantics, LTNs are not inherently capable of neural-symbolic generative modelling (Section <ref>).
For a broader overview of the field of neural-symbolic AI, we refer the reader to a series of survey papers that have been published in recent years [Garcez et al., 2019, Marra et al., 2021, Garcez et al., 2022].

From a probabilistic programming perspective, is related to languages that handle discrete and continuous random variables such as BLOG [Milch, 2006], Distributional Clauses [Gutmann et al., 2011] and Anglican [Tolpin et al., 2016], which have all been given declarative semantics, i.e., the meaning of the program does not depend on the underlying inference algorithm. However, these languages have the drawback of non-differentiability.
This drawback stands in stark contrast to end-to-end (deep) probabilistic programming languages such as Pyro [Bingham et al., 2019] or Tensorflow Probability [Dillon et al., 2017], but these have only been equipped with operational semantics and do not support logical constraints.
not only introduces the ability to express such logical constraints in the form of PCFs to construct challenging posterior distributions, but

does so in an end-to-end differentiable fashion.

Finally, our gradient estimate can be related to relaxation-based methods like REBAR [Tucker et al., 2017] or RELAX [Grathwohl et al., 2018], but without the REINFORCE [Williams, 1992] inspired component. Instead, we utilise the differentiability of knowledge compilation to obtain exact gradients of discrete variables.
Since our inference scheme innately requires knowledge compilation, the use of other discrete gradient estimators like [Niepert et al., 2021] does not directly apply to .
Moreover, we exploit the structure of our problem by directly relaxing comparison formulae in a sound manner [Petersen et al., 2021], in contrast to introducing an artificial relaxation of the whole problem [Grathwohl et al., 2018].

§ EXPERIMENTAL EVALUATION

We illustrate the versatility of by tackling three different problems. Section <ref> discusses the detection of handwritten dates without location supervision. In Section <ref> a hybrid Bayesian network with conditional probabilities dependent on the satisfaction of certain logical constraints will be optimised. Finally, Section <ref> introduces neural-symbolic variational auto-encoders, inspired by Misino et al., 2022.

The details of our experimental setup, including the precise programs, coolness annealing schemes and hyperparameters used for the neural networks are given in Appendix <ref>.

§.§ Neural-symbolic attention

A problem that cannot yet be solved to a satisfactory degree by purely neural or other neural-symbolic systems is detecting handwritten years.
Given a single image with a handwritten year, the task is to predict the correct year as a sequence of 4 digits together with the location of these digits (Figure <ref>, left).
This year can be anywhere in the image and the only supervision is in the digits of the year, not where these digits are in the image. In other words, the problem is equivalent to object detection without bounding box supervision.

Solving such a problem seems to be out of scope for current methods.
On the one hand, existing neural approaches are often complex pipelines of neural components that break end-to-end differentiability [Seker and Ahn, 2022].
On the other hand, current neural-symbolic methods lack sufficient spatial reasoning capabilities in order to perform the necessary image segmentation.

We exploit probabilistic programming by modelling the location of a digit as a deep generalised normal distribution [Nadarajah, 2005].
That is, we use a convolutional neural network to regress the parameters of four generalised normal distributions, one for each digit of a year.
Then, we take inspiration from the spatial transformer literature [Carion et al., 2020] and convert the distribution of each location to an attention map (Figure <ref>, right).

On the left, an example of a handwritten year.

On the right, the attention map for the digit 8 as a generalised normal distribution.

Intuitively, we can view generalised normal distributions as differentiable bounding boxes. This allows gradients to flow from a downstream classification network to the regression component.

In our experimental validation we compare to a neural baseline and logic tensor networks.
The neural baseline applies the four probabilistic attention maps, one for the location of each of the four digits, to the input image.
The resulting four attenuated images are then passed on to a classification network without additional reasoning.
The network is simply required to predict the digits in the right order.

With , we encode that a year is a sequence of digits, i.e., the order matters, by enforcing an explicit order on the digit locations.
Doing so requires spatial reasoning, i.e., reasoning which digit is at which location.

For LTNs, we encode the same information.
However, as LTNs lack a proper distribution semantics, they can only reason on the level of the expected values of the generalised normal distributions.

In our experiment, the sets of years appearing in the training, validation and test data are all disjoint.
Moreover, the sets of handwritten digits used to generate those years are also disjoint.
Partitioning the data in such a way leads to a challenging learning problem;
the difficulty lies in out-of-distribution inference, as the years and handwritten digits in the validation and test set have never been seen during training.

We evaluate all methods in terms of accuracy and Intersection-over-Union (IoU).
For the accuracy, we compare the sequence of predicted digits to the correct sequence of digits constituting a year.
A prediction is correct if all digits are correctly predicted in the right order.
For the IoU, we map each predicted generalised normal distribution to a bounding box by using the mean as the centre and the scale parameter as the width of the box. The IoU is then given by the overlap between this box and the true location of the handwritten digit.

Mean accuracy and IoU with standard error for classifying the correct year, taken over 10 runs.

Method     2cResults

acc.     IoU

2 - 3
$93.77\pm 0.57\phantom{0}$     $17.69\pm 0.23$

LTN      $76.50\pm 12.10$     $10.73\pm 1.69$

Neural Baseline      $54.71\pm 14.33$     0$6.26\pm 1.77$

We present our results in Table <ref>.
The most striking observation is the poor performance and large variance of the neural baseline.
It fails to predict the location of the digits in the right order, as can be seen from the lower IoU values.
Since classification depends on the predicted locations, these lower values also explain the lack in accuracy.
We can conclude that the neural baseline struggles to generalise to out-of-distribution data.
While LTNs fare better, the high standard error on the accuracy indicates that their continuous reasoning capabilities are insufficient to guarantee consistent solutions.
distinguishes itself by a higher and more consistent accuracy.
The reason is also clear; exploits the entire domain of the distribution of each location. This then leads to a higher IoU value that in turn results in a higher accuracy.

§.§ Neural hybrid Bayesian networks

Hybrid Bayesian networks [Lerner, 2003] are probabilistic graphical models that combine discrete and continuous random variables.
allows for the introduction of optimisable neural components and logical constraints to such models, as shown in Example <ref>.

We further extend this example (Figure <ref>) and specify the datasets that form the input to the various neural networks.
The temperature is predicted from a real meteorological dataset [Cho et al., 2020] and we use CIFAR-10 images as proxies for observing clouds and humidity.

Moreover, dependencies on a number of constraints are added, which goes beyond the capabilities of traditional probabilistic programming.

Graphical model of Enjoying the weather (E).(E) holds when Depressed (D) is not true and there is Good weather (G).
A person has a higher probability of being depressed when it is Cloudy (C), while the degree of good weather is beta distributed depending on various logical constraints on Temperature (T) and Rain (R).

Finally, rain is probable when it is both Cloudy and Humid (H).

Our neural Bayesian model was optimised by only giving probabilistic supervision on whether E was true or false, i.e., the weather was enjoyed or not.
Given our model, such distant supervision only translates into a learning signal on different ranges of temperature values that satisfy different PCFs. We will see that 's reasoning over the full domain of the temperature distribution allows it to perform meaningful density estimation from such a signal.

The optimised Bayesian model can be evaluated in two ways. First, the accuracy on CIFAR-10 of the networks utilised in Cloudy and Humid, which were

$95.24\pm 3.32$ and $98.96\pm 0.11$, respectively.
Second, we measure the quality of the density estimation on Temperature by looking at the MSE between the true and predicted mean values, which was

$0.1799 \pm 0.0139$

Importantly, was able to approximate the standard deviation of Temperature from just the distant supervision, deviating by only

$0.60\pm 0.22$.

§.§ Neural-symbolic variational auto-encoder

Probabilistic programming is well-suited to generative tasks, but it can not perform generation conditioned on logical constraints.
Inspired by the work of Misino et al., 2022, we showcase how extends the generative power of probabilistic programming to such constraints.
To this end, we will consider the task of learning to generate 2 images of digits given the value of their subtraction.

A diagrammatic overview of our program is given in Figure <ref>.
It uses a conditional variational auto-encoder (CVAE) [Sohn et al., 2015] to generate images conditioned on a digit value.
finds those digit values from a given subtraction result by logical reasoning.
It can also condition generation on other variables in the CVAE latent space as this space is an integral part of 's deep, relational model.

We will exploit this property later on when we extend the task to generating digits in the same writing style as a given image without any additional optimisation.

Given example pairs of images and the value of their subtraction, e.g., $(\inlineimage{Imagery/2_6_6_7_original.png}, \inlineimage{Imagery/5_3_29_2_original.png})$ and $3$, the CVAE encoder vae_latent first encodes each image into a multivariate normal NDF (latent) and a latent vector. The latter is the input of a categorical NDF digit, completing the CVAE latent space. Supervision is dual; generated images are compared to the original ones in a probabilistic reconstruction loss, while both digits need to subtract to the given value.

Both the CVAE and digit classifier are successfully trained jointly. Example generations of images that satisfy the subtraction result
\inlineimage{Imagery/generation.pdf}
\inlineimage{Imagery/generation.pdf}
= 5
can be seen below. In general, finds all possible digits that subtract to a given value and generates images for each correct combination. Below, we left out 2 such combinations for clarity of exposition.

While our program is inspired by the VAEL architecture of Misino et al., 2022, conceptual differences exist.
Most notably, for VAEL, the image generation resides outside the probabilistic logic program. Conversely, the CVAE, including its latent space, is explicitly declared and accessible in .
This difference allows to generalise to conditional generative queries that differ significantly from the original optimisation task.
For example, we can zero-shot query the program to fill in the blank in

\inlineimage{Imagery/imageconditionalgeneration_input2_grey.png}
\inlineimage{Imagery/generation.pdf}
= \text{\probloginline{Diff}}$
instead of the two blanks of the learning task
\inlineimage{Imagery/generation.pdf}
\inlineimage{Imagery/generation.pdf}
= \text{\probloginline{Diff}}$.
Even more, we can enforce that the generated digit is in the same writing style as the given digit by conditioning the generation on the latent space of the given image (Figure <ref>).

Four random images of right digits (top row) and their generated left digits for 3 given random difference values (bottom row). Note the preservation of the style of the given minuends.

§ CONCLUSION

We presented , a novel neural-symbolic probabilistic logic programming language that integrates hybrid probabilistic logic and neural networks. Inference is dealt with efficiently through approximate weighted model integration while learning is facilitated by reparametrisation and continuous relaxations of non-differentiable logic components. Our experiments illustrate how is capable of intricate probabilistic modelling allowing for meaningful weak supervision while maintaining strong out-of-distribution performance. Moreover, they show how hybrid probabilistic logic can be used as a flexible structuring formalism for the neural paradigm that can effectively optimise and reuse neural components in different tasks.

[Abadi, 2016]
Martín Abadi.
Tensorflow: learning functions at scale.
In International Conference on Functional Programming, 2016.

[Badreddine et al., 2022]
Samy Badreddine, Artur d'Avila Garcez, Luciano Serafini, and Michael Spranger.
Logic tensor networks.
Artificial Intelligence, 2022.

[Barrett and Tinelli, 2018]
Clark Barrett and Cesare Tinelli.
Satisfiability modulo theories.
In Handbook of model checking. Springer, 2018.

[Belle et al., 2015]
Vaishak Belle, Andrea Passerini, and Guy Van den Broeck.
Probabilistic inference in hybrid domains by weighted model
In IJCAI, 2015.

[Bingham et al., 2019]
Eli Bingham, Jonathan P Chen, Martin Jankowiak, Fritz Obermeyer, Neeraj
Pradhan, Theofanis Karaletsos, Rohit Singh, Paul Szerlip, Paul Horsfall, and
Noah D Goodman.
Pyro: Deep universal probabilistic programming.
The Journal of Machine Learning Research, 2019.

[Bishop, 1994]
Christopher M Bishop.
Mixture density networks.
Technical report, Aston University, 1994.

[Carion et al., 2020]
Nicolas Carion, Francisco Massa, Gabriel Synnaeve, Nicolas Usunier, Alexander
Kirillov, and Sergey Zagoruyko.
End-to-end object detection with transformers.
In Computer Vision–ECCV. Springer, 2020.

[Chavira and Darwiche, 2008]
Mark Chavira and Adnan Darwiche.
On probabilistic inference by weighted model counting.
Artificial Intelligence, 2008.

[Cho et al., 2020]
Dongjin Cho, Cheolhee Yoo, Jungho Im, and Dong-Hyun Cha.
Comparative assessment of various machine learning-based bias
correction methods for numerical weather prediction model forecasts of
extreme air temperatures in urban areas.
Earth and Space Science, 2020.

[Darwiche and Marquis, 2002]
Adnan Darwiche and Pierre Marquis.
A knowledge compilation map.
Journal of Artificial Intelligence Research, 2002.

[De Raedt and Kimmig, 2015]
Luc De Raedt and Angelika Kimmig.
Probabilistic (logic) programming concepts.
Machine Learning, 2015.

[De Raedt et al., 2007]
Luc De Raedt, Angelika Kimmig, and Hannu Toivonen.
Problog: A probabilistic prolog and its application in link
In IJCAI. Hyderabad, 2007.

[De Raedt et al., 2021]
Luc De Raedt, Sebastijan Dumančić, Robin Manhaeve, and Giuseppe
From statistical relational to neural-symbolic artificial
In IJCAI, 2021.

[Diligenti et al., 2017]
Michelangelo Diligenti, Marco Gori, and Claudio Saccà.
Semantic-based regularization for learning and inference.
Artif. Intell., 2017.

[Dillon et al., 2017]
Joshua V Dillon, Ian Langmore, Dustin Tran, Eugene Brevdo, Srinivas Vasudevan,
Dave Moore, Brian Patton, Alex Alemi, Matt Hoffman, and Rif A Saurous.
Tensorflow distributions.
arXiv, 2017.

[Evans and Grefenstette, 2018]
Richard Evans and Edward Grefenstette.
Learning explanatory rules from noisy data.
Journal of Artificial Intelligence Research, 2018.

[Fierens et al., 2015]
Daan Fierens, Guy Van den Broeck, Joris Renkens, Dimitar Shterionov, Bernd
Gutmann, Ingo Thon, Gerda Janssens, and Luc De Raedt.
Inference and learning in probabilistic logic programs using weighted
boolean formulas.
Theory and Practice of Logic Programming, 2015.

[Flanders, 1973]
Harley Flanders.
Differentiation under the integral sign.
The American Mathematical Monthly, 1973.

[Garcez et al., 2019]
Artur d'Avila Garcez, Marco Gori, Luis C Lamb, Luciano Serafini, Michael
Spranger, and Son N Tran.
Neural-symbolic computing: An effective methodology for principled
integration of machine learning and reasoning.
arXiv, 2019.

[Garcez et al., 2022]
Artur d’Avila Garcez, Sebastian Bader, Howard Bowman, Luis C Lamb, Leo
de Penning, BV Illuminoo, Hoifung Poon, and Gerson Zaverucha.
Neural-symbolic learning and reasoning: A survey and interpretation.
Neuro-Symbolic Artificial Intelligence: The State of the Art,

[Garcez et al., 2002]
Artur S d'Avila Garcez, Krysia Broda, Dov M Gabbay, et al.
Neural-symbolic learning systems: foundations and
Springer Science & Business Media, 2002.

[Grathwohl et al., 2018]
Will Grathwohl, Dami Choi, Yuhuai Wu, Geoff Roeder, and David Duvenaud.
Backpropagation through the void: Optimizing control variates for
black-box gradient estimation.
In ICLR, 2018.

[Griewank and Walther, 2008]
Andreas Griewank and Andrea Walther.
Evaluating derivatives: principles and techniques of
algorithmic differentiation.
SIAM, 2008.

[Gutmann et al., 2011]
Bernd Gutmann, Ingo Thon, Angelika Kimmig, Maurice Bruynooghe, and Luc
De Raedt.
The magic of logical inference in probabilistic programming.
Theory and Practice of Logic Programming, 2011.

[Hoffman et al., 2014]
Matthew D Hoffman, Andrew Gelman, et al.
The no-u-turn sampler: adaptively setting path lengths in hamiltonian
monte carlo.
J. Mach. Learn. Res., 2014.

[Huang et al., 2021]
Jiani Huang, Ziyang Li, Binghong Chen, Karan Samel, Mayur Naik, Le Song, and
Xujie Si.
Scallop: From probabilistic deductive databases to scalable
differentiable reasoning.
NeurIPS, 2021.

[Jang et al., 2017]
Eric Jang, Shixiang Gu, and Ben Poole.
Categorical reparametrization with gumble-softmax.
In ICLR, 2017.

[Kingma and Ba, 2015]
Diederik P Kingma and Jimmy Ba.
Adam: A method for stochastic optimization.
In ICLR, 2015.

[Lerner, 2003]
Uri Nahum Lerner.
Hybrid Bayesian networks for reasoning about complex systems.
stanford university, 2003.

[Maddison et al., 2017]
Chris J Maddison, Andriy Mnih, and Yee Whye Teh.
The concrete distribution: A continuous relaxation of discrete random
In ICLR, 2017.

[Manhaeve et al., 2021]
Robin Manhaeve, Sebastijan Dumančić, Angelika Kimmig, Thomas
Demeester, and Luc De Raedt.
Neural probabilistic logic programming in deepproblog.
Artificial Intelligence, 2021a.

[Manhaeve et al., 2021]
Robin Manhaeve, Giuseppe Marra, and Luc De Raedt.
Approximate inference for neural probabilistic logic programming.
In KR, 2021b.

[Marra et al., 2021]
Giuseppe Marra, Sebastijan Dumančić, Robin Manhaeve, and Luc
De Raedt.
From statistical relational to neural symbolic artificial
intelligence: a survey.
arXiv, 2021.

[Milch, 2006]
Brian Christopher Milch.
Probabilistic models with unknown objects.
PhD thesis, University of California, Berkeley, 2006.

[Misino et al., 2022]
Eleonora Misino, Giuseppe Marra, and Emanuele Sansone.
Vael: Bridging variational autoencoders and probabilistic logic
In NeurIPS, 2022.

[Morettin et al., 2021]
Paolo Morettin, Pedro Zuidberg Dos Martires, Samuel Kolb, and Andrea Passerini.
Hybrid probabilistic inference with logical and algebraic
constraints: a survey.
In IJCAI, 2021.

[Nadarajah, 2005]
Saralees Nadarajah.
A generalized normal distribution.
Journal of Applied statistics, 2005.

[Niepert et al., 2021]
Mathias Niepert, Pasquale Minervini, and Luca Franceschi.
Implicit mle: backpropagating through discrete exponential family
NeurIPS, 2021.

[Paszke et al., 2019]
Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory
Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, et al.
Pytorch: An imperative style, high-performance deep learning library.
NeurIPS, 2019.

[Petersen et al., 2021]
Felix Petersen, Christian Borgelt, Hilde Kuehne, and Oliver Deussen.
Learning with algorithmic supervision via continuous relaxations.
NeurIPS, 2021.

[Riguzzi and Swift, 2013]
Fabrizio Riguzzi and Terrance Swift.
Well–definedness and efficient inference for probabilistic logic
programming under the distribution semantics.
Theory and practice of logic programming, 2013.

[Rocktäschel and Riedel, 2017]
Tim Rocktäschel and Sebastian Riedel.
End-to-end differentiable proving.
NeurIPS, 2017.

[Ruiz et al., 2016]
Francisco J. R. Ruiz, Michalis K. Titsias, and David M. Blei.
The generalized reparameterization gradient.
NeurIPS, 2016.

[Sato, 1995]
Taisuke Sato.
A statistical learning method for logic programs with distribution
In ICLP, 1995.

[Schwartz, 1957]
Laurent Schwartz.
Théorie des distributions à valeurs vectorielles. i.
In Annales de l'institut Fourier, 1957.

[Seker and Ahn, 2022]
Ahmet Cagatay Seker and Sang Chul Ahn.
A generalized framework for recognition of expiration dates on
product packages using fully convolutional networks.
Expert Systems with Applications, 2022.

[Sohn et al., 2015]
Kihyuk Sohn, Honglak Lee, and Xinchen Yan.
Learning structured output representation using deep conditional
generative models.
NeurIPS, 2015.

[Sterling and Shapiro, 1994]
Leon Sterling and Ehud Y Shapiro.
The art of Prolog: advanced programming techniques.
MIT press, 1994.

[Tolpin et al., 2016]
David Tolpin, Jan-Willem van de Meent, Hongseok Yang, and Frank Wood.
Design and implementation of probabilistic programming language
In Symposium on the Implementation and Application of
Functional programming Languages, 2016.

[Tran et al., 2017]
Dustin Tran, Matthew D. Hoffman, Rif A. Saurous, Eugene Brevdo, Kevin Murphy,
and David M. Blei.
Deep probabilistic programming.
In ICLR, 2017.

[Tucker et al., 2017]
George Tucker, Andriy Mnih, Chris J Maddison, John Lawson, and Jascha
Rebar: Low-variance, unbiased gradient estimates for discrete latent
variable models.
NeurIPS, 2017.

[Williams, 1992]
Ronald J Williams.
Simple statistical gradient-following algorithms for connectionist
reinforcement learning.
Machine learning, 1992.

[Winters et al., 2022]
Thomas Winters, Giuseppe Marra, Robin Manhaeve, and Luc De Raedt.
Deepstochlog: Neural stochastic logic programming.
In AAAI, 2022.

[Xu et al., 2018]
Jingyi Xu, Zilu Zhang, Tal Friedman, Yitao Liang, and Guy Van den Broeck.
A semantic loss function for deep learning with symbolic knowledge.
In ICML, 2018.

[Yang et al., 2020]
Zhun Yang, Adam Ishay, and Joohyung Lee.
Neurasp: Embracing neural networks into answer set programming.
In IJCAI, 2020.

[Zuidberg Dos Martires, 2020]
Pedro Zuidberg Dos Martires.
From Atoms to Possible Worlds: Probabilistic Inference in the
Discrete-Continuous Domain.
PhD thesis, KU Leuven, 2020.

[Zuidberg Dos Martires et al., 2019]
Pedro Zuidberg Dos Martires, Anton Dries, and Luc De Raedt.
Exact and approximate weighted model integration with probability
density functions using knowledge compilation.
In AAAI, 2019.

§ SPECIAL CASES OF

The syntax and semantics of generalise a number of probabilistic logic programming dialects. For instance, if we assume no dependency of the distributional facts on input data or external neural functions, we obtain a language equivalent to [Gutmann et al., 2011]'s  Distributional Clauses (DC) [Gutmann et al., 2011] when restricted to distributional facts. Finally, if we allow for data dependent neural functions in the NDFs but restrict them to Bernoulli and categorical distributions, we obtain [Manhaeve et al., 2021]'s  [Manhaeve et al., 2021] as a special case.
[strictly generalises ]propositiondeepproblogspecial

is a strict subset of where the set of comparison predicates is restricted to $\{ \text{\probloginline{=:=}} \}$, comparisons involve exactly one random variable and the measure $\differential P_{\dfacts}$ factorizes as a product of independent Bernoulli measures
\prod_{i: \text{\probloginline{x@$_i$@~b@$_i$@}}\in \dfacts}
\differential P_{b_i}
The subscript on $\differential P_{b_i}$ explicitly identifies the measure as the $i^{\text{th}}$ Bernoulli measure and the indices of the product go over all the (Bernoulli) random variables defined in the set of distributional facts $\dfacts$.

We prove Proposition <ref> by showing that applying the restrictions on the constraints and measure in a program leads to possible worlds that have the same probability of being true as in .
First we write down the definition of the probability $P(\world_{\compsubset})$ of a possible world in a program

\begin{align}
\int
\left[
\left(\prod_{c_i\in \compsubset} \indicator(c_i) \right)
\left( \prod_{c_i \in \compset\setminus \compsubset}  \indicator(\bar{c}_i) \right)
\right]
\ \differential P_{\dfacts}.
\label{eq:proof_deepproblogspecial_1}
\end{align}

Now observe that, since there are only Bernoulli distributions, we only need to consider two possible outcomes of a random variable x@$_i$@, either zero or one. Therefore, only two kinds of comparisons are present in the program, x@$_i$@=:=0 or x@$_i$@=:=1 (remember that we restrict ourselves to univariate comparisons). Now note that the following equivalence
\text{\probloginline{x@$_i$@=:=1}}\leftrightarrow \neg (\text{\probloginline{x@$_i$@=:=0}})
holds, which means that we can arbitrarily limit comparisons to one of the two possible outcomes of a random variable, e.g., x@$_i$@=:=0.

This equivalence can be used to replace the constraints $c_i$ in Equation <ref> by equality constraints involving comparisons to the zero outcome, i.e., $P(\world_{\compsubset})$ is equal to

\begin{align}
\label{eq:integralofproduct}
\int
% \left[
&\left(\prod_{i:c_i\in \compsubset} \indicator(x_i{=}0) \right)
\cdot \\
&\left( \prod_{i:c_i \in \compset\setminus \compsubset}  \indicator(x_i{\neq}0) \right)
% \right]
\prod_{i: \text{\probloginline{x@$_i$@~b@$_i$@}}\in \dfacts}
\ \differential P_{b_i},\nonumber
\end{align}

where the factorisation of the measure was also applied. Next, we introduce the following notation for the random variables present in the set of constraints $\compsubset$ and $\compset\setminus \compsubset$:

\begin{align}
\variables{x}^+ &\coloneqq {x_i: c_i
\in \compsubset}
\\
\variables{x}^- &\coloneqq {x_i: c_i
\in \compset\setminus \compsubset}
\end{align}

Note that we only need to consider the case where $\variables{x}^+ \cap \variables{x}^-=\emptyset$, as otherwise the probability of the possible world would simply be zero and would not contribute to the overall probability of the query atom. Because of this, we can further factorize the measure as

\begin{align}
&\prod_{i: \text{\probloginline{x@$_i$@~b@$_i$@}}\in \dfacts}
\ \differential P_{b_i}
\\
\underbrace{
\left(
\prod_{i:x_i\in \variables{x}^+}
\ \differential P_{b_i}
\right)
\eqqcolon \differential P^+
\underbrace{
\left(
\prod_{i:x_i \in \variables{x}^-} \ \differential P_{b_i}
\right)
\eqqcolon \differential P^-
\end{align}

so the integral of a product in Equation <ref> can be rewritten as the product of integrals

\begin{align}
= \\
\int
\left(
\prod_{i:c_i\in \compsubset}
\indicator(x_i{=}0)
\ \differential P^+
\right)
\right]
\cdot
\nonumber
\\
\int
\left(
\prod_{i:c_i \in \compset\setminus \compsubset}
\indicator(x_i{\neq}0)
\ \differential P^-
\right)
\right].
\nonumber
\end{align}

We have two integrals with integrands that are a product of univariate comparisons. In other words, the factors are all independent. Furthermore, we have a Bernoulli product measure, which means that we can again push the integral inside the product to yield

\begin{align}
= \\
\prod_{i:c_i\in \compsubset}
\left(
\int
\indicator(x_i{=}0)
\ \differential P^+
\right)
\right]
\cdot
\nonumber
\\
\prod_{i:c_i \in \compset\setminus \compsubset}
\left(
\int
\indicator(x_i{\neq}0)
\ \differential P^-
\right)
\right].\nonumber
\end{align}

At this point we can simply perform the integrations and obtain

\begin{align}
\prod_{i:c_i\in \compsubset}
\prod_{i:c_i \in \compset\setminus \compsubset}
\end{align}

which coincides with the probability of a possible world in <cit.>.

Proposition <ref> can easily be extended to also allow for measures of finite categorical distributions, which then translates to (neural) annotated disjunctions.
Consequently, as is a strict superset of ProbLog [Fierens et al., 2015], also strictly generalises ProbLog.

§ PROOF OF PROPOSITION <REF>

is in essence a subset of the probabilistic logic programming language defined by Gutmann et al., 2011 – the only difference being that the parameters on the right-hand side of a neural distributional fact are not limited to numerical constants any more but can be arbitrary numeric terms.
Under the condition that all NDFs and PCFs are valid, this does, however, not violate any of the assumptions made in <cit.> (proving the measurability of a program).
We can, hence, conclude that a valid program induces a probability measure for $q$.

Note that, similar to ProbLog and , the semantics for are only defined for so-called sound programs [Riguzzi and Swift, 2013], which means that all programs become ground eventually when queried.

§ PROOF OF PROPOSITION <REF>

First, let us consider the indices of the two product expressions in Equation <ref>. We define

\begin{align*}
\negcompsubset \coloneqq \left\{\bar{c}_i\ |\ c_i \in \compset {\setminus} \compsubset\right\}
\end{align*}

such that Equation <ref> can be rewritten as

\begin{align}
P(\world_{\compsubset}) =
\int
\left(\prod_{c_i \in \compsubset {\cup} \negcompsubset}  \indicator(c_i(\variables{x}) \right)
\ \differential P_{\dfacts}
\end{align}

Furthermore, decomposing the measure into a probability density function $\weight(\variables{x})$ and a differential $\differential\variables{x}$ of the integration variables yields

\begin{align}
\int
\left(\prod_{c_i\in \compsubset {\cup} \negcompsubset}  \indicator(c_i(\variables{x})) \right)
\cdot \weight(\variables{x})
\ \differential \variables{x}.
\end{align}

We can now plug this last expression into Equation <ref> resulting in

\begin{align}
P(q) =
\int
\sum_{\substack{C_M \subseteq \mathcal{C}_M:\\ q \in \omega_{\mathcal{C}_M}}}
\left(
\prod_{c_i\in \compsubset {\cup} \negcompsubset}  \indicator(c_i(\variables{x}))
\right)
\cdot \weight(\variables{x})
\ \differential \variables{x}.
\label{eq:proof_inference_as_wmi}
\end{align}

Note that we changed the order of the integration and summation. This operation was shown to be valid in Zuidberg Dos Martires et al., 2019 using de Finetti's theorem. Zuidberg Dos Martires et al., 2019 also showed that the expression in Equation <ref> is indeed a weighted model integral as defined by Belle et al., 2015. Specifically, line P2 in the proof of Theorem 2 in Zuidberg Dos Martires et al., 2019 corresponds to C.3, which is shown to be equal to an instance of WMI.

§ DETAILS ON DERIVATIVE ESTIMATE

To give further details on estimating the derivative we will write the expression $\deriv P_{\Neuralparams}(q)$ in terms of indicator functions

\begin{align}
&\deriv P_{\Neuralparams}(q)
\\
\deriv \int \amc(\variables{x}) \cdot \weight_{\Neuralparams}(\variables{x})\ \partial \variables{x} \nonumber
\\
&= \nonumber
\deriv \int
\sum_{\substack{C_M \subseteq \mathcal{C}_M:\\ q \in \omega_{\mathcal{C}_M}}}
\left(
\prod_{c_i\in \compsubset {\cup} \negcompsubset}  \indicator(c_i(\variables{x}))
\right)
\cdot \weight_{\Neuralparams}(\variables{x})
\ \differential \variables{x}
\end{align}

where the dependency of the probability on the neural parameters $\Neuralparams$ is again made explicit.
Reparametrising the distribution $\weight_{\Neuralparams}(\variables{x})$ yields

\begin{align}
\label{eq:detailedreparamder}
&\deriv P_{\Neuralparams}(q)
= \\
\sum_{\substack{C_M \subseteq \mathcal{C}_M:\\ q \in \omega_{\mathcal{C}_M}}}
\left(
\prod_{c_i\in \compsubset {\cup} \negcompsubset}  \indicator(c_i(\reparam)
\right)
\cdot p(\variables{u})
\ \mathrm{d}\boldsymbol{u}. \nonumber
\end{align}

Explicitly writing out the indicators clearly illustrates the non-differentiability of $\amc(\variables{x})$, which prevents us from applying Leibniz' integral rule [Flanders, 1973] to swap the order of integration and differentiation. To obtain the necessary differentiability of the integrand, the continuous relaxations introduced by Petersen et al., 2021 are utilised. These relaxations allow for comparison formulae of the form

\begin{align}
(g(\boldsymbol{x}) \bowtie 0),
\quad \text{with}
\bowtie\ \in\ \left\{<, \leq, >, \geq, =, \neq\right\}
\end{align}

to be relaxed. We write the continuous relaxation of an indicator function
\indicator(g_i(\variables{x})\bowtie 0)$ as
\relaxation_i(\variables{x})
Four specific cases of relaxations arise, depending on the comparison operator used. Specifically, we define

\begin{align}\label{eq:adaptedsoft}
&\relaxation_i(\variables{x}) =
\\ \nonumber
\sigma(\coolness_i \cdot g_i(\variables{x}))
\text{if $\bowtie\ \in\ \left\{>, \geq \right\}$}, \\
\sigma(-\coolness_i \cdot g_i(\variables{x}))
\text{if $\bowtie\ \in\ \left\{ <, \leq\right\}$}, \\
\prod_{}
\sigma(\coolness_i \cdot g_i(\variables{x})) \cdot \sigma(-\coolness'_i \cdot g_i(\variables{x}))
\text{if $\bowtie\ \in\ \left\{=\right\}$},
\\
1 - \sigma(\coolness_i \cdot g_i(\variables{x})) \cdot \sigma(-\coolness'_i \cdot g_i(\variables{x}))
\text{if $\bowtie\ \in\ \left\{\neq\right\}$},
\end{cases}
\end{align}

where $\coolness_i$ and $\coolness'_i$  are the coolness parameters of the continuous relaxations and $\sigma$ denotes the sigmoid function.
Note that all four cases originate from the root choice of approximating the step function as a sigmoid function. Additionally, this choice is sound as we have that

\begin{equation}
\lim_{\coolness_i \rightarrow +\infty} \sigma(\coolness_i \cdot g_i(\variables{x})) = \indicator(g_i(\variables{x}) \geq 0).
\end{equation}

Continuously relaxing indicator functions using the definition of Equation <ref> renders the integrand differentiable, allowing the application of Leibniz' integral rule and yielding

\begin{align}
&\deriv P_{\Neuralparams}(q)
\approx \\
\nonumber
&\int  \deriv
\sum_{
\substack{\compsubset \subseteq \compset:
\\
q \in \omega_{\compset}}
\left(
\prod_{i: c_i\in \compsubset {\cup} \negcompsubset}  \relaxation_i(\reparam)
\right)
\cdot p(\variables{u})
\ \differential \variables{u}.
\end{align}

The derivative $\deriv P_{\Neuralparams}(q)$ can now be computed using off-the-shelf automatic differentiation software such as PyTorch [Paszke et al., 2019] or TensorFlow [Abadi, 2016], which entails that estimating the gradient  $\nabla_{\Neuralparams} P(q) = (\deriv P(q))_{\lambda \in \variables{\Neuralparams}}$ is computationally as expensive as computing the probability itself, up to a constant factor [Griewank and Walther, 2008].

§ PROOF OF PROPOSITION <REF>

First we express $P(q)$ using Equation <ref>, which we then rewrite without loss of generalisation using only Heaviside distributions[Here we use the term distribution in the sense of a generalised function [Schwartz, 1957] and not in the sense of a probability distribution.].

\begin{align}\label{eq:prob_query_heavi}
&P(q) \\
&= \nonumber
\int
\sum_{\substack{C_M \subseteq \mathcal{C}_M:\\ q \in \omega_{\mathcal{C}_M}}}
\left(
\prod_{c_i\in \compsubset {\cup} \negcompsubset}  \indicator(c_i(\variables{x}))
\right)
\cdot \weight(\variables{x})
\ \differential \variables{x} \\
&= \nonumber
\int
\sum_{\substack{C_M \subseteq \mathcal{C}_M:\\ q \in \omega_{\mathcal{C}_M}}}
\left(
% \prod_{c_i \in \compsubset {\cup} \negcompsubset}
\prod_{g_i \in \Sigma_{\compsubset {\cup} \negcompsubset}}
\right)
\cdot p(\variables{u})
\ \mathrm{d}\boldsymbol{u}
\end{align}

In the Equation above, $H(x)$ denotes the Heaviside distribution and $\Sigma_{\compsubset {\cup} \negcompsubset}$ is the set of all sigmoid functions involved in the continuous relaxations of the set $\compsubset {\cup} \negcompsubset$.

This rewrite is possible as the indicator function of any PCF $c(\variables{x})$ is either a step function or decomposes into a product of step functions. Indeed, if $c(\variables{x})$ is of the form $g(\boldsymbol{x}) \geq 0$, then $\indicator(c(\variables{x})) = H(g(\boldsymbol{x}))$. If it is of the form $g(\boldsymbol{x}) = 0$, then $\indicator(c(\variables{x})) = H(g(\boldsymbol{x})) \cdot H(- g(\boldsymbol{x}))$. The other cases with different comparison operators follow from these two.

Differentiating in a distributional sense and applying Leibniz' integral rule [Flanders, 1973] then yields

\begin{align}\label{eq:heavisideder}
\sum_{\substack{C_M \subseteq \mathcal{C}_M:\\ q \in \omega_{\mathcal{C}_M}}}
\sum_{g_j \in \Sigma_{\compsubset {\cup} \negcompsubset}}
\int
&\deriv H(g_j(\reparam)) \cdot \\ \nonumber
&\prod_{i \neq j} H(g_i(\reparam))
\cdot p(\variables{u})
\ \mathrm{d}\boldsymbol{u}.
\end{align}

We can reduce the discussion by considering each term in this equation separately, because of the linearity of the integral. In other words, to prove our statement, it suffices to show that

\begin{align}\label{eq:heavisideterm}
\int
&\deriv H(g_j(\reparam)) \cdot \\ \nonumber
&\prod_{i \neq j} H(g_i(\reparam))
\cdot p(\variables{u})
\ \mathrm{d}\boldsymbol{u},
\end{align}

is equal to

\begin{align}
\lim_{\coolness_1, \dots, \coolness_n\rightarrow +\infty}
\int
&\deriv \sigma(\coolness_j \cdot g_j(\reparam)) \cdot \\ \nonumber
&\prod_{i \neq j} \sigma(\coolness_i \cdot g_i(\reparam))
\cdot p(\variables{u})
\ \mathrm{d}\boldsymbol{u}.
\label{eq:sigmaterm}
\end{align}

For brevity's sake, we will write the products

\begin{align}
\prod_{i \neq j} H(g_i(\reparam))
\end{align}

\begin{align}
\prod_{i \neq j} \sigma(g_i(\reparam)),
\end{align}

as $\pi_j(\variables{u})$ and $\pi_j^\sigma(\variables{u})$, respectively.
Next, using distributional notation, Equation <ref> can be further simplified as

\begin{align}
&\left\langle \deriv (H \circ g_j \circ r),\ \pi_j \cdot p \right\rangle
= \\ \nonumber
&\left\langle \delta \circ g_j \circ r,\ \deriv \left(g \circ r\right) \cdot \pi_j \cdot p \right\rangle.
\end{align}

Note that this expression utilises the assumption that $\deriv (g_j \circ r) \in L^1_{loc}(\mathbb{R}^k)$, i.e., $\deriv (g_j \circ r)$ is locally integrable over $\mathbb{R}^k$.
This asssumption is not very demanding, since distributions (generalised functions) are only well-defined when acting on functions that are at least locally integrable. Equation <ref> can similarly be rewritten and simplified to obtain the equality

\begin{align}
&\phantom{={}} \lim_{\coolness_1, \dots, \coolness_n\rightarrow +\infty} \left\langle \deriv (\sigma \circ g_j \circ r),\ \pi_j^\sigma \cdot p \right\rangle\\
&= \nonumber
\lim_{\coolness_1, \dots, \coolness_{j - 1}, \coolness_{j+1}, \dots, \coolness_n\rightarrow +\infty} \\
\nonumber
\left\langle \delta \circ g_j \circ r,\ \deriv (g\circ r) \cdot \pi_j^\sigma \cdot p \right\rangle.
\end{align}

More explicitly,

\begin{align}
\lim_{\coolness_1, \dots, \coolness_n\rightarrow +\infty} \\ \nonumber
\int
\deriv \sigma(\coolness_j \cdot g_j(\reparam)) \cdot
\pi_j^\sigma(\variables{u}) \cdot p(\variables{u})
\ \differential \variables{u} \\
&= \nonumber
\lim_{\coolness_1, \dots, \coolness_n\rightarrow +\infty} \\ \nonumber
\int
\frac{l\cdot e^{-g(\reparam)\cdot \coolness_j}}{(1 + e^{-g(\reparam)\cdot \coolness_j})^2}\cdot \\
&\phantom{={}\int} \nonumber
\deriv g_j(\reparam) \cdot
\pi_j^\sigma(\variables{u}) \cdot p(\variables{u})
\ \differential \variables{u} \\
&= \nonumber
\lim_{\coolness_1, \dots, \coolness_{j - 1}, \coolness_{j+1}, \dots, \coolness_n\rightarrow +\infty} \\ \nonumber
\int
\delta(g_j(\reparam)) \cdot \\
&\phantom{={}\int} \nonumber
\deriv g_j(\reparam)\cdot
\pi_j^\sigma(\variables{u}) \cdot p(\variables{u})
\ \differential \variables{u}.
\end{align}

The last transition uses the fact that

\begin{equation}
\lim_{\coolness_j \rightarrow +\infty} \frac{\coolness_j \cdot e^{-g(\reparam)\cdot \coolness_j}}{(1 + e^{-g(\reparam)\cdot \coolness_j})^2}
= \delta(g(\reparam)),
\end{equation}

in the distributional sense. In addition, we also have (again in distributional sense) that

\begin{align}
\lim_{\coolness_i \rightarrow +\infty} \sigma(\coolness_i \cdot g_i(\reparam))
= H(g_i(\reparam)).
\end{align}

This final equation allows us to replace $\pi_j^\sigma(\variables{u})$ in the final line of Equation <ref> with $\pi_j(\variables{u})$ by repeating the above steps for each index $i$ separately. Hence, we can conclude that our relaxation of $\deriv P(q)$ is indeed unbiased in the infinite coolness limit.

§ EXPERIMENTAL DETAILS

This section will give detailed programs, neural network architectures and elaborated figures for each of the experiments present in the main body of the paper.

All experiments were run on an RTX 3080 Ti coupled with a Intel Xeon Gold 6230R CPU @ 2.10GHz and 256 GB of RAM, except the LTN results of Section <ref>.
Note that the optimisation of any hyperparameters, such as learning rate or number of training epochs, was done via a grid search on a separate validation set.

§.§ NeSy attention

Setup details and program.
The full program for the detection of handwritten years is given in Listing <ref>.
The query year is optimised for a different number of samples depending on the experiment. For , we have 28 000 training samples while there are 4000 validation and 8000 test samples. The set of years in the training, validation and test set are disjoint. For , the size of validation and test set is the same as in the case of , but with a training set of 40 000 samples. Here, the set of years in validation and test set are disjoint, but both are a subset of the set of years of the training set.

box(Params, B)
digit(Im, Loc, D)
categorical(classifier([Im, Loc]),
[0, ..., 9]).
year(Im, Year1, Year2, Year3, Year4) :-
region(Im, [Y1, Y2, Y3, Y4]),
ordered_output([Y1, Y2, Y3, Y4]),
box(Y1, B1), box(Y4, B4),
x_diff(0.0, B1, B1diff),
B1diff < 0,
x_diff(1.0, B4, B4diff),
0 < B4diff,
digit(Im, Y1, D1), digit(Im, Y2, D2),
digit(Im, Y3, D3), digit(Im, Y4, D4),
Year1 =:= D1, Year2 =:= D2),
Year3 =:= D3, Year4 =:= D4.

ordered_output([[Mu, Sigma]]).
ordered_output([[Mu, Sigma], H2 | T]) :-
box([Mu, Sigma], B1), box(H2, B2),
x_diff(B1, B2, Bdiff),
Bdiff < 0,
ordered_output([H2 | T]).

There is one continuous NDF, box, which represents a bounding box as a generalised normal distribution with mean and scale being the center and width of the box, respectively. digit is a discrete NDF that denotes the categorical distribution of the digit classifications made by the network classifier.
region is the detection network that predicts the 4 bounding boxes, i.e., the parameters of four instances of box.
Given these parameters, the predicate ordered_output will enforce the spatial constraints that region predicts its boxes in order from left to right on the image. It does so by taking the difference of the $x$ coordinate of each subsequent bounding box, which is a 2-dimensional random variable, and employing a $<$ PCF.
Finally, the supervision on the digits of the year is given to the correct bounding box.

Parameters and neural architectures.
A schematic overview of the neural architecture used for all different methods can be seen in Figure <ref>.
The neural baseline simply outputs the four predictions of the classification network and optimises them by minimising the categorical cross-entropy on each digit of the year.
In the case of the neural-symbolic methods, the output of both the regression and classification components are used in the logic. optimises a binary cross-entropy on the probability of year, while LTN optimises the MAX-sat objective function.
As optimiser, we utilised Adamax [Kingma and Ba, 2015] with its default learning rate of $10^{-3}$. and LTNs were run for 10 epochs, while the neural baseline was given 20 epochs, all with a batch size of 10. This number of epochs proved sufficient for all methods to converge.
Interestingly, no special annealing scheme was necessary for this experiment as constant value of $50$ for the coolness parameters lead to satisfactory results.
All these hyperparameters were determined through a grid search on the validation set.

Overall neural architecture for the dates experiment. Following the orange arrows first, the parameters of 4 generalised normal distributions are predicted for each image. Then, following the green arrows, the images is attenuated separately by each of the 4 distributions and then classified as a digit between 0 and 9 to give the total overall year prediction as an ordered tuple of 4 digits.

Additional results and interpretations.

Roughly speaking, every 100 iterations took about 25 seconds for while the neural baseline took around 15 seconds. Given the results of Table <ref> and 's satisfactory solution to the problem, the additional computational cost of adding probabilistic logic is worthwile in this case.

§.§ Neural hybrid Bayesian network

Setup details and program.
Our encoding of the neural hybrid Bayesian network of experiment <ref> is given in Listing <ref>. The goal is to optimise the neural networks responsible for the classification of humid and cloudy conditions, as well as the network that predicts the temperature value. Additionally, we explicitly model the noise present on the true temperature labels as a learnable parameter. To achieve this, a set of 1200 triples (Im1, Im2, X) are used as training set, where Im1 is a CIFAR-10 image belonging to one of the first three classes, while Im2 belongs to the last two classes. In other words, we use CIFAR-10 images as proxies for real imagery data. X is a set of 25 numerical meteorological features sampled from a publicly available Kaggle dataset [Cho et al., 2020]. The label of each triple is the probability that the weather, as described by the correct labels of humid, cloudy and temperature, is good. Computing this probability label is non-trivial in itself. We utilised a large set of 1000 samples to approximate the correct underlying distributions and to obtain an approximate probability label.

humid(Im, H)
cloudy(Im, C)
[0, 1, 2]).

temperature(X, T)
normal(temperature_detector(X), t(_)).
snowy_pleasant   beta(11, 7).
rainy_pleasant   beta(1, 9)
cold_sunny_pleasant   beta(1, 1).
warm_sunny_pleasant   beta(9, 2).

rainy(I1, I2) :-
cloudy(I1, C), C =0̄,
humid(I2, H), H =:= 1.

good_weather(I1, I2, X) :-
rainy(I1, I2),
temp(X, T), T < 0,
snowy_pleasant > 0.5.
good_weather(I1, I2, X) :-
rainy(I1, I2),
temp(X, T), T >= 0,
rainy_pleasant > 0.5.
good_weather(I1, I2, X) :-
rainy(I1, I2),
temp(X, T), T > 15,
warm_sunny_pleasant > 0.5.
good_weather(I1, I2, X) :-
rainy(I1, I2),
temp(X, T), T <= 15,
cold_sunny_pleasant > 0.5.

P :: depressed(I1) :-
cloudy(I1, C), C =:= N,
P is N * 0.2.

enjoy_weather(I1, I2, X) :-
good_weather(I1, I2, X).

The NDFs humid and cloudy classify a given image as describing humid and cloudy conditions, respectively. temp takes a set of 25 numerical features and predicts a mean temperature from those. Note that t(_) is ProbLog notation for a single optimisable parameter. Depending on the value of the temperature, 4 different cases of weather and their degree of pleasantness are described by beta distributions. We define good_weather as being true if the degree of pleasantness of any case is larger than 0.5. Finally, a person can be depressed with probability 0.2 or 0.4 depending on the degree of cloudy. Both then determine whether a person can enjoy the weather, if they are not depressed and good_weather is the case.

Parameters and neural architectures.
We utilise simple classifiers (Figure <ref>) in the NDFs cloudy and humid, while the network in the neural predicate temperature has three dense layers of size 35, 35 with ReLU activations and 1 with linear activation. Both classifiers share a common set of convolutional layers, requiring the learning of features that generalise to both classification problems. Additionally, the noise on the temperature prediction is modelled explicitly as a learnable TensorFlow variable with an initial value of 10. This choice is not arbitrary, as the initial neural parameter estimate will hover around the middle of the possible temperature values and a choice of 10 as initial standard deviation allows covering the entire range of temperature values with a non-insignificant probability mass. In this way, gradient information across the entire temperature domain can be accumulated during learning. Finally, was trained for 10 epochs using Adamax with learning rate $10^{-3}$ and batch size of 10.

Ideally, simple 0-1 labels of enjoy_weather would be more intuitive, as we often do not observe the probability of an event but single cases where it is either true or false. However, our experiments have showed that our small dataset is insufficient to find an optimal solution using such labels in conjunction with the very distant supervision. To show that is still able to find solutions in cases where the supervision is slightly less distant using only 0-1 labels, we added a different neural hybrid Bayesian network experiment in Section <ref> based on the well-known burglary-alarm example of probabilistic logic.

Additional results and interpretations.
We want to stress that learning to predict the right mean temperature from the distant supervision is not straightforward. The only learning signal for the temperature has to pass through PCFs with a very wide range, meaning they do not specify the exact temperature value. Additionally, these PDFs still do not directly influence the supervision of enjoy_weather, only good_weather.
The Gaussian noise that renders the temperature into a continuous random variable only further convolutes the task. We conclude that can extract meaningful learning signals from reasonably distant supervision.

§.§ Neural-symbolic variational autoencoder

Setup details and programs.
Each data sample consists of 2 regular MNIST digits and the result of their subtraction. The first digit takes the place of the minuend while the second one is interpreted as the subtrahend. The training, validation and test sets had 30 000, 1 000 and 1 000 samples of this form, respectively. Encoding a VAE without additional logic in is straightforward (Listing <ref>), while adding logic involves more engineering freedom (Listing <ref>). We opted for the simplest use of a conditional variational auto-encoder by only using the classified digit as additional input to the decoder. Note that during optimisation, both the VAE and digit classifier are trained jointly.

prior(P)   normal(0, 1).
latent(Im, L)

good_image(Image) :-
prior(P), latent(Im, L),
P =:= L,
decoder_net(L, G),
soft_unification(G, Image).

Prototypical implementation of a Gaussian VAE in .
A normal prior prior is used to regularise a Gaussian latent space modelled by the second NDF by expressing that they should be equal.
The decoder component of the VAE is given by decoder_net and returns a generated image G by sampling the latent space. This generation is self-supervised by soft unifying it with the given image.

prior(ID, P)   normal(0, 1).
digit(Emb, D)
[0, ..., 9]).
latent(Im, L)

good_subtraction(Im1, Im2, Diff) :-
prior(1, P1), prior(2, P2),
latent(Im1, L1), latent(Im2, L2),
L1 =:= P1, L2 =:= P2,
embedding(Im1, E1),
embedding(Im2, E2),
digit(E1, D1), digit(E2, D2),
Diff =:= D1 - D2,
concat(L1, D1, ConditionalL1),
concat(L2, D2, ConditionalL2),
decoder_net(ConditionalL1, G1),
decoder_net(ConditionalL2, G2),
soft_unification(G1, Image1),
soft_unification(G1, Image1).

Combining subtraction logic with a VAE in .
Each image is encoded into a Gaussian latent space and embedded into a lower-dimensional real space. The latent space is regularised by the standard normal prior while the embedding forms the input to a digit classifier to find which digit is on the image. The two classified digits, which follow a categorical distribution, should subtract to the given value of Diff. Finally, the Gaussian latent space and the categorical digits are concatenated into the conditional latent space of the CVAE. The decoder network again samples from this space to construct a generation for both images, which should softly unify with the original images.

Parameters and neural architectures.  The NeSy VAE has two main neural components (Figure <ref>), one for the VAE itself and another that handles the digit classification used in the subtraction logic. A small set of 256 samples with direct supervision on the digit labels is used to pre-train the classification portion of the overall network to avoid degenerate solutions. All training utilised Adam as optimiser with a learning rate of $\cdot 10^{-3}$ and took 20 epochs using a batch size of 10. The pre-training was given 1 epoch with a batch size of 4.

VAE encoder-decoder architecture. The decoder is equivalent to the transpose of the encoder. All layers use ReLU activation functions, except the final convolutional one, which applies a hyperbolic tangent.

Regular Gaussian VAE optimisation has two components: a Kullback-Leibler (KL) divergence term and a reconstruction loss term. Since requires probabilistic values, i.e., between 0 and 1, a probabilistic translation of these terms is necessary for optimisation in . The KL divergence term compares the latent distribution of the VAE to a standard normal prior and can as such be replaced by a =:= comparison in the logic. The reconstruction loss is chosen to be the exponentiation of a negated average $L_1$ loss function, as it yields a value between 0 and 1 that can be interpreted as the probability that two images match. Specifically, the loss between two such images $\boldsymbol{I_1}, \boldsymbol{I_2} \in \mathbb{R}^{768}$ is given by

\begin{equation}
\exp(-\frac{1}{768}\sum_{i=1}^{768} \left|I_{1i} - I_{2i}\right|).
\end{equation}

The latter can be interpreted as a form of soft unification [Rocktäschel and Riedel, 2017], which is why we denote it by the predicate soft_unification.

Additional results and interpretations.
Emphasis has to be put on the flexibility of generation in , as the generation of digits can be carried out in a range of different contexts without further optimisation. One only needs to write a query describing that logical context. The query that yields an image of both a left and right digit that subtract to a given value is given in Listing <ref>. The conditional query that generates an image of a right digit given an image of the left digit and their difference value is given in Listing <ref>.

generate_subtraction(G1, G2, Diff) :-
member(D1, [0, ..., 9]),
member(D2, [0, ..., 9]),
prior(1, P1), prior(2, P2),
Diff =:= D1 - D2,
concat(P1, D1, ConditionalL1),
concat(P2, D2, ConditionalL2),
decoder_net(ConditionalL1, G1),
decoder_net(ConditionalL2, G2).

The logic finds all possible combinations for D1 and D2 that meet the subtraction evidence Diff and concatenates these to a standard normal prior component into the conditional latent space. The decoder then generates images from a sample of this space.

generate_left(RightIm, Diff, LeftG) :-
member(D1, [0, ..., 9]),
embedding(RightIm, RightE),
digit(RightE, RightD)
Diff =:= LeftD - RightD,
latent(RightIm, RightL),
concat(RightL, LeftD, LeftCondL),
decoder_net(LeftCondL1, LeftG).

Given an image of the right digit and a difference value, we generate an image of a the left digit. The right's image is classified such that the logic can find the value of LeftD that meets the given difference. By attaching that value to the Gaussian latent space of the right digit, the VAE can generate an image of the correct left digit in the style of the right one.

§ ADDITIONAL EXPERIMENT

An additional experiment was performed to show the promise of discrete-continuous neural probabilistic logic programming.
It is similar to the neural hybrid Bayesian network of Section <ref>, but with more practical 1-0 query supervision.

§.§ Neural-continuous burglary alarm

Setup details and program.
The neural-continuous burglary alarm (Listing <ref>) extends the classic example from Bayesian network literature (Listing <ref>).

0.1 :: earthquake.
0.3 :: burglary.
0.9 :: hears.

0.7 :: alarm :- earthquake.
0.9 :: alarm :- burglary.

calls :- alarm, hears.

Classical burglary-alarm ProbLog program. Three probabilistic facts earthquake, burglary and hears are given with their probabilities. A neighbour calls when hearing an alarm, while an alarm can go off because of an earthquake or a burglary.

Each data sample is a triple $(E, B, L)$, where $E$ can be an MNIST digit 0, 1 or 2 while $B$ can be an MNIST 8 or 9. Values for $E$ of 0, 1 and 2 correspond to no earthquake, a mild earthquake or a heavy earthquake respectively. If $B$ is an MNIST 8, then there is no burglary. If it is 9, then there is a burglary. $L$ can have either the value 0 or 1, indicating whether the neighbour called or not. Our dataset contains 12 000 such triples for training, while having 1 000 for validation and 2000 for testing purposes. Obtaining the weak supervision $L$ is done by sampling according to the true probability of calling given the input To compute this true probability, a single sample is taken from the neighbour's true distribution. This true distribution has respective means of 6 and 3 for the horizontal and vertical Gaussian while both directions have a standard deviation of 3. Additionally, there are two possible ways to express that the distance of the neighbour should be smaller than 10 distance steps before hearing the alarm. One can use either the squared distance or the true distance in the rule hears. A separation is often maintained in the weighted model integration literature [Zuidberg Dos Martires et al., 2019] between comparison formulae that are polynomial and those that are generally non-polynomial. To illustrate that can deal with both classes of formulae, we will perform experiments for both the squared distance (polynomial, Listing <ref>) and the true distance (non-polynomial, Listing <ref>). Both these functions are implemented in Python and allows them to be easily imported as built-in predicates.

earthquake(Im, E)
[0, 1, 2]).
burglary(Im, B)
[8, 9]).

normal([t(@$\mu_x$@), t(@$\mu_y$@)],
[t(@$\sigma_x$@), t(@$\sigma_y$@)]).

hears :-
squared_distance(0, N, D),
D < 100.

P :: alarm(Im1, _) :-
earthquake(Im1, E), E =:= N,
P is N * 0.35.
0.9 :: alarm(_, Im2) :-
burglary(Im2, B), B =:= 9.

calls(Im1, Im2) :-
alarm(Im1, Im2), hears.

Our extension of the burglary alarm example has two categorical NDFs that model the chance of an earthquake and a burglary given an image. Additionally, whether the neighbour can hear the alarm if it goes off depends on their spatial distribution, which is modelled as a two-dimensional Gaussian distribution. This distribution is randomly initialised and its parameters need to be optimised.

hears :-
distance(0, N, D),
D < 10.

Using the true distance in the hears predicate as a case of a non-polynomial comparison formula.

Parameters and neural architectures. The complete neural architecture of both the earthquake and burglary classifiers is given in Figure <ref>. In addition to the neural parameters in these networks, four independent parameters are present in the program. These are used as the means and standard deviations for the neighbour's spatial distribution and are randomly initialised. Specifically, the means are sampled uniformly from the interval $\left[0, 10\right]$ while the standard deviations were sampled from $\left[2, 10\right]$. All optimisation was performed using regular stochastic gradient descent with a learning rate of $8\cdot 10^{-2}$ for two epochs using a batch size of 10.

Overview of the architecture of the earthquake and burglary networks. Both share two convolutional layers, but each specific network applies its own final convolutional layer followed by a global average-pooling operation with softmax activation. All other activation functions are ReLUs.

Because of the difference in nature between the parameters in the neural networks and the four independent parameters in the Gaussian distribution, the latter required a boosted learning rate to provide consistent convergence. Specifically, the gradients for these four parameters were multiplied by a value of 20, which was found by a hyperparameter optimisation on the validation set.

Results and interpretation.
Initial learning progress of the neural networks seems volatile (Figure <ref>), which is likely due to the unoptimised state of the neighbour's spatial distribution. Two epochs of training proves to be sufficient to optimise both the neural detectors and the neighbour's distribution.

The 4 parameters of the neighbour's distribution do not converge to the true values, which is to be expected as their supervision is underspecified. However, they do converge to values that result in PCF probabilities that are close to those of the true underlying distribution. All in all, three conclusions can be drawn. First, this experiment indicates that is capable of jointly optimising neural parameters and independent, distributional parameters. Second, seems to be able to fully exploit both polynomial and more general non-polynomial comparison formulae. It shows the strength of our approximate approach, as exact methods often fail to efficiently deal with non-polynomial formulae [Zuidberg Dos Martires et al., 2019]. Third, can deduce meaningful probabilistic information from weak labels. Indeed, in order to optimise the neural detectors and the neighbour's distribution, has to aggregate meaningful update signals from the 0-1 labels across the given training data set to approximate the underlying probability of calls.

Evolution of the training loss and validation accuracy of the neural earthquake and burglary detectors. For both squared (left) and true distance (right), the discrete supervision seems to be sufficient to facilitate meaningful learning.

§ LIMITATIONS

The main limitation of is one that it inherits from probabilistic logic in general, computational tractability. Efficiently representing a probabilistic logic program is done via knowledge compilation, which is $\#P$-hard. Once the probabilistic program is knowledge compiled, evaluating the compiled structure is linear in the size of this structure.
Inference remains linear in the size of the compiled structure after the addition of continuous random variables as all samples can be run in parallel with the current inference algorithm.

Although our sampling strategy is efficient in the sense that it is linear in the number of samples, uses the advanced inference techniques of Tensorflow Probability to effectively sample higher dimensional distributions, and it can be executed in parallel for each sample, it remains ignorant of the comparison formulae that are approximated. More intricate inference strategies exist within the field of weighted model integration [Morettin et al., 2021], yet they currently lack the differentiability property to be integrated in 's gradient-based optimisation. Conversely, our examples illustrate that our rather naive strategy is sufficient to solve basic tasks.

It is still an open question how to perform successful joint inference and gradient-based learning under general comparisons.

Orthogonal to the estimation of the integral during inference, exact knowledge compilation also prevents the scaling of to larger problem instances. Approximate knowledge compilation is the field of research that deals with tackling this issue. While it contains interesting recent work [Fierens et al., 2015, Huang et al., 2021, Manhaeve et al., 2021], it was highlighted by [Manhaeve et al., 2021] that the introduction of the neural paradigm does lead to further complications. As such, we opted for exact knowledge compilation, but it has to be noted that we will be able to benefit from any future advances in the field of approximate inference. Alternatively, different semantics [Winters et al., 2022] can simplify inference, but they lead to a degradation of expressivity of the language.

A potential future avenue for scaling up inference would be the use of further continuous relaxation schemes. More specifically, replacing discrete random variables with relaxed categorical variables [Maddison et al., 2017, Jang et al., 2017] might allow us, for instance, to forego the knowledge compilation step while still being able to pass around training signals.