# Curved Space-Filling Tiles Using Voronoi Decomposition with Line, and Curve
Segments Closed Under Wallpaper Symmetries

Haard Panchal
Visualization Department,
Texas A&M University, College Station, TX, 77831
<EMAIL_ADDRESS>
Ergun Akleman
Visual Computing & Computational Media,
Texas A&M University, College Station, TX, 77831
<EMAIL_ADDRESS>
Vinayak Krishnamurthy
Department of Mechanical Engineering,
Texas A&M University, College Station, TX, 77831
<EMAIL_ADDRESS>Tolga Talha Yildiz
Computer Science and Engineering Department,
Texas A&M University, College Station, TX, 77831
<EMAIL_ADDRESS>Varda Grover
Visual Computing & Computational Media,
Texas A&M University, College Station, TX, 77831
<EMAIL_ADDRESS>
Joint with Computer Science and Engineering Department.Joint with Computer
Science and Engineering Department.

###### Abstract

In this paper, we present a new approach to obtain symmetric tiles with curved
edges. Our approach is based on using higher-order Voronoi sites that are
closed under wallpaper symmetries. The resulting Voronoi tessellations provide
us with symmetric tiles with curved edges. We have developed a web application
that provides real-time tile design. Our application can be found at
https://voronoi.viz.tamu.edu/. One of our key findings in this paper is that
not all symmetry operations are useful for creating curved tiles. In
particular, all symmetries that use mirror operation produce straight lines
that are useless for creating new tiles. This result is interesting because it
suggests that we need to avoid mirror transformations to produce unusual
space-filling tiles in 2D and 3D using Voronoi tessellations.

(a) An example of curved tiles with p3 symmetry.

(b) Same tiles with Escher-like hand-drawn faces.

(c) Same tiles interlocked together.

(d) Same tiles filling 2D space.

Figure 1: Using our system artists can create Escher-like curved congruent
tiles. In this example, we demonstrated how a congruent tile was created in
our system using p3 symmetry can later be painted to obtain the look of
Escher-style animal drawings that can tesselate the 2D space.

##  1 Introduction and Motivation

With the development of the internet and computer graphics, many interactive
web-based systems to design symmetric patterns have been developed. Most web-
based symmetric pattern design systems are really painting systems, that do
not produce tiles. One alternative approach is to develop them as drawing
systems. There currently exist two types of drawing systems: (1) The stand-
alone systems that allow for manipulating polygons based on Escher’s method of
creating space-filling tiles such as Tesselmaniac Tessellations (2015) and
Tesselation Exploration; and (2) the web-based systems that allow drawing
graphs based on 17 wallpaper symmetries, such as Kali Amenta (1995), Symmetric
Graph Drawing system Akleman et al. (2000), and EscherSketch Levskaya (2017).
The advantage of the first approach is that it guarantees watertight
tessellations without gaps and overlaps. However, it only allows for a change
of the geometric shape of the polygonal tiles Tessellations (2015).

(a) Patterns with p1 symmetry.

(b) Patterns with p2 symmetry.

(c) Patterns with pg symmetry.

(d) Patterns with p4 symmetry.

Figure 2: Examples of curved tiles created by using single lines that are
closed under p1, p2, pg, and p4 symmetries as Voronoi sites. The lines that
are used as Voronoi sites are shown as light gray. These are the only ones of
square-based symmetries (total 12)that can be used effectively to obtain
curved boundaries with Voronoi decomposition. Other eight wallpaper
symmetries, i.e. pm, cm, pmg, pgg, cmm, p4g, p4m, and pmm, include mirror
operation and therefore, produce tiles with some straight boundary edges. In
other words, these eight symmetries do not produce desired curved edges.

(a) Patterns with p3 symmetry.

(b) Patterns with p6 symmetry.

Figure 3: Examples of curved tiles created by using single lines that are
closed under p3, and p6 symmetries as Voronoi sites. The lines that are used
as Voronoi sites are shown as light gray. These are the only ones of hexagon-
based symmetries (total 5) that can be used effectively to obtain curved
boundaries with Voronoi decomposition. Other wallpaper symmetries, p31m, p3m1,
and p6m, include mirror operation and therefore, do not produce interesting
tiling.

(a) Patterns with p2 symmetry.

(b) Pattern with p4 symmetry.

(c) Pattern with p3 symmetry.

Figure 4: More complicated examples using Voronoi sites that are a union of a
set of lines. Using multiple lines more complex tiles can be obtained,

(a) Example pattern with pm symmetry.

(b) Example pattern with cm symmetry.

(c) Example pattern with pmg symmetry.

(d) Example pattern with pgg symmetry.

Figure 5: Examples that show the four symmetries that are not very useful
since they cannot make all boundaries curved. Note that these four symmetries
always produce in straight infinite lines regardless of how we choose Voronoi
sites.

(a) Example pattern with cmm symmetry.

(b) Example pattern with p4g symmetry.

(c) Example pattern with p31m symmetry.

Figure 6: Examples that show the three symmetries that provide a decomposition
of polygons. Two tiles in cmm create squares. Four tiles in p4g also create
squares. Three tiles in p31 symmetry form equilateral triangles.

(a) Example pattern with p4m symmetry.

(b) Example pattern with pmm symmetry.

(c) Example pattern with p13m symmetry.

(d) Example pattern with p6m symmetry.

Figure 7: Examples that show four symmetries that always create the same
polygon regardless of how many lines we use as Voronoi sites. The p4m symmetry
only creates Right Isosceles Triangles. The pmm symmetry does not create
anything beyond square packing. The p13m symmetry creates only Equilateral
Triangles. The p6m symmetry creates right triangles that divide regular
hexagons into 12 congruent shapes.

The advantage of the second approach is that it can allow changes to the
number of sides of polygonal tiles and allow the creation of multiple polygons
Amenta (1995); Akleman et al. (2000). Because of this flexibility, this
approach appears to be the appropriate choice for designing symmetric tiles.
However, tiles are not just graphs, which are a set of arbitrary vertices,
lines, and polygons. The tiles are cellular space-filling structures that come
together without gaps and overlaps. If we do not guarantee the avoidance of
gaps and overlaps, the user must be cautious not to include any gap or overlap
when designing symmetric tiles.

In this work, we present a new approach to developing flexible symmetric
drawing systems that can guarantee no gaps and overlaps. Our approach is based
on the Voronoi tessellations of 2D space by using the Voronoi sites closed on
symmetry operations. This approach ensures that we obtain tiles that come
together without gaps and overlaps or special attention from the user. This is
really a generalization of Delaunay’s original intention for the use of the
Delaunay diagrams. Delaunay was the first to use symmetry operations on points
and the Voronoi tessellations to produce space-filling polyhedra, which he
called Stereohedra Delaunay and Sandakova (1961); Schmitt (2016).

In this paper, we have extended his approach using higher-order Voronoi sites,
i.e., lines, polylines, and curves, to construct space-filling tiles in 2D
with curved boundaries (See Figures 2. and 3 for examples of using lines that
are closed under some symmetry operations as Voronoi sites). Generally, 2D
space is divided into regions so that each point $P$ in the region is closest
to the Voronoi site $S$, which is also a point. Here, the distance is the
length of the line that joins the point $P$ and the Voronoi site $S$. In our
case, each Voronoi site is essentially a line segment $AB$, and space is
tessellated using the closest distance of a point $P$ in the region to the
line segment $AB$ instead of the conventional point. We assume the distance of
a point $P$ to a line segment $AB$ is the length of the shortest line that can
be drawn that connects the line segment $AB$ and the point $P$. We obtain
Voronoi decompositions with polyline and curve Voronoi sites by using unions
of line-based distance functions.

##  2 Related Work

A space-filling shape is a cellular structure whose replicas together can fill
all of the space watertight, i.e., without having any voids between them Loeb
(1991). Equivalently, a space-filling shape is a cellular structure that can
generate a tessellation of space Grünbaum and Shephard (1980). The 2D space-
filling shapes have been widely used in Architecture since antiquity. They are
frequently used by almost every civilization as wallpapers, wall decorations,
ceilings, floor tiles, street pavements, and even facades of buildings Akleman
et al. (2000). The symmetric patterns in Alhambra, Granada are probably the
most well-known architectural usage of 2D space-filling shapes Du Sautoy
(2008).

In this work, we developed an approach to provide a web-based space-filling
symmetric tile design system. Two interesting cases of 2D tilings relevant to
our approach are those presented by Kaplan Kaplan (2000) showing a wide
variety of artistic patterns using specific Voronoi site configurations and
Rao Rao (2017) that show a systematic construction of 2D pentagonal tilings.
In a sense, our work expands on these two works to move beyond polygonal
tilings in 2D space to a rich design of curved symmetric tilings.

This work can also be a specialized version of recent works on the
construction of 2.5D and 3D space-filling structures using higher-order
Voronoi sites, namely Delaunay Lofts Subramanian et al. (2019), Generalized
Abeille Tiles Akleman et al. (2020), Bi-axial woven tiles Krishnamurthy et al.
(2021); Yildiz et al. (2023), and Voronoodles Ebert et al. (2003); Mullins et
al. (2022). modular These more powerful approaches that are designed for
higher dimensional tiles Krishnamurthy et al. (2023) are not reachable for
most people. The systems require higher dimensional Voronoi decomposition and
it is hard to turn them into an interactive web-based system. This particular
work provides an interactive 2D tile design with a real-time interface. The
system is available at https://voronoi.viz.tamu.edu/.

The theoretical foundation for our curved tiles comes from the concept of
symmetric Delone sets, attributed to Boris Delone (Delaunay). The original
(non-symmetric) Delone sets essentially consist of a well-spaced set of
points. Let $S$ denote a set of points in 2D Euclidean space, $\Re^{2}$. The
$S$ is called a Delone set if it is both uniformly discrete and relatively
dense Delone et al. (1970). Formally, let $r_{1}>r_{0}>0$ be two positive
numbers. $S$ is uniformly discrete if each ball of radius $r_{0}$ contains at
most one point in $S$. $S$ is relatively dense if every ball of radius $r_{1}$
contains at least one point of $S$ Delone et al. (1976). If we use the points
in $S\in\Re^{2}$ as Voronoi sites, one would obtain 2-Honeycombs that contain
similar-sized convex polyhedra as Voronoi cells. Owing to this property,
Delone sets, and related Meyer sets have been used to define tessellations
Lagarias (1999a, b); Lagarias and Pleasants (2003); Lagarias (1996).

The symmetric Delone sets in 2D are invariant with respect to one of the 17
wallpaper groups Dolbilin (2015). Therefore, an ideal space filling 2D tiles
can be described by the Delone sets in 2D Euclidean space along with a
wallpaper group of Euclidean isometries acting on this point Delaunay and
Sandakova (1961). A Delone set $S$ is symmetric in 2D if, for every two points
$\mathbf{p},\mathbf{q}\in S$, there exists a rigid motion of 2D space that
takes $S$ to $S$ and $\mathbf{p}$ to $\mathbf{q}$. The standard mathematical
model of an ideal tesselation also involves symmetric Delone sets Dolbilin
(2015); Dolbilin et al. (2021). The congruent polygonal tiles with straight
edges are guaranteed to be obtained by the Voronoi decomposition of points in
a symmetric Delone set. The underlying principle for tesselations is that if a
2D arrangement of given a set of Voronoi sites is symmetric Delone, then the
Voronoi tessellation results in a unique repeatable space-filling polygon with
straight edges.

##  3 Delone sets with 1-manifolds with boundaries in 2D

Our key observation in this paper the Voronoi sites in symmetric Delone sets
need not be points. They can be 1-manifolds or 1-manifolds with boundaries in
2D space such as line segments, open or closed polylines, and open and closed
curves. In this paper, we only experiment with 1-manifolds with boundaries, in
other words, line segments, polylines, and curve segments. We also assume that
these polylines and curve segments are simple. In other words, these shapes do
not self-intersect with themselves.

We now need to generalize the concept of the original Delone sets with points
to Delone sets with such “simple” 1-manifolds with boundaries. For Delone
sets, we essentially need a well-spaced set of 1-manifolds with boundaries
filling 2D space. Now, let $\mathbf{C}$ denote a set of 1-manifolds with
boundaries with exactly the same geometric shape in 2D Euclidean space,
$\Re^{2}$. Note that for such simple 1-manifolds with boundaries, the distance
to the shape is well-defined. Let $\mathbf{c}\in\mathbf{C}$ denote a simple
1-manifold with boundary and $p\in\Re^{2}$, then let the function
$D(p,\mathbf{c}$ denote the distance of point $p$ to the $\mathbf{c}$. Now,
assume we define a generalized ball of radius $r$ as follows:

$B=\\{p=(x,y)|D(p,\mathbf{c}\leq r.$

For lines, this ball will have the shape of a capsule as shown in Figure 9(a).
Piecewise linear 1-manifolds that consist of multiple lines will produce a
union of these capsules as generalized balls. Using these generalized balls,
we can redefine the terms uniformly discrete and relatively dense for Delone
sets with such “simple” 1-manifolds with boundaries. Similarly, Let
$r_{1}>r_{0}>0$ be two positive numbers. $S$ is uniformly discrete if there
exist an $r_{0}$ such that each ball of radius $r_{0}$ contains at most shape
in $\mathbf{C}$. Similarly, $\mathbf{C}$ is relatively dense if every ball of
radius $r_{1}$ contains at least one shape of $\mathbf{C}$. We call
$\mathbf{C}$ a Delone set if it is both uniformly discrete and relatively
dense. These non-symmetric Delone sets essentially consist of a well-spaced
set of the same type of 1-manifolds with boundaries.

We also need to define symmetric Delone sets for ‘simple” 1-manifolds with
boundaries. For symmetry, the Delone set $\mathbf{C}$ must be invariant with
respect to one of the 17 wallpaper groups. Formally, the Delone set
$\mathbf{C}$ is symmetric if, for every two points
$\mathbf{c_{0}},\mathbf{c_{1}}\in\mathbf{C}$, there exists a rigid motion of
2D space that takes $\mathbf{C}$ to $\mathbf{C}$ and $\mathbf{c_{0}}$ to
$\mathbf{c_{1}}$. We can then obtain congruent polygonal tiles with curved
edges by the Voronoi decomposition of shapes in a symmetric Delone set
$\mathbf{C}$. In other words, the Voronoi tessellations of $\mathbf{C}$ can
result in unique repeatable space-filling polygons with curved edges.

With this in view, the main objectives of our framework are to (1) produce an
arrangement of higher dimensional Voronoi sites based on a given wallpaper
symmetry group, and (2) satisfy the conditions such that the arrangement will
be symmetric Delone and (3) partition the space using Voronoi decomposition.

Note that, we focus only on non-mirror symmetry groups based on isometries of
a cube. The main reason for this restriction is that operations that create
chiral shapes such as mirrors are not practically useful. In multi-mirror
operations, such as pmm and cmm in 2D and pmmm and cmmm in 3D, there is no
difference between using high-dimensional shapes and a single point as Voronoi
sites. Even if a mirror operation is used just once, we end up with planar
(i.e. not-curved) interfaces regardless of the shape of Voronoi sites. We
develop a framework to obtain all elements of symmetries induced by cube
isometries.

###  3.1 Line Segments as Voronoi Sites

Let a line segment be represented by two endpoints $\textbf{p}_{0}$ and
$\textbf{p}_{1}$ as follows

$\textbf{p}=\textbf{p}_{0}(1-t)+\textbf{p}_{0}t$

where $t\in[0,1]$. All lines in a symmetric Delone set can be obtained by
transforming the original two points $\textbf{p}_{0}$ & $\textbf{p}_{1}$, and
creating a line between them. When all the Voronoi sites are line segments,
the boundaries can be either paraboloid curves or lines. Note that the
parabola is already defined as the locus (or collection) of points equidistant
from a given point (the focus) and a given infinite line (the directrix) Jones
(2023) (See Figure 8(a)). If the directrix line is not infinite, i.e. a line
segment, Voronoi decomposition produces only a segment of a parabola (See
Figure 8(b)). Voronoi decomposition with randomly oriented multiple line
segments can produce even more complex boundary that consists of multiple
parabolic curves (see Figure 8(c)) Subramanian et al. (2019). On the other
hand, two mirrored line segments can only produce the original mirror line as
Voronoi boundary as shown in Figure 8(d).

(a) Infinite Line as Directix produces a complete parabola.

(b) A line segment as Directix produces a piecewise parabolic curve.

(c) Two random line segments.

(d) Two mirrored line segments.

Figure 8: These figures show how parabolic Voronoi boundaries are obtained
when using line segments as Voronoi sites. Note that a line segment as
Directix produces a piecewise curve that consists of a parabolic segment along
with two lines. Multiple line segments can produce more complex boundary that
consists of multiple parabolic curves. On the other hand, if two line segments
are mirrored, they only produce the original mirror line as the Voronoi
boundary.

Mirrored shapes (not only line segments) are a major issue to obtain with
curved boundaries. If the symmetry operations include a mirror, we cannot
avoid getting straight boundaries. Figures 5, 6, and 7 show three types of
symmetry operations that can produce straight lines. This problem is something
that is not critical for Delone set-based methods since original Delone sets
use only points to obtain Voronoi decomposition. The results, therefore, will
always be a straight line. It is also important to note that p4m, pmm, p13m,
and p6m are the worst cases since they can produce the same boundaries
regardless of the shapes of Voronoi sites in Figure 7. These four cases are
also undesirable for the original point-based Delone set methods since they
can not provide any control to creating tiles other than four polygons shown
in Figure 7.

##  4 Image Based Voronoi Tessellation with Lines as Voronoi Sites

For Voronoi Tessellation with lines as Voronoi sites can be computed in image
space by observing that the distance functions can be represented as 3D shapes
that consist of cones and triangular prisms as shown in Figure 9(b) Hoff et
al. (1999). Using this property, Voronoi decomposition of multiple lines can
be obtained by creating these shapes in 3D using different colors as shown in
Figure 9(c) and looking at them from the top view as shown in Figure 9(e).
This property makes it easier to implement Voronoi decomposition in graphics
hardware using lines as Voronoi sites. Multiple lines and curves can also be
implemented using union of the distance to a set of line segments.

(a) The shape of a ball that is defined as equidistant to a line.

(b) 3D Shape of the distance to a single line.

(c) Intersection of two 3D distance functions.

(d) The top view of the intersection.

(e) Voronoi partition.

Figure 9: The essence of the method developed by Hoff et al.Hoff et al.
(1999). The top view of the intersection of two 3D shapes that come from 3D
distance functions directly provides Voronoi partition.

###  4.1 Polylines and Curves as Voronoi Sites

Using only single lines as Voronoi sites, we can produce non-convex tiles as
shown in Figures 2, and 3. To construct such complicated tile shapes, one can
use multiple lines as a single Voronoi site. As shown in Figure 4 using
multiple lines as Voronoi sites can produce significantly more complicated
tiles with curved boundaries.

The only constraint with multiple lines is that the Voronoi sites are still
essentially line segments. Therefore, the boundaries can only consist of
parabolic curves. To obtain more general boundary curves, it is necessary to
extend Voronoi sites to curves. Bezier and B-spline curves Bartels et al.
(1995) are appropriate since they can efficiently represent a few number of
line segments through subdivision processes such as DeCasteljau, Chaikin, and
Catmull-Clark subdivision algorithms Boehm and Müller (1999); Chaikin (1974);
Catmull and Clark (1998). The Figures 10, 11, 12, and 13 show examples of
curves as Voronoi sites.

(a) Single Cubic curve with p1 symmetry.

(b) Two Cubic curves with p1 symmetry.

Figure 10: Two examples of P1 symmetry by using Hermitian curves.

(a) Single Cubic curve with pg symmetry.

(b) Two Cubic curves with pg symmetry.

Figure 11: Two examples of Pg symmetry by using Hermitian curves.

(a) Single Cubic curve with p3 symmetry.

(b) Two Cubic curves with p3 symmetry.

Figure 12: Two examples of P3 symmetry by using Hermitian curves.

(a) Single Cubic curve with p6 symmetry.

(b) Two Cubic curves with p6 symmetry.

Figure 13: Two examples of P6 symmetry by using Hermitian curves.

##  5 implementation

For real-time and interactive Voronoi decomposition of 2D space, we have
implemented the image-based Voronoi tessellation of 2D space approach
developed by Hoff et al.Hoff et al. (1999) using graphics hardware. This
approach utilizes the efficiency of the 3D graphics pipeline for discrete
Voronoi diagram computation on a rectangular grid of points (pixels).

For each primitive, we provide the distance function in the form of a 3D
polygonal mesh Hoff et al. (1999). For a point primitive 3D polygonal mesh is
a cone and for a line primitive the polygonal mesh consists of a ”tent” with a
cone at each of the two ends (See Figure 9). Our system consists of an
effective combination of three existing methods: (1) The Voronoi tessellation
of 2D space using graphics hardware; (2) Symmetry operations over the Voronoi
sites; and (3) Coloring and Detecting Boundaries.

These polygonal meshes are rendered at the position of each primitive. This
process defines the Voronoi region for each primitive. The scene is rendered
using a polygon scan conversion and Z-buffer depth comparison engine. The
resulting top view of the scene is an image that is the Voronoi tessellation
of the 2D space.

We have implemented a web-based version and it is available at
https://voronoi.viz.tamu.edu/. In this web-based version, we used the open-
source code repository, EscherSketch Levskaya (2017), to generate the
seventeen planar symmetries for implementing symmetry operations over the
Voronoi sites. We also used the implementation of the drawing interface from
the repository, with minor changes. In addition to the web-based version, we
also implemented this system as a plug-in in Blender, mostly using Blender
functions.

For coloring and detecting boundaries (in both versions), we assign a unique
color to each region of the Voronoi tessellation by computing pseudo-random
values in the HSV color space to avoid two neighboring regions being assigned
the same color. The discrete boundaries of the regions are computed using the
Sobel Operator to detect edges Kanopoulos et al. (1988) in the rendered image
as the change in color would imply a change in the Voronoi region Hoff et al.
(1999).

We have also implemented curves to obtain more complicated boundaries beyond
parabolas. For a simple interface, we use the Hermitian form. Users define
positions and tangents to design curves. We then transform these curves into
the third-degree Bezier polynomials. The advantage of the Bezier form, the
curves can be constructed by using the DeCasteljau subdivision algorithm. The
main advantage of the DeCasteljau subdivision is that it gives a good
approximation of the curve only with 12 or 24 line segments.

To simplify the curve description for users, we also used the Catmull-Rom
formulation that allows to use only positions of points to describe Hermitian
curves by implicitly computing tangent vectors Catmull and Rom (1974); Yuksel
et al. (2011). Since the Catmull-Rom curves are also piecewise cubic curves,
it is also straightforward to transform the Catmull-Rom curves into the third-
degree Bezier polynomials.

##  6 Conclusion and Future Work

In this paper, we present an approach for the development of real-time systems
to interactively design symmetric tiles with curved edges. Our approach is an
extended version of the image-based Voronoi decomposition method Hoff et al.
(1999). We made two extensions: (1) We use higher-order shapes such as curves
as Voronoi sites instead of just points and line segments, and (2) we use
Voronoi sites that are closed under wallpaper symmetries. These two extensions
allow us to create Voronoi tessellations with 2D space-filling tiles with
curved edges. We have developed both a web application and a stand-alone
application that can allow the design of space-filling tiles with curved edges
interactively in real time.

We have also demonstrated that using Voronoi-based methods not all symmetry
operations are useful for creating curved tiles: all symmetries that use
mirror operation produce straight lines that are useless for creating new
tiles. This result is interesting because it suggests that we need to avoid
mirror transformations to produce unusual space-filling tiles in 2D and 3D
using Voronoi tessellations.

This work can be extended for the construction of 2.5D space-filling
structures such as Delaunay Lofts Subramanian et al. (2019), Generalized
Abeille Tiles Akleman et al. (2020), and VoroNoodles Ebert et al. (2003);
Mullins et al. (2022). We currently work on the development of interactive
systems to provide a generalization of all of these methods as layer-by-layer
curves.

###  6.1 Individual Contributions

Haard Panchal implemented the web-based system. Ergun Akleman and Vinayak
Krishnamurthy developed the main approach in this paper. The paper is written
mostly Ergun Akleman with initial text by Haard Panchal. Tolga Talha Yildiz
implemented a Blender-based interactive system. Varda Grover implemented
Catmull-Rom curve interfaces into web-based system. Varda Grover and Tolga
Talha Yildiz are currently implementing 2.5D systems for web-based and stand-
alone systems.

## References

  * Tessellations [2015] Tessellations.  Tesselmaniac!  http://www.tessellations.com/TesselManiac.html, 2015.
  * Amenta [1995] Nina Amenta.  Kali: Drawing symmetric tiles.  http://www.geom.uiuc.edu/apps/kali/about.html, 1995.
  * Akleman et al. [2000] Ergun Akleman, Jianer Chen, and Burak Meric.  Web-based intuitive and effective design of symmetric tiles.  In _Proceedings of the 2000 ACM workshops on Multimedia_ , pages 1–4, 2000.
  * Levskaya [2017] Anselm Levskaya.  Eschersketch.  Github: https://github.com/levskaya/eschersketch, 2017.
  * Delaunay and Sandakova [1961] Boris Nikolaevich Delaunay and Nina Nikolaevna Sandakova.  Theory of stereohedra.  _Trudy Matematicheskogo Instituta imeni VA Steklova_ , 64:28–51, 1961.
  * Schmitt [2016] Moritz W Schmitt.  _On Space Groups and Dirichlet–Voronoi Stereohedra_.  PhD thesis, Berlin: Freien Universität Berlin, 2016.
  * Loeb [1991] Arthur L Loeb.  Space-filling polyhedra.  In _Space Structures_ , pages 127–132. Springer, 1991.
  * Grünbaum and Shephard [1980] Branko Grünbaum and Geoffrey C Shephard.  Tilings with congruent tiles.  _Bulletin of the American Mathematical Society_ , 3(3):951–973, 1980.
  * Du Sautoy [2008] Marcus Du Sautoy.  _Symmetry: A journey into the patterns of nature_.  Harper Collins, 2008.
  * Kaplan [2000] Craig S. Kaplan.  Voronoi diagrams and ornamental design.  In _in The First Annual Symposium of the International Society for the Arts, Mathematics, and Architecture_ , pages 277–283, 2000.
  * Rao [2017] Michaël Rao.  Exhaustive search of convex pentagons which tile the plane.  arXiv preprint arXiv:1708.00274, 2017.
  * Subramanian et al. [2019] Sai Ganesh Subramanian, Mathew Eng, Vinayak R Krishnamurthy, and Ergun Akleman.  Delaunay lofts: A biologically inspired approach for modeling space filling modular structures.  _Computers & Graphics_, 82:73–83, 2019.
  * Akleman et al. [2020] Ergun Akleman, Vinayak R Krishnamurthy, Chia-An Fu, Sai Ganesh Subramanian, Matthew Ebert, Matthew Eng, Courtney Starrett, and Haard Panchal.  Generalized abeille tiles: Topologically interlocked space-filling shapes generated based on fabric symmetries.  _Computers & Graphics_, 89:156–166, 2020.
  * Krishnamurthy et al. [2021] Vinayak R Krishnamurthy, Ergun Akleman, Sai Ganesh Subramanian, Matthew Ebert, Jiaqi Cui, Chia-an Fu, and Courtney Startett.  Geometrically interlocking space-filling tilingbased on fabric weaves.  _Accepted: IEEE Transactions of Visualization and Computer Graphics_ , 2021.
  * Yildiz et al. [2023] Tolga Yildiz, Ergun Akleman, Vinayak Krishnamurthy, and Matthew Ebert.  A modular approach for creation of any bi-axial woven structure with congruent tiles.  _Computers & Graphics_, 2023.
  * Ebert et al. [2003] Matthew Ebert, Ergun Akleman, Vinayak Krishnamurthy, Roman Kulagin, and Yuri Estrin.  Voronoodles: Topological interlocking with helical layered 2-honeycombs.  _Advanced Engineering Materials_ , 2003.
  * Mullins et al. [2022] Cassie Mullins, Matthew Ebert, Ergun Akleman, and Vinayak Krishnamurthy.  Voronoi spaghetti & voronoodles: Topologically interlocked, space-filling, corrugated & congruent tiles.  In _SIGGRAPH Asia 2022 Technical Communications_ , pages 1–4. ACM Siggraph, 2022.
  * Krishnamurthy et al. [2023] Vinayak Raman Krishnamurthy, Ergun Akleman, and Sai Ganesh Subramanian.  Systems and methods for creating space-filling solid models, January 12 2023.  US Patent App. 17/782,583.
  * Delone et al. [1970] Boris Nikolaevich Delone, Nikolai Petrovich Dolbilin, SS Ryškov, and MI Štogrin.  A new construction in the theory of lattice coverings of an n-dimensional space by equal spheres.  _Mathematics of the USSR-Izvestiya_ , 4(2):293, 1970.
  * Delone et al. [1976] Boris Nikolaevich Delone, Nikolai Petrovich Dolbilin, Mikhail Ivanovich Shtogrin, and Ravil V Galiulin.  A local criterion for regularity of a system of points.  _Dokl. Akad. Nauk SSSR_ , 227(1):19–21, 1976.
  * Lagarias [1999a] Jeffrey C Lagarias.  Geometric models for quasicrystals i. delone sets of finite type.  _Discrete & Computational Geometry_, 21(2):161–191, 1999a.
  * Lagarias [1999b] Jeffrey C Lagarias.  Geometric models for quasicrystals ii. local rules under isometries.  _Discrete & Computational Geometry_, 21(3):345–372, 1999b.
  * Lagarias and Pleasants [2003] Jeffrey C Lagarias and Peter AB Pleasants.  Repetitive delone sets and quasicrystals.  _Ergodic Theory and Dynamical Systems_ , 23(3):831–867, 2003.
  * Lagarias [1996] Jeffrey C Lagarias.  Meyer’s concept of quasicrystal and quasiregular sets.  _Communications in mathematical physics_ , 179(2):365–376, 1996.
  * Dolbilin [2015] Nikolay Dolbilin.  Delone sets: local identity and global symmetry.  In _Geometry and Symmetry Conference_ , pages 109–125. Springer, 2015.
  * Dolbilin et al. [2021] Nikolay Dolbilin, Alexey Garber, Undine Leopold, Egon Schulte, and Marjorie Senechal.  On the regularity radius of delone sets in 3d.  _Discrete and Computational Geometry_ , 66(3):996–1024, 2021.
  * Jones [2023] James Jones.  Development of a parabola from the definition.  https://people.richland.edu/james/lecture/m116/conics/paradef.html, 2023.
  * Hoff et al. [1999] Kenneth E. Hoff, John Keyser, Ming Lin, Dinesh Manocha, and Tim Culver.  Fast computation of generalized voronoi diagrams using graphics hardware.  In _Proceedings of the 26th Annual Conference on Computer Graphics and Interactive Techniques_ , SIGGRAPH ’99, page 277–286, USA, 1999. ACM Press/Addison-Wesley Publishing Co.  ISBN 0201485605.  doi: 10.1145/311535.311567.  URL https://doi.org/10.1145/311535.311567.
  * Bartels et al. [1995] Richard H Bartels, John C Beatty, and Brian A Barsky.  _An introduction to splines for use in computer graphics and geometric modeling_.  Morgan Kaufmann, 1995.
  * Boehm and Müller [1999] Wolfgang Boehm and Andreas Müller.  On de casteljau’s algorithm.  _Computer Aided Geometric Design_ , 16(7):587–605, 1999.
  * Chaikin [1974] George Merrill Chaikin.  An algorithm for high-speed curve generation.  _Computer graphics and image processing_ , 3(4):346–349, 1974.
  * Catmull and Clark [1998] Edwin Catmull and James Clark.  Recursively generated b-spline surfaces on arbitrary topological meshes.  In _Seminal graphics: pioneering efforts that shaped the field_ , pages 183–188. ACM, Siggraph, 1998.
  * Kanopoulos et al. [1988] N Kanopoulos, N Vasanthavada, and R L Baker.  Design of an image edge detection filter using the sobel operator.  _IEEE J. Solid-State Circuits_ , 23(2):358–367, April 1988.
  * Catmull and Rom [1974] Edwin Catmull and Raphael Rom.  A class of local interpolating splines.  In _Computer aided geometric design_ , pages 317–326. Elsevier, 1974.
  * Yuksel et al. [2011] Cem Yuksel, Scott Schaefer, and John Keyser.  Parameterization and applications of catmull–rom curves.  _Computer-Aided Design_ , 43(7):747–755, 2011.


The AI community is increasingly focused on merging logic with deep learning to create Neuro-Symbolic (NeSy) paradigms and assist neural approaches with symbolic knowledge.
A significant trend in the literature involves integrating axioms and facts in loss functions by grounding logical symbols with neural networks
and operators with fuzzy semantics.
Logic Tensor Networks (LTN) is one of the leading representatives in this category, known for its simplicity, efficiency, and versatility.
However, it has been previously shown that not all fuzzy operators perform equally when applied in a differentiable setting.
Researchers have proposed several configurations of operators, trading off between effectiveness, numerical stability, and generalization to different formulas.
This paper presents a configuration of fuzzy operators for grounding formulas end-to-end in the logarithm space.
Our goal is to develop a configuration that is more effective than previous proposals, able to handle any formula, and numerically stable.
To achieve this, we propose semantics that are best suited for the logarithm space and introduce novel simplifications and improvements that are crucial for optimization via gradient-descent.
We use LTN as the framework for our experiments, but the conclusions of our work apply to any similar NeSy framework.
Our findings, both formal and empirical, show that the proposed configuration outperforms the state-of-the-art and that each of our modifications is essential in achieving these results.

§ INTRODUCTION

Recently, there has been an increasing interest in combining logic and neural networks in Neuro-Symbolic (NeSy) integrations.
The goal of such systems is often to guide the learning of neural networks using symbolic knowledge, allowing them to reason at a higher level of abstraction.
Much of the recent progress in this area has focused on developing differentiable approaches for knowledge representation and reasoning.

A trend of approaches involves grounding logical symbols using neural networks and relaxing logical operators into continuous operations using fuzzy semantics.
The resulting formulas, such as $\forall x \exists y P(x,y) \lor R(x)$, are associated with a truth degree in the interval $[0,1]$ that represents their level of satisfiability.
This satisfiability can then be derived with respect to the parameters of the neural networks that ground the symbols,
and is incorporated in the loss function of said neural networks to act as an additional supervision when training.
In this study, we conduct experiments and analyses using Logic Tensor Networks (LTN), a well-established framework for differentiable fuzzy logics.

Previous research has highlighted that not all fuzzy operators are appropriate for this type of application.
Different configurations of operators have been proposed in the literature, each with varying degrees of effectiveness, numerical stability, and applicability across different formulas.
However, as of yet, no configuration has met all of these requirements simultaneously.
The goal of this paper is to develop a configuration of operators that is superior to previous proposals and capable of handling any formula.
To achieve this, we propose operators in the logarithm space, which is known to address certain issues.
We build upon existing findings and introduce novel improvements that are crucial for optimization through gradient-descent.
We call our new solution , and release it on the official github repository for LTN.

The foundation of the LTN framework is explained in Section <ref>.
In Section <ref>, we provide an in-depth examination of the semantics in the logarithm space, with limitations in section <ref>.
Our main contribution is in Section <ref>, which includes all the key simplifications and computational techniques, along with their formal justifications, that improve the derivability of the framework.

In Section <ref>, we experimentally confirm that our proposition surpasses state-of-the-art configurations and, using ablation studies, that each of our modifications plays a critical role in achieving these results.

Our findings are expected to help all differentiable frameworks that rely on fuzzy semantics.

§ BACKGROUND ON LTN

§.§ Real Logic concepts

LTN is built on Real Logic, a first-order language that allows to specify relational knowledge about the world.
For example, the formula $\Friend(\alice,\bob)$ states that $\alice$ is a friend of $\bob$,
and the formula $\forall u \forall v (\Friend(u,v) \imp \Friend(v,u))$ states that $\Friend$ is a symmetric relation,
where $u$ and $v$ are variables, $\alice$ and $\bob$ are individuals, and $\Friend$ is a predicate.

In Real Logic, a grounding $\G_\theta$ associates mathematical, real-valued semantics to every logical symbol depending on a set of parameters $\btheta$.
Individuals are grounded with vectors of real values.
Often, the vectors come from real-world features and data.
A variable is grounded with a finite batch of individuals from a domain.
Finally, relations are grounded using mathematical functions (generally, neural networks) that map to the truth domain $[0,1]$.

Complex formulas are constructed using the usual logical connectives and quantifiers $\land$, $\lor$, $\imp$, $\lnot$, $\forall$, $\exists$.
The connectives are grounded using t-norms fuzzy logic:
$\land$ is grounded using a t-norm $\tnorm$,
$\lor$ using a t-conorm $\snorm$,
$\imp$ using a fuzzy implication $\fzimp$,
and $\lnot$ using a negation $\fzneg$.
The quantifiers are grounded using aggregators $\Agg^\forall$ and $\Agg^\exists$.

Examples of fuzzy operators are
the standard negation $\Ns(x)=1-x$,
the product t-norm $\Tprod(x,y)=xy$,
and its dual t-conorm $\Sprod(x,y)=x+y-xy$.
$\Ns$ is inspired by the negation of a probability.
$\Tprod$ is inspired by the intersection probability of two independent events.
$\Sprod$ is the dual t-conorm derived from the other two operators using De Morgan's laws.

For brevity, let us denote $\mathbf{x} = (x_1, \dots, x_n)$ as a vector of $n$ values.
An example of a universal aggregator is $\Atprod(\mathbf{x}) = \prod_{i=1}^{n} x_i$, which is equivalent to the conjunction of $n$ events.

Notice that all the operators function within the usual interval $[0,1]$.

In LTN, the parameters $\btheta$ are learned using maximal satisfiability of a knowledgebase $\K$.
The satisfaction of a formula $\phi$ is its evaluation $\Gt(\phi)$, which returns a truth-value in $[0,1]$.
Let $\K$ define a collection of formulas.
The satisfaction of $\K$ is defined as the aggregation of the satisfactions of each $\phi \in \K$.
The result depends on the choice of aggregate operator, denoted by $\Sat$ (typically, the same operator as the universal aggregator).

The optimal set of parameters maximizes the objective function $  \btheta^\ast = \argmax_{\btheta}\ \Sat_{\phi\in\K}\Gt(\phi)$.
The following loss function is used to find that objective via gradient descent:

\begin{equation}
\mathcal{L}_{\LTN}(\Gt, \K) = - \Sat_{\phi \in \K} \Gt(\phi)
\end{equation}

For more intuition, we give a concrete example in Appendix <ref>.

§.§ Appropriate Operators for Gradient Optimization

The ability to find an optimum satisfying a formula greatly depends on the choice of operators that ground the logical connectives.
[21] demonstrate that some fuzzy logic operators are unsuitable in a differentiable setting.

For example, the Łukasiewicz t-norm $\Tluk(x,y) = \max(x+y-1,0)$ has vanishing gradients when $x+y-1<0$.

The authors show that the Product Real Logic configuration is the most suitable for grounding the logical connectives.
It uses the product t-norm, its dual t-conorm and the standard negation.
For the universal aggregator, it avoids the potential underflow issues of multiplying many small numbers together by working with the log-product $(\log \circ \Atprod)(\mathbf{x}) = \sum_{i=1}^{n} \log(x_i)$.
Because this configuration mixes operators in the usual and logarithm spaces, it has limitations in expressivity and cannot handle certain formulas (e.g. $(\forall u P(u)) \lor (\forall v Q(v))$).
We discuss this further in Appendix <ref>.

In this paper, we aim to explore a configuration that
1) performs better than Product Real Logic,
2) can handle any formula,
and 3) is numerically stable.

§ INTRODUCING

We present , a specification of LTN with end-to-end semantics in the logarithm space.
Section <ref> introduces operators that can manipulate appropriately log truth degrees.
Section <ref> shows how to modify these operators to perform well in a differentiable setting.

§.§ Semantics

We employ the product t-norm $\Tprod(x,y) = xy$ and the maximum t-conorm $\Smax(x,y) = \max(x,y)$.
These operators are known to simplify easily in the logarithm space and are commonly used in the log probability literature.
We also use the standard negation operator $\Ns(1-x)$.
Implications are replaced using the material implication rule $(\phi \imp \psi) \equiv (\lnot\phi \lor \psi)$,
which means we rewrite every implication using $\fzimp(x,y) = \snorm(\fzneg(x), y)$.

The universal aggregator is defined as the conjunction of $n$ events $\Atprod(\mathbf{x}) = \prod_{i=1}^{n} x_i$,
and the existential aggregator is defined as the disjunction of $n$ events $\Asmax(\mathbf{x}) = \max_{i=1}^n (x_i)$.

§.§.§ Logarithm space

We denote $(\log\circ\Gt)(\phi)$ as the log-grounding of a formula and $(\log \circ\ \mathrm{FuzzyOp})$ as the log-grounding of an operator.
Note that maximizing the log grounding of a formula is equivalent to maximizing its grounding as logarithms are monotone increasing functions.
The log-grounding of $\Tprod$, $\Smax$, and their generalizations in aggregators, simplify easily.

\begin{align}
\label{eq:log_semantics_conj}    (\log\circ \Tprod) (x,y) & = \log(x) + \log(y) \\
\label{eq:log_semantics_disj}    (\log\circ \Smax) (x,y) & = \max(\log(x),\log(y))\\
\label{eq:log_semantics_forall} (\log \circ \Atprod)(\mathbf{x}) & = \sum_{i=1}^n \log(x_i)\\
\label{eq:log_semantics_exst}(\log \circ \Asmax)(\mathbf{x}) & = \max_{i=1}^n (\log(x_i))
\end{align}

Expressing $\log\circ\Ns(x)$ as a function with a logarithmic input
requires the computation of an exponent and a logarithm.
This means that the operator cannot easily take an input in the logarithm space,
for example, in $\lnot(A \land B)$.
To overcome this, we write formulas in negative normal form (NNF).
A formula is in NNF
when the scope of each negation operator only applies to atoms (predicates), not to complex formulas,
and when the formula does not contain any implication or equivalence symbols.
For example, if $A$ and $B$ are two atoms, $\lnot A \land \lnot B$ is in NNF but $\lnot(A \lor B)$ is not.

§.§ Optimizing in

§.§.§ Numerical stability of log negations

Let $f(x)$ be the output of a neural predicate in the interval $[0,1]$ depending on a mathematical variable $x$.
Converting a value to the logarithm space is a risky operation in a computational graph, as both $\log(f(x))$ and $\frac{\partial \log(f(x))}{\partial x}= \frac{1}{f(x)} \frac{\partial f(x)}{\partial x}$
can cause overflow errors when $f(x)$ tends to $0$.

In NeSy AI, predicates are typically grounded using a final sigmoid or softmax layer to normalize outputs in $[0,1]$.
Fortunately, the computation and differentiation of the logarithm of a sigmoid or softmax simplifies to a stable expression (refer to Appendix <ref>).
For this reason, most frameworks for automatic differentiation, such as TensorFlow or PyTorch,
offer built-in and all-in-one layer implementations of the log sigmoid and log softmax functions.

These should be used when log-grounding a predicate to avoid unstable gradients.

However, the same issue arises when log-grounding the negation of a predicate, $\log(1-f(x))$,
Fortunately, we show how to reformulate the log-negation of a sigmoid or softmax predicate to numerically stable expressions.

The log-negation of a sigmoid function $S(x)=\frac{1}{1+e^{-x}}$, $x \in \R$, simplifies as

\begin{equation}
\label{eq:log_not_sigmoid}
(\log\circ\Ns)(S(x)) = \log(S(x)) - x
\end{equation}

Proof in Appendix <ref>.

The log-negation of a softmax function $\softmax(\mathbf{z})_i = \frac{e^{z_i}}{\sum_{j=1}^{K} e^{z_j}}$,
where $\mathbf{z} =(z_{1},\dotsc ,z_{K})\in \mathbb {R} ^{K}$ is a vector of $K$ real values,
and $i=1,\dots,K$, simplifies as

\begin{equation}
\label{eq:log_not_softmax}
(\log \circ \Ns) (\softmax(\mathbf{z})_i) = \log(\softmax(\mathbf{z})_i) + \log (\sum_{\substack{j=1\\j\neq i}}^K e^{z_j}) - z_i
\end{equation}

Proof in Appendix <ref>.

These two proposed reformulations have numerically stable implementations.
The first uses the logarithm of a sigmoid and a linear term.
The second uses the logarithm of a softmax function, a linear term, and a logarithm of a sum of exponentials, also known as LogSumExp.
LogSumExp also has a stable implementation and its derivative is a softmax function.

Below, we briefly show the stability advantage of our reformulation for sigmoid by comparing its output with a naive definition $\log(1-\sigmoid(x))$.
The results are obtained in TensorFlow with float32 precision.
The same can be reproduced with the softmax reformulation.

§.§.§ Relaxation of the disjunctions

The maximum operator in equations (<ref>) and (<ref>) is unsuitable in a differentiable setting as it has single-passing gradients.
This means that it only propagates gradients to one input at a time, the one with the highest value.
Intuitively, let the formula $\exists x \ P(x)$ be a constraint used to optimize a neural predictor $P$.
If several individuals in the batch $x$ tend to verify $P(x)$, $\max$ will have non-zero gradients for only one of them.
This can be inefficient in practice as it will push the predictor to overfit that single individual in $x$ and ignore the others. Also, it is particularly sensitive to initial conditions.

A common solution is to use a smooth approximation of the maximum operator.
A popular candidate in the logarithm space is the LogSumExp ($\LSE$) operator, defined as:

\begin{align}
\LSE(\mathbf{x} \mid \alpha, C) = \frac{1}{\alpha} \left( C + \log(\sum_{i=1}^n e^{\alpha x_i - C})\right)\\
\end{align}

$C=\max(\alpha \mathbf{x})$ is a constant that does not change the result of the expression but prevents overflow errors in the exponential terms.
$\alpha$ is a hyperparameter that scales the bounds of $\LSE$ according to the following inequality:

\begin{equation}
\label{eq:bounds_LSE_traditional}
\max(\mathbf{x}) \leq \LSE(\mathbf{x} \mid \alpha, C) \leq \max(\mathbf{x}) + \frac{\log(n)}{\alpha}
\end{equation}

Here, we identify an issue in that LogSumExp approaches the maximum value via a higher bound.
This is problematic, as truth degrees are bound to the interval $[0,1]$,
and log truth degrees should be bound in the interval $\left[-\infty,0\right]$.
[The edge case $\log(0) = -\infty$ can be avoided by add a small real value $\epsilon>0$ to zero truth degrees.
However, this is rarely a problem in practice as sigmoid and softmax layers output values in $]0,1[$.]
With $\LSE$, the output of a log-disjunction can exceed these bounds and become non-negative.

To address this issue, we propose the use of a LogMeanExp operator $\LME$:

\begin{equation}
\label{eq:LME}
\LME(\mathbf{x} \mid \alpha, C) = \frac{1}{\alpha} \left( C + \log(\frac{\sum_{i=1}^n e^{\alpha x_i-C}}{n} ) \right)
\end{equation}

It approaches the maximum operator from below values (proof in Appendix <ref>):

\begin{equation}
\label{eq:bounds_LSE_lower}
\max(\mathbf{x}) - \frac{\log(n)}{\alpha} \leq \LME(\mathbf{x} \mid \alpha, C) \leq \max(\mathbf{x})
\end{equation}

This operator is numerically stable, well-bounded, and suitable for derivation.
We use it to ground disjunctions and existential quantifications in .
For best practice, the parameter $\alpha$ that scales the smooth maximum should be scheduled over time to balance exploration and exploitation [3].

§.§.§ Batch-size invariance for the universal aggregation

We have improved the derivability
of negations, disjunctions, and existential quantifiers in the logarithm space.
Here, we identify an issue with the universal quantification.

Consider a knowledge base with two rules $\phi_1 = \forall u P(u)$ and $\phi_2 = \exists v Q(v)$.
Let $\G(u) = [x_1, \dots, x_m]$ and $\G(v) = [y_1, \dots, y_n]$ be two batches of individuals.
Let us develop the groundings of the rules:

\begin{align}
\label{eq:batch_inv_forall}     (\log\circ\Gt)(\phi_1) &
% = (\log\circ\G)(\forall u P(u))
= \sum_{i=1}^m (\log\circ\Gt)(P)(x_i) \\
\label{eq:batch_inv_exists}     (\log\circ\Gt)(\phi_2) &
% = (\log\circ\G)(\exists v Q(v))
= \max_{i=1}^n (\log\circ\Gt)(Q)(y_i)
\end{align}

In the loss $\mathcal{L} = - (\log\circ\Gt)(\phi_1 \land \phi_2) = - (\log\circ\Gt)(\phi_1) - (\log\circ\Gt)(\phi_2)$,
the optimization will tend to overfit the rule with the universal quantifier and ignore the existential rule for large batch sizes.
This is due to the fact that the first sums $m$ log truth degrees, whereas the second only takes one log truth degree as a maximum.
In terms of differentiability,
the universal rule weights more on the gradient updates,
as $\sum_{i=1}^m \frac{\partial \mathcal{L}}{\partial \log\circ\G(P)(x_i)} = \sum_{i=1}^m -1 = -m$,
whereas $\sum_{i=1}^n\partialderivative{\mathcal{L}}{\log\circ\G(Q)(y_i)} = \sum_{i=1}^n - \delta_{ij} = -1$ given $j=\argmax_{j=1}^n \G(Q)(y_j)$.
The problem remains with the smooth maximum LogMeanExp, whose gradients are a softmax function summing to 1 as well.

This weighing problem can also arise when comparing two universal quantifiers, such as $\forall u P(u)$ and $\forall v Q(v)$.
If the batches for $u$ and $v$ have varying sizes, the optimization algorithm will tend to overfit the rule that has more examples of individuals and ignore the other.

To solve this problem, we propose to use a mean instead of a sum as a weighting scheme to balance universal quantifiers:

\begin{equation}
\label{eq:forall_invariant}
(\log\circ\G)(\forall u P(u)) = \sum_{i=1}^m \frac{(\log\circ\G)(P)(x_i)}{m}
\end{equation}

By averaging log truth degrees instead of summing them, we obtain a batch-size invariant aggregator.
The weight of the gradients becomes $\sum_{i=1}^m \frac{\partial \mathcal{L}}{\partial \log\circ\G(P)(x_i)} = -1$.
This ensures every formula weighs equally in the loss function.

In the normal space, the universal quantifiers then correspond to geometric means instead of products.
This trick alters the objective and search space of the task, but we find that it is crucial to approach good solutions in our experiments.

§.§.§ Summary

By implementing all the aforementioned modifications,
we arrive at the log-grounding routine presented in Algorithm <ref>.

Compute $(\log\circ\Gt)(\phi)$ for derivability

[1]
Step 1  Rewrite $\phi$ in negative normal form
Step 2  Log-ground predicates and their negations
Use the log-negation simplifications with sigmoid/softmax layers (Equations <ref> and <ref>)

Step 3  Compute connectives in the logarithm space
$\land$ becomes $+$
$\lor$ becomes $\LME$
$\forall$ becomes $\mean$
$\exists$ becomes $\LME$

§ LIMITATIONS

§.§ De Morgan's Laws and NNF

The Negative Normal Form (NNF) does not preserve equivalence with the logarithmic semantics introduced in Section <ref>.
To transform a formula in NNF, one must push the negations in front of atoms using De Morgan's laws, but the laws do not hold given that $\Tprod(x,y)=xy$ and $\Smax(x,y) = \max(x,y)$ are not fuzzy dual operators.
However, we can prove the following De Morgan's inequalities:

Let $P$ and $Q$ be two formulas. We can show that

\begin{align}
\G(\lnot(P \land Q)) & \geq \G(\lnot P \lor \lnot Q) \\
\G(\lnot(P \lor Q)) & \geq \G(\lnot P \land \lnot Q) \\
\G(\lnot(\forall u P(u))) & \geq \G(\exists u \lnot P(u)) \\
\G(\lnot(\exists u P(u))) & \geq \G(\forall u \lnot P(u))
\end{align}

Proof in Appendix <ref>.

These results stem from the fact that $\max(x,y)$ is a lower-bound to other t-conorms including the dual product t-conorm.
In Appendix <ref>, we analyze the tightness of these bounds.
Given that NNF is obtained by repeatedly applying De Morgan's laws, we can infer the following property:

Let $\phi$ be any formula and $\phi'$ be a NNF formula derived syntactically from $\phi$ using De Morgan's laws.
Then, $\Gt(\phi') \leq \Gt(\phi)$, meaning that the satisfaction of the NNF formula $\phi'$ is a lower bound for the satisfaction of the formula $\phi$.

This is particularly useful as, if we convert a formula into NNF and find a parametric grounding that satisfies it, we know that the original formula is at least as satisfied.

§ EXPERIMENTS

§.§ Task 1: Clustering

The first experiment is a clustering problem based on the gene expression cancer RNA-Seq benchmark from the UCI ML datasets repository [9].
The dataset has 801 samples of 20531 features, which we reduce to 16 features using PCA.
The task is to divide samples into five clusters, which roughly correspond with five ground truth cancer types.
We train a neural predictor $C(x,c)$ that returns the belief of a point $x$ belonging in a given cluster $c$ using these three constraints:

\begin{gather}
\label{eq:clustering1}    \forall x \exists c\ C(x,c)\\
\label{eq:clustering2}    \forall c \exists x\ C(x,c)\\
\label{eq:clustering3}    \forall (c,x,y: \lvert x-y \rvert < \mathrm{th}_{2.5})\ C(x,c) \imp C(y,c)
\end{gather}

$x$ and $y$ are grounded with the batch of 801 points.
$c$ is a variable that ranges over five cluster ids.
$C(x,c)$ outputs beliefs using a softmax output layer that ensures mutual exclusivity of clusters.
(<ref>) uses the concept of "guarded quantification" introduced by [3].
It means that the quantification only retains the individuals verifying the condition $\lvert x-y \rvert < \mathrm{th}_{2.5}$,
where $\mathrm{th}_{2.5}$ is the $2.5$-th percentile of the euclidean distances between all pairs of points.
Intuitively, the constraint states that for any pair of points that are very close, if one belongs to a cluster, the other must belong in the same cluster.

We use this task, inspired by the toy example from [3] and extended on real-world data, as it is one of the rare NeSy tasks with existential clauses ranging over many individuals (here, 801 individuals).

§.§ Task 2: MNISTAdd

The second problem is the MNISTAdd task from [13].
We experiment on the 2-digits number variant of the task.
In this problem, we learn to recognize the sum of two numbers of two digits using only the result of the sum as a training label.
For example, a dataset sample would be $([\mnistthree,\mnisteight], [\mnistnine,\mnisttwo], 130)$.

The task is modeled using a digit classifier $\isdigit(\mnisteight,d)$ which predicts beliefs for the MNIST image being the digit $d=0\dots 9$.
Whereas we only provide labels for the final result of the addition, NeSy methods use prior knowledge about intermediate labels (possible digits used in the addition) to propagate ground truth information to the digit classifier.
Given the MNIST images and sum result $([x_1,x_2],[y_1, y_2],n)$, we use the LTN constraint from [3]:

\begin{align}
& \exists d_1,d_2,d_3,d_4 : 10d_1+d_2+10d_3+d_4=n \\
& \quad (\isdigit(x_1,d_1) \land \isdigit(x_2,d_2) \land \isdigit(y_1,d_3) \land \isdigit(y_2,d_4)) \notag
\end{align}

The loss signal is a universal aggregation of the constraint over minibatches of labeled examples.
We use the same neural network for $\isdigit$ as [13].
This is a basic experiment with a single training constraint.
We assess it because many probabilistic NeSy methods use it as a standard for comparison.

§.§ Task 3: Semantic PASCAL-Part

The third experiment is a semi-supervised semantic image interpretation task on the semantic PASCAL-Part dataset from [8].
The goal is to train a type classifier $\istype(x,\bottle)$, $\istype(x,\captype)$, etc., that predicts the type of an object within a bounding box $x$,
and to train a relation predictor $\ispartof(x,y)$ that determines if one bounding box $x$ is part of another bounding box $y$.
An example of such bounding boxes is presented in Figure <ref>.

Example of bounding boxes from PASCAL-Part.

Training is guided by three constraints based on ground truth examples,
one for labeled type examples, one for pairs of positive examples for $\ispartof$, and one for pairs of negative examples.
The ground truth labels are made available for only $5\%$ of the training data.
However, training is carried on the unlabeled data using mereological constraints that relate to the types and their meanings, for example:

\begin{gather}
\forall x,y \ \istype(x,\bottle) \land \ispartof(y,x) \imp (\istype(y,\captype) \lor \istype(x,\body)) \label{eq:constraint_pascalpart} \\
\forall x,y \ \istype(x,\captype) \land \ispartof(x,y) \imp \istype(y,\bottle)
\end{gather}

[8] grounded the bounding boxes using predictions produced by an object detector trained on PASCAL-Part.
This means that LTN was only used to correct the predictions of the detector.
We increase the difficulty of the task by implementing the bounding boxes with a latent vector of 1024 features output by a pre-trained FasterRCNN backbone.
That is, LTN has to learn all the final layers of the object detector and its specialization on PASCAL-Part.

We release our version of the dataset on <https://github.com/sbadredd/semantic-pascal-part>.
More details on the experiment are available in Appendix <ref>.

Out of the three tasks, this is by far the largest with a total 59 object types and 60 corresponding constraints.
It showcases the power of LTN and its capability to simply integrate many constraints in a loss function.
We evaluate the type classification using balanced accuracy, $\ispartof(x,y)$ using the area under precision-recall curves,
and the semantic interpretation by reporting the number of false positives that violate the mereological constraints – for example, a bottle is predicted to be part of a cap.

§.§ Baselines

We compare  with the following baselines.
Note that we skip the ablation of the log-negation simplifications, as their numerical practicability is already illustrated in Section <ref>.

Product Real Logic was identified by [21] as the best performing operator semantics for differentiable fuzzy logics.
It uses the product t-norm $\Tprod(x,y)=xy$ and its dual t-conorm $\Sprod(x,y)=x+y-xy$.
The universal quantifier uses the log-product aggregator $(\log \circ \Atprod)(\mathbf{x}) = \sum_{i=1}^{n} \log(x_i)$
and the existential quantifier uses a smooth maximum.
Its combination of operators both in the usual and logarithm space makes it difficult to handle certain formulas.
We discuss this issue further in appendix <ref>.

Stable Product Real Logic [3] is a modification of that uses a smooth minimum for the universal aggregator,
such that all operators perform in the usual space.

A limitation of the smooth minimum is that it depends on a smoothing hyperparameter $p$, which we show to greatly influence the results.

the configuration introduced in this paper, performing fully in the logarithm space.

an ablation of logLTN using a sum instead of a batch-size invariant mean for universal aggregations. See Section <ref>

an ablation logLTN that uses a non-relaxed maximum operator for existential aggregations. See Section <ref>.

an ablation of logLTN that uses a traditional LogSumExp operator for existential aggregations instead of LogMeanExp. See Section <ref>.

§ RESULTS

The code for our experiments is available at <https://github.com/sbadredd/logltn-experiments>.
We perform experiments with runs of 1000 training steps for the clustering problem, runs of 20 epochs for MNISTAdd with two different dataset sizes, and runs of 1000 training steps for Semantic PASCAL-Part.
The results are summarized in Tables <ref> and <ref>.
Because Semantic PASCAL-Part is more computationally demanding, we conducted ablation studies only on the two first experiments.
For additional implementation details, please refer to Appendix <ref> covering training configuration and baseline hyperparameters.
Across all metrics,  consistently achieved the best or second best performance.

performs poorly in the clustering task.
In a qualitative analysis of the cluster assignments (Figure <ref>), we observe that disregards constraint (<ref>) stating that each cluster contains at least one point.
This is due to the batch-variant log-product.
As that constraint is quantified over five clusters ($\forall c$), it has relatively less weight compared to the other ones aggregated over all points.
Also, in Semantic PASCAL-Part, while exhibits good results,  still demonstrates superior performance by avoiding on average 34% more mereological violations and having less deviation across all metrics than .

We tested with the smooth minimum parameter $p=2$ and $p=6$.
When $p=2$, the smooth minimum is less strict and corresponds to a Mean-Squared Error (MSE) aggregator.
This leads to low accuracy in the PASCAL-Part problem, as the constraint aggregator focuses on satisfying the 57 logical constraints rather than the three ground truth constraints which it treats as "outliers".
The result is a predictor classifying all objects into barely constrained types (e.g. background) and all $\ispartof(x,y)$ as false negatives in order to reach low mereological violations.
With $p=6$, the aggregator is more strict but can overfit outliers and exhibit instability in other experiments.
Despite its name, we find to be too dependent on the hyperparameter $p$ and unstable.

Regarding the ablations, generally performed poorly due to inadequate gradient propagation.
showed similar performance to  overall, except for deviating results in MNISTAdd, possibly due to the unbounded maximum breaking at an edge case.
Also, even on well-performing problems, we find that grounding the knowledgebase with $\logLTNLSE$ can lead to log truth degrees reaching values as high as $(\log\circ\G)(\K)=5$.
That corresponds to a truth degree of approximately $(\log\circ\G)(\K)=150$.
Since fuzzy truth degrees should be within the range of $[0,1]$, the higher values generated by make it unusable in many cases, making  a more suitable option.

Finally, exhibits behavior akin to in clustering due to its batch-variant aggregator.
In MNISTAdd, outperforms , but the only difference between the two baselines is a constant factor in the loss function due to taking a mean over the minibatch of samples instead of a sum.
We assume that scaling the learning rate accordingly would yield comparable results with .

In Table <ref>, we compare our MNISTAdd results with those reported by popular probabilistic frameworks [13, 23, 15].
In their study, [3] showed that LTN managed to solve the MNISTAdd problem, but the outcomes varied significantly due to instability during initialization.
We show that by training with  in the logarithm space, we resolved this issue and achieved standard state-of-the-art performance in the task.

2D PCA plots of typical cluster assignments by each baseline

Clustering     2cMNISTAdd

1,500 samples     15,000 samples

$0.37 \pm 0.18$     $ 88.39 \pm 0.97 ^*$     $ 95.39 \pm 0.33 $

(p=2)     $0.71 \pm 0.10$     $ 62.77 \pm 37.16 $     $ 80.34 \pm 32.82 $

(p=6)     $0.62 \pm 0.05$     $ 30.51 \pm 34.19 $     $ 95.12 \pm 0.62 $

$0.71 \pm 0.06 ^*$     $ 88.29 \pm 0.78 $     $ 95.61 \pm 0.51 ^* $

$ {\bf 0.72\pm 0.05 } $      $ 75.35 \pm 29.96 $     $ 95.27 \pm 0.53 $

$0.35\pm 0.17$     $ 48.39 \pm 36.80 $     $ 80.81 \pm 32.64 $

$0.36\pm 0.22$     $ {\bf 88.49 \pm 0.88}  $     $ {\bf 95.62 \pm 0.42}  $

Results on Clustering (Adjusted Rand Index scores – averaged on 10 runs) and MNISTAdd (test accuracy – averaged on 5 runs).

1,500 samples     15,000 samples

DeepProbLog     $ 87.21 \pm 1.92^* $     $ 95.16 \pm 1.70 $

DeepStochLog     NA     $ {\bf 96.40 \pm 0.10} $

NeuPSL     $ 87.05 \pm 1.48 $     $ 93.91 \pm 0.37 $

$ {\bf 88.29 \pm 0.78} $     $ 95.61 \pm 0.51^* $

Reported test accuracy on MNISTAdd by probabilistic baselines.
DeepStochLog does not report results on the training set size 1,500.

PartOf AUC     Type Accuracy     # Mereological Violations

${\bf 71.38 \pm 10.38} $     $53.40 \pm 1.81^*$     $16,312.8 \pm 4,180.2 $

(p=2)     $18.62 \pm 28.16$     $1.90 \pm 0.46$     ${\bf 1,196.0 \pm 2,674.3}$

(p=6)     $64.85 \pm 3.05$     $34.53 \pm 1.62$     $24,324.5 \pm 14,052.3$

$69.57 \pm 4.33 ^*$     ${\bf 55.54 \pm 1.00}$     $10,708.2 \pm 1,455.9 ^*$

Test results on Semantic PASCAL-Part averaged on 5 runs.

§ RELATED WORK

The field of combining logic and neural networks in NeSy integrations is gaining interest, as outlined by [10]. For an overview of the approaches and challenges, see [12]. To understand the prevalence of these systems, refer to [17].

A family of approaches converts logical connectives into differentiable operations using fuzzy semantics.
Systems that employ this approach include
LTN [18, 3],
KALE [11],
SBR [7],
and LRNN [19]

among others.
Unlike probabilistic logics [13, 23], fuzzy approaches change logic semantics and are less common in proof reasoning.
Nevertheless, fuzzy frameworks excel in knowledge-aided learning and offer simplicity compared to probabilistic methods which must often solve the exponentially complex model counting problem.

Fuzzy frameworks have been used in a wide range of applications and fields in recent years.
These include but are not limited to semantic image interpretation [8],
natural language processing [4],
reinforcement learning [2],
query answering over knowledge graphs [1, 6],
or open-world reasoning [22].

Our paper aligns with the research stream of [21] as it strives to improve the performance of all these related works by providing mathematical and computational cues for fuzzy semantics.

§ CONCLUSIONS

Many NeSy approaches rely on fuzzy operator semantics to ground knowledge in loss functions.
However, it is clear that not all semantics are suitable for gradient descent optimization algorithms.
In this paper, we propose a set of semantics that can be used to train logic end-to-end in the logarithm space.
We demonstrate that the proposed configuration outperforms semantics previously considered state-of-the-art in such NeSy systems.

We propose the solution, which we refer to as logLTN, as an additional set of semantics for LTN and implement it in the repository of the framework.
Each of our findings can also be applied separately to any framework that works with logic in the logarithm space.
In summary, our recommendations for such systems include computing log-negations using Equations (<ref>) and (<ref>), relaxing disjunctions using Equation (<ref>), and making universal quantifications batch size-invariant using Equation (<ref>).
Our research is expected to improve the performance of all NeSy approaches that rely on fuzzy operator semantics.

[1]
Arakelyan, Erik, Daza, Daniel, Minervini, Pasquale, & Cochez, Michael. 2021.
Complex Query Answering with Neural Link Predictors.
In: 9th International Conference on Learning
Representations, ICLR 2021, Virtual Event, Austria, May 3-7, 2021.

[2]
Badreddine, Samy, & Spranger, Michael. 2019.
Injecting Prior Knowledge for Transfer Learning into Reinforcement
Learning Algorithms using Logic Tensor Networks.
In: Doran, Derek, d'Avila Garcez, Artur S., &
Lécué, Freddy (eds), Proceedings of the 2019 International
Workshop on Neural-Symbolic Learning and Reasoning (NeSy 2019), Annual
workshop of the Neural-Symbolic Learning and Reasoning Association, Macao,
China, August 12, 2019.

[3]
Badreddine, Samy, d'Avila Garcez, Artur, Serafini, Luciano, & Spranger,
Michael. 2022.
Logic Tensor Networks.
Artificial Intelligence, 303(Feb.), 103649.

[4]
Bianchi, Federico, Palmonari, Matteo, Hitzler, Pascal, & Serafini, Luciano.
Complementing Logical Reasoning with Sub-symbolic
Pages  161–170 of: Fodor, Paul, Montali, Marco, Calvanese,
Diego, & Roman, Dumitru (eds), Rules and Reasoning.
Lecture Notes in Computer Science.
Cham: Springer International Publishing.

[5]
Chen, Xianjie, Mottaghi, Roozbeh, Liu, Xiaobai, Fidler, Sanja, Urtasun, Raquel,
& Yuille, Alan L. 2014.
Detect What You Can: Detecting and Representing Objects Using
Holistic Models and Body Parts.
Pages  1979–1986 of: 2014 IEEE Conference on Computer
Vision and Pattern Recognition, CVPR 2014, Columbus, OH, USA, June 23-28,
IEEE Computer Society.

[6]
Chen, Xuelu, Hu, Ziniu, & Sun, Yizhou. 2022.
Fuzzy Logic Based Logical Query Answering on Knowledge Graphs.
Proceedings of the AAAI Conference on Artificial Intelligence,
36(4), 3939–3948.

[7]
Diligenti, Michelangelo, Roychowdhury, Soumali, & Gori, Marco. 2017.
Integrating Prior Knowledge into Deep Learning.
Pages  920–923 of: 2017 16th IEEE International
Conference on Machine Learning and Applications (ICMLA).
Cancun, Mexico: IEEE.

[8]
Donadello, Ivan, Serafini, Luciano, & d'Avila Garcez, Artur S. 2017.
Logic Tensor Networks for Semantic Image Interpretation.
Pages  1596–1602 of: Sierra, Carles (ed), Proceedings of
the Twenty-Sixth International Joint Conference on Artificial Intelligence,
IJCAI 2017, Melbourne, Australia, August 19-25, 2017.

[9]
Dua, Dheeru, & Graff, Casey. 2017.
UCI Machine Learning Repository.

[10]
Garcez, Artur d'Avila, & Lamb, Luis C. 2020 (Dec.).
Neurosymbolic AI: The 3rd Wave.
arXiv:2012.05876 [cs].

[11]
Guo, Shu, Wang, Quan, Wang, Lihong, Wang, Bin, & Guo, Li. 2016.
Jointly Embedding Knowledge Graphs and Logical Rules.
Pages  192–202 of: Su, Jian, Carreras, Xavier, & Duh, Kevin
(eds), Proceedings of the 2016 Conference on Empirical Methods in
Natural Language Processing, EMNLP 2016, Austin, Texas, USA, November 1-4,
The Association for Computational Linguistics.

[12]
Hitzler, P., & Sarker, M.K. 2022.
Neuro-symbolic Artificial Intelligence: The State of the
Frontiers in artificial intelligence and applications.
IOS Press.

[13]
Manhaeve, Robin, Dumancic, Sebastijan, Kimmig, Angelika, Demeester, Thomas, &
De Raedt, Luc. 2018.
DeepProbLog: Neural Probabilistic Logic Programming.
In: Bengio, S., Wallach, H., Larochelle, H., Grauman, K.,
Cesa-Bianchi, N., & Garnett, R. (eds), Advances in Neural Information
Processing Systems,  vol. 31.
Curran Associates, Inc.

[14]
Manigrasso, Francesco, Miro, Filomeno Davide, Morra, Lia, & Lamberti,
Fabrizio. 2021.
Faster-LTN: a neuro-symbolic, end-to-end object detection
CoRR, abs/2107.01877.

[15]
Pryor, Connor, Dickens, Charles, Augustine, Eriq, Albalak, Alon, Wang, William,
& Getoor, Lise. 2022.
NeuPSL: Neural Probabilistic Soft Logic.

[16]
Ren, Shaoqing, He, Kaiming, Girshick, Ross B., & Sun, Jian. 2017.
Faster R-CNN: Towards Real-Time Object Detection with Region
Proposal Networks.
IEEE Trans. Pattern Anal. Mach. Intell., 39(6),

[17]
Sarker, Md. Kamruzzaman, Zhou, Lu, Eberhart, Aaron, & Hitzler, Pascal. 2021.
Neuro-Symbolic Artificial Intelligence: Current Trends.
CoRR, abs/2105.05330.

[18]
Serafini, Luciano, & d'Avila Garcez, Artur S. 2016.
Logic Tensor Networks: Deep Learning and Logical Reasoning from Data
and Knowledge.
In: Besold, Tarek R., Lamb, Luís C., Serafini, Luciano,
& Tabor, Whitney (eds), Proceedings of the 11th International Workshop
on Neural-Symbolic Learning and Reasoning (NeSy'16) co-located with the Joint
Multi-Conference on Human-Level Artificial Intelligence (HLAI 2016), New
York City, NY, USA, July 16-17, 2016.
CEUR Workshop Proceedings, vol. 1768.

[19]
Sourek, Gustav, Aschenbrenner, Vojtech, Zelezny, Filip, Schockaert, Steven, &
Kuzelka, Ondrej. 2018.
Lifted Relational Neural Networks: Efficient Learning of
Latent Relational Structures.
Journal of Artificial Intelligence Research, 62(May),

[20]
van Krieken, Emile, Acar, Erman, & van Harmelen, Frank. 2019.
Semi-supervised Learning using Differentiable Reasoning.
FLAP, 6(4), 633–652.

[21]
van Krieken, Emile, Acar, Erman, & van Harmelen, Frank. 2022.
Analyzing Differentiable Fuzzy Logic Operators.
Artificial Intelligence, 302(Jan.), 103602.

[22]
Wagner, Benedikt, & d'Avila Garcez, Artur S. 2022.
Neural-Symbolic Reasoning Under Open-World and Closed-World
In: Martin, Andreas, Hinkelmann, Knut, Fill, Hans-Georg,
Gerber, Aurona, Lenat, Doug, Stolle, Reinhard, & van Harmelen, Frank (eds),
Proceedings of the AAAI 2022 Spring Symposium on Machine Learning and
Knowledge Engineering for Hybrid Intelligence (AAAI-MAKE 2022), Stanford
University, Palo Alto, California, USA, March 21-23, 2022.
CEUR Workshop Proceedings, vol. 3121.

[23]
Winters, Thomas, Marra, Giuseppe, Manhaeve, Robin, & Raedt, Luc De. 2022.
DeepStochLog: Neural Stochastic Logic Programming.
Proceedings of the AAAI Conference on Artificial Intelligence,
36(9), 10090–10100.
Number: 9.