# Pattern Complexity of Aperiodic Substitutive Subshifts

Etienne Moutot<EMAIL_ADDRESS>
Aix-Marseille Université, CNRS, LIS, France  Coline Petit-Jean  coline.petit-
<EMAIL_ADDRESS>
ENS de Lyon, France

###### Abstract

This paper aims to better understand the link better understand the links
between aperiodicity in subshifts and pattern complexity. Our main
contribution deals with substitutive subshifts, an equivalent to substitutive
tilings in the context of symbolic dynamics. For a class of substitutive
subshifts, we prove a quadratic lower bound on their pattern complexity.
Together with an already known upper bound, this shows that this class of
substitutive subshifts has a pattern complexity in $\Theta(n^{2})$. We also
prove that the recent bound of Kari and Moutot, showing that any aperiodic
subshift has pattern complexity at least $mn+1$, is optimal for fixed $m$ and
$n$.

##  1 Introduction

One of the most fascinating aspects of tilings of the plane is the existence
of _aperiodic tilesets_. That is, a set of tiles that tessellates the plane
but only in a non-periodic manner. Aperiodic tilesets even exit in the
specific setting of Wang tilings, where tiles are square with colored edges
and cannot be rotated nor reflected. Two tiles can then be placed next to each
other if the colors of the matching edges are the same. Interestingly,
aperiodic tilesets do not exist for a similar model in dimension one. If one
considers tilings of the infinite line $\mathbb{Z}$ by bi-color Wang dominoes,
then any set of dominoes tiling the line also tiles it periodically. This has
a consequence on the decidability of the most fundamental problem about
tilings: the domino problem. This problem just asks if a given tileset can or
cannot tile the plane. Considering this problem on a line or a plan changes
everything: the domino problem for bi-infinite line is decidable in polynomial
time, as it is equivalent to the existence of a cycle in a finite graph;
whereas it is undecidable if one considers tilings of the plane. From now on,
“tiling” will always mean “Wang tiling”, except stated otherwise.

The first aperiodic tileset is due to Berger in 1966 [3], who disproved Wang’s
conjecture stating that aperiodic tilesets of the plane could not exist. His
tileset was expressed as a set of Wang tiles, a model introduced by Wang to
study fragments of first order logic [26]. Each tile is a unit square with a
color on each edge. Two tiles can then be placed next to each other if the
colors of the matching edges are the same. Then from a finite set of tiles,
one tries to tile the plane by copies of these tiles without rotations. Berger
set of tiles was initially made of 20426 tiles, reduced to 104 in his PhD
thesis [2]. This started the search for the smallest aperiodic set of Wang
tiles. Among many, one can cite Robinson’s tileset [24] which has only 56
tiles. With a radically different idea, Kari built a 14 tiles aperiodic
tileset [14]. His construction was quickly improved by Culik to 13 tiles [6],
which held the record for many years. Finally, in 2015 Jeandel and Rao found
an 11 tiles aperiodic tileset, and proved that it was the smallest possible
[13].

If it is interesting to have aperiodic tilesets with few tiles, the number of
tiles is not a good estimation of how “complex” the resulting tilings can be.
Indeed, with the same number of tiles one can end up with very simple tilings
as well as extremely complex (encoding for example the space-time diagram of a
Turing machine). A better way of estimating how complex a tiling is by its
_pattern complexity_. For a tileset $T$, we denote by $\mathcal{L}_{m,n}(T)$
the set of all $m\times n$ rectangular patterns appearing in at least one
valid tiling by $T$. The pattern complexity of $T$ is then the cardinal of
this set. Besides to carrying more information about how “complicated” a
tiling is, this definition of complexity can be seen as a generalization of
the number of tiles, as the latter is simply the number of patterns of size
$1\times 1$. Computing the exact complexity of a tileset if often complicated,
and a more studied quantity is the topological entropy of a tileset. However,
topological entropy is a very rough estimation of the complexity, and provide
only an asymptotic estimation of the growth of the complexity function. In
this paper, we are interested in bounding the exact complexity of aperiodic
tilesets. More precisely, we want to understand what is the minimal complexity
that can be achieved by an aperiodic tileset; or equivalently, how “simple”
can an aperiodic tileset be. In [15], it is shown that any aperiodic tileset
have complexity at least $mn+1$. The first result of the current paper shows
that this bound is optimal for fixed $m$ and $n$.

One of the most common techniques to build an aperiodic tileset is to use two
dimensional substitutions. In most of the cases, the resulting tilings are
close to being fixpoints of well-chosen two-dimensional substitutions, leading
to their aperiodicity. This is the case of Berger’s tileset, Robinson’s, and
even Jeandel and Rao’s tileset was recently found to have a substitutive
structure [16]. Their substitutive nature also ensures that they have zero
topological entropy. Studying the complexity of substitutive tilesets is
therefore a way of understanding the complexity of many aperiodic tilesets.
One of the rare exceptions is Kari’s aperiodic tileset, which is known not to
be substitutive, as it was shown to have positive topological entropy [7].
This also shows that it has an exponential pattern complexity.

More generally, _substitutive subshifts_ have been extensively studied in
dimension one [8, 10]. A substitutive subshift is a set of colorings of the
infinite line $\mathbb{Z}$ generated by infinite iterations of a one-
dimensional substitution. In dimension one, Pansiot fully classified
substitutions in terms of factor complexity [21].

In the two-dimensional case, substitutions have been mostly studied in the
context of geometrical tilings [1, 9, 12, 22, 25]. In this paper, we focus on
discrete substitutions over two-dimensional words. Our second main result is a
quadratic lower bound on the pattern complexity of a large class of two
dimensional substitutive subshifts. Together with Robinson’s upper bound on
the complexity of two-dimensional substitutions [23, Theorem 7.16], this shows
that the pattern complexity of subshifts from this class is in
$\Theta(n^{2})$. Such a bound paves the way towards a classification of two-
dimensional substitutions in terms of pattern complexity, as Pansiot did for
one-dimensional substitutions.

After introducing definitions and basic properties in Section 2, we show in
Section 3 that the $mn+1$ complexity bound is optimal for fixed $m,n$. Then we
prove our lower bound on our class of substitutive subshifts in Section 4.
Finally we conclude by some remarks on possible directions for future work in
Section 5.

##  2 Preliminaries

In this section, we introduce all the useful definitions and properties.

We will use the following conventions for our notations. Capital $A$
designates a finite alphabet and $a$ and $b$ elements of it. Capital $X$ is a
subshift, whose configurations are usually denoted by the letter $c$, and
patterns by $p$. Vectors of $\mathbb{Z}^{2}$ are bold
$\mathbf{u},\mathbf{v}\in\mathbb{Z}^{2}$. $\sigma$ is a substitution, and
$i,j,k,l,m,n$ designate integers.

###  2.1 Subshifts

Tilings have an equivalent definition as subshifts of finite type. Let $A$ be
a finite alphabet. A _pattern_ of support $D\subset\mathbb{Z}^{2}$ is a
coloring of $D$ with colors of $A$, i.e. $p\in A^{D}$, and a _configuration_
$c\in A^{\mathbb{Z}^{2}}$ is a coloring of $\mathbb{Z}^{2}$. We say that a
pattern $p$ of support $D$ _appears_ in a pattern or a configuration
$p^{\prime}$ if there exists $\mathbf{u}\in\mathbb{Z}^{2}$ such that for all
$\mathbf{v}\in D$, $p_{\mathbf{v}}=p^{\prime}_{\mathbf{v}-\mathbf{u}}$, in
this case we write $p\sqsubseteq p^{\prime}$ (and $p\sqsubset p^{\prime}$ when
$p\neq p^{\prime}$). Let $F$ be a set of patterns, a _subshift_ $X_{F}$ is a
set of configurations in which none of the patterns of $F$ appear:

$X_{F}:=\\{c\in A^{\mathbb{Z}^{2}}\mid\forall p\in F,p\not\sqsubset c\\}.$

If $F$ is finite, $X_{F}$ is called a _subshift of finite type_ (SFT for
short). Subshifts have an equivalent definition as subsets of
$A^{\mathbb{Z}^{2}}$ that are shift-invariant and topologically closed. In
other words, any shift of a configuration of a subshift $X$ is still in $X$,
and the limit of a sequence of shifts of configurations of $X$ is in $X$.

Let $\tau_{\mathbf{u}}$ denote the shift action by the vector
$\mathbf{u}\in\mathbb{Z}^{2}$: let $p$ be a pattern or configuration of
support $D$ (which is $\mathbb{Z}^{2}$ if $p$ is a configuration), then for
all $\mathbf{v}\in
D+\mathbf{u},\tau_{\mathbf{u}}(p)_{\mathbf{v}}=p_{\mathbf{v}-\mathbf{u}}$. A
configuration $c$ is said to be _periodic_ if there exists
$\mathbf{u}\in\mathbb{Z}^{2}\backslash\\{(0,0)\\}$ such that
$c=\tau_{\mathbf{u}}(c)$. A subshift is _aperiodic_ if it is not empty and
contains no periodic configuration.

The set of valid tilings by a set of Wang tiles mentioned in the introduction
is an SFT. These models are actually equivalent, as any SFT can be made into a
set of Wang tilings by application of a single morphism. As subshifts are the
most natural model to deal with substitutions, we chose to use this formalism
in this paper instead of Wang tilesets.

###  2.2 Pattern complexity

The set of all patterns appearing in a configuration $c$ is called the
language of $c$, $\mathcal{L}(c)$; and we denote by
$\mathcal{L}_{m,n}(c):=\mathcal{L}(c)\cap
A^{\\{0,\ldots,m-1\\}\times\\{0,\ldots,n-1\\}}$ the set of all $m\times n$
rectangular patterns appearing in $c$. Similarly,
$\mathcal{L}(X)=\bigcup_{c\in X}\mathcal{L}(c)$ is the set of all patterns of
the subshift $X$ and $\mathcal{L}_{m,n}(X):=\bigcup_{c\in
X}\mathcal{L}_{m,n}(c)$ is the set of $m\times n$ patterns of $X$. The
_pattern complexity_ (or just _complexity_) of a configuration (resp. a
subshift) is the number of such rectangular patterns $|\mathcal{L}_{m,n}(c)|$
(resp. $|\mathcal{L}_{m,n}(X)|$).

For configurations of dimension one, periodicity and complexity are linked:
the Morse-Hedlund theorem states that a one-dimensional configuration $w\in
A^{\mathbb{Z}}$ is periodic if and only if there exists n such that
$P_{c}(n)\leq n$ [18]. For one-dimensional SFTs, the question is trivial as
every non-empty SFT contains a periodic configuration (see for example [17]),
and thus the complexity of the subshift cannot have any impact on its
aperiodicity.

In dimension two, the link between complexity and periodicity becomes much
more involved. Nivat conjectured in 1997 an implication similar to Morse and
Hedlund’s: any two-dimensional configuration with $P_{c}(m,n)\leq m,n$ for
some $m,n$ must be periodic [20]. The other direction does not hold, as there
exists periodic configurations with exponential complexity [4]. Since two
dimensional aperiodic SFTs do exist, it is natural to try to link their
complexity with their aperiodicity, and a bound relating the two has been
recently found:

######  Theorem (Kari, Moutot [15]).

Let $X$ be a subshift such that $\exists c\in X,\exists
m,n,{|\mathcal{L}_{m,n}(c)|\leq mn}$. Then there exists $d\in X$ that is
periodic.

Which can be rephrased as a lower bound on the complexity of aperiodic
subshifts.

######  Corollary 1.

Let $X$ be an aperiodic subshift. Then for all $c\in X$ and all $m,n$,

${|\mathcal{L}_{m,n}(c)|\geq mn+1}.$

###  2.3 Substitutions

Substitutions on words have been extensively studied, see for example [10] for
a good reference about them. In dimension two, they have been studied mostly
in the context of substitution of $\mathbb{R}^{2}$, that is, substitutions
producing geometrical tilings. In the current paper, we restrict ourselves to
discrete uniform two-dimensional substitutions. In this context, a
_substitution_ of size $m\times n$ (or support
$\\{0,\ldots,m-1\\}\times\\{0,\ldots,n-1\\}$) is a map $\sigma:A\rightarrow
A^{\\{0,\ldots,m-1\\}\times\\{0,\ldots,n-1\\}}$. It is said to be _primitive_
if there exists $k$ such that for all $a,b\in A,$ we have
$a\sqsubseteq\sigma^{k}(b)$, and _invertible_ if no two letters have the same
image. The _subshift of a substitution_ , also called a substitutive subshift,
is the set of configurations that can be made from applying the substitution
to letters of $A$:

$X^{\sigma}:=\left\\{c\in A^{\mathbb{Z}^{2}}\mid\forall p\sqsubset c,\exists
a\in A,\exists k\in\mathbb{N},p\sqsubseteq\sigma^{k}(a)\right\\}.$

In general this subshift is only a sofic subshift and not an SFT, as shown by
Mozes [19].

Substitutions are a very convenient way of generating aperiodicity, by example
by considering fixpoints configurations. In the context of geometrical
tilings, Solomyak characterized when the subshift of a substitution is
aperiodic or not [25]. We recall here his result in the specific case of
uniform rectangular substitutions.

######  Definition 1 (Unique desubstitution).

A $\mathcal{D}$-pattern $p\in A^{\mathcal{D}}$ can be _uniquely desubstitued_
by a substitution $\sigma$ of size $(k,l)$ if there exists
$\mathbf{t}\in\mathbb{Z}^{2}$ which is unique modulo $(k,l)$ such that there
exists $c\in X^{\sigma},\tau_{\mathbf{t}}(\sigma(c))|_{\mathcal{D}}=p$.

######  Theorem 2 (Solomyak [25]).

Let $\sigma$ be a primitive substitution. Then $X^{\sigma}$ is aperiodic if
and only if there exists $\rho>0$ such that all patterns $p$ whose support
fits in $\\{1,\ldots,\rho\\}^{2}$ can be uniquely desubstituted by $\sigma$.

A quadratic upper bound is already known for substitutive subshifts, due to
Hansen and Robinson.

######  Theorem (Hansen, Robinson [23, Theorem 7.17]).

Let $\sigma$ be a primitive and invertible substitution with square support.
Then there exists $K>0$ such that

$|\mathcal{L}_{n,n}(X)|\leq Kn^{2}.$

Our lower bound will not deal with all invertible substitutions, but only with
the ones having what we call a _determining position_ , which is a position in
the image that allows us to uniquely recover the antecedent.

######  Definition 2 (Determining position).

A substitution $\sigma$ has a _determining position_ $(i,j)$ if for all
$a,b\in A$:

$\sigma(a)_{i,j}=\sigma(b)_{i,j}\Rightarrow a=b.$

Equivalently, this means that $\left\\{\sigma(a)_{i,j}\mid a\in A\right\\}=A$.
See Fig. 1 for an illustration of this property.

Figure 1: Illustration of Definition 2: $\sigma_{1}$ has a determining
position $(0,0)$ whereas $\sigma_{2}$ has no determining position.

######  Remark 1.

If $|A|=2$, all non-constant substitutions have a determining position.
Indeed, as the images of the two letters must be different, they must differ
in at most one position, which is then a determining position.

##  3 Optimality of the $mn+1$ bound

Corollary 1 shows that all configurations in aperiodic subshifts must have
complexity at least $mn+1$ for all $m,n$.

It is not known whether this bound is optimal in general, however if $m$ and
$n$ are fixed, it is possible to prove its optimality.

######  Theorem 3.

For all $m,n\in\mathbb{N}$, there exists an aperiodic SFT $X$ such that for
all $c\in X,$ we have $|\mathcal{L}_{m,n}(c)|=mn+1$.

###### Proof.

Let $m,n\in\mathbb{N}$ be fixed and $Y\subseteq A^{\mathbb{Z}^{2}}$ be an
aperiodic SFT on some alphabet $A$ (for example the set of all Robinson
tilings [24]). We will “blow up” $Y$ by encoding each of its colors into a
large rectangle over alphabet $\\{0,1\\}$.

Let $k$ be such that $2^{k^{2}}\geq|A|$. Then, for each $a\in A$, it is easy
to build a $k\times k$ square $S(a)\in\\{0,1\\}^{\\{0,\ldots,k-1\\}^{2}}$
encoding $a$ in such a way that $S$ is an isomorphism. $S$ can be naturally
extended to patterns or configurations of $A^{\mathbb{Z}^{2}}$.

Let $R_{0}$ and $R_{1}$ be the $m\times n$ rectangles with respectively a $0$
and a $1$ in position $(0,0)$ and 0s elsewhere. We define the substitution

$\sigma:\begin{cases}0\mapsto R_{0}\\\ 1\mapsto R_{1}\end{cases}.$

Figure 2: Illustration of the sparse encoding of Theorem 3 with $k=3,m=4$ and
$n=2$. We first apply $S$ to some $s\in A$, and then apply $\sigma$ to this
image.

Let $X:=\sigma(S(Y))=\\{\sigma(S(c))\mid c\in Y\\}$. To prove that $X$ is an
SFT, let $F_{0}$ be the finite set of fobidden patterns of $Y$. Let

$F_{1}=\left\\{p\in\\{0,1\\}^{\\{0,\ldots,2km-1\\}\times\\{0,\ldots,2kn-1\\}}\mid\forall
q\in
A^{\\{0,\ldots,3\\}\times\\{0,\ldots,3\\}},p\not\sqsubset\sigma(S(q))\right\\}$

bet the set of $2km\times 2kn$ patterns not appearing in any $2km\times 2kn$
patterns of the form $\sigma(S(q))$. Let $F=F_{1}\cup\sigma(S(F_{0}))$. Then,
$X=X_{F}$. Indeed, $F_{1}$ ensures that all configurations of $X_{F}$ are
desubstituable by $\sigma(S)$, and together with $\sigma(S(F_{0}))$, we have
$c\in X_{F}\Leftrightarrow c\in\sigma(S(Y))$ thanks to the fact that
$\sigma(S)$ is one-to-one. As $F$ is finite, $X$ is an SFT.

Now, since every $c\in X$ can be written $\sigma(d)$ with $d\in S(X_{0})$,
every $m\times n$ rectangular pattern of $c$ contains at most one $1$,
therefore $|\mathcal{L}_{m,n}(c)|\leq mn+1$. But because $X$ is aperiodic,
Corollary 1 ensures that $|\mathcal{L}_{m,n}(c)|\geq mn+1$, so
$|\mathcal{L}_{m,n}(c)|=mn+1$. ∎

One might hope to improve the $mn+1$ bound of Corollary 1 by using a result of
Cassaigne:

######  Theorem (Cassaigne, direct consequence of [5]).

Let $c\in A^{\mathbb{Z}^{2}}$ be such that for all
$m,n,{|\mathcal{L}_{m,n}(c)|=mn+1}$. then there exists
$d\in\overline{\mathcal{O}(c)}$ which is uniform. In particular, the subshift
$\overline{\mathcal{O}(c)}$ is not periodic.

However, it gives an information about configurations with
$|\mathcal{L}_{m,n}(c)|=mn+1$ for _all_ $m,n$. In order to improve Corollary 1
to $mn+2$, one would need a similar result with the hypothesis of the
existence of $m,n$ such that $|\mathcal{L}_{m,n}(c)|=mn+1$. Together with
Theorem 3, this emphasizes the importance of the quantifiers and their order
when studying the complexity of subshifts.

The SFT built in Theorem 3 depends on $m$ and $n$, so it does not rule out the
possibility of improving the uniform bound in $m$ and $n$. In the next section
we show that the uniform lower bound can indeed be improved for a large class
of substitutive subshifts.

##  4 Lower bound for substitutive subshifts

In this section, we show our lower bound for a class of two-dimensional
substitutive subshifts. First, we need a classical lemma on substitutive
subshifts, whose (short) proof is included for sake of completeness.

######  Lemma 4.

Let $\sigma$ be a primitive substitution and $c\in X^{\sigma}$. Then for all
$p\sqsubset c,\sigma(p)\sqsubset c$.

###### Proof.

Let $p\sqsubset c$. By definition of $X^{\sigma}$, there exists $a\in A$,
$k\in\mathbb{N}$ such that $p\sqsubset\sigma^{k}(a)$. This implies that
$\sigma(p)\sqsubset\sigma^{k+1}(a)$. Moreover, $\sigma$ is primitive: there
exists some $l\in\mathbb{N}$ such that for all $b\in A$,
$a\sqsubset\sigma^{l}(b)$ and consequently
$\sigma^{k+1}(a)\sqsubset\sigma^{l+k+1}(b)$.

By definition of $X^{\sigma}$, for a big enough pattern $p^{\prime}$ appearing
in c, there is some $b$ such that $\sigma^{l+k+1}(b)\sqsubset p^{\prime}$.
Eventually, we have:

$\sigma(p)\sqsubset\sigma^{k+1}(a)\sqsubset\sigma^{l+k+1}(b)\sqsubset
p^{\prime}\sqsubset c.$

∎

Then, we take a look at the determining position property for composition of
substitutions.

######  Lemma 5.

If $\sigma$ and $\sigma^{\prime}$ have determining positions, then
$\sigma\circ\sigma^{\prime}$ has a determining position.

###### Proof.

Let $\sigma$ and $\sigma^{\prime}$ be two substitutions of size $(m,n)$ having
determining positions $(k,l)$ and $(k^{\prime},l^{\prime})$ respectively. Let
$a,a^{\prime}\in A$. Then,

$\displaystyle(\sigma\circ\sigma^{\prime})(a)_{(mk+k^{\prime},nl+l^{\prime})}=(\sigma\circ\sigma^{\prime})(a^{\prime})_{(mk+k^{\prime},nl+l^{\prime})}$
$\displaystyle\Rightarrow$
$\displaystyle\sigma\left(\sigma^{\prime}(a)_{(k,l)}\right)_{(k^{\prime},l^{\prime})}=\sigma\left(\sigma^{\prime}(a^{\prime})_{(k,l)}\right)_{(k^{\prime},l^{\prime})}$
$\displaystyle\Rightarrow$
$\displaystyle\sigma^{\prime}(a)_{(k,l)}=\sigma^{\prime}(a^{\prime})_{(k,l)}$
$\displaystyle\Rightarrow$ $\displaystyle a=a^{\prime}.$

∎

An immediate consequence is the following.

######  Corollary 6.

If $\sigma$ has a determining position, then for all $k\geq 1,\sigma^{k}$ has
a determining position.

The next technical lemma is the key of our lower bound. Essentially, it
provides us a lower bound on the complexity based on the number of times $k$
that a patterns can be uniquely desubstituted.

######  Lemma 7.

Let $\sigma$ be a substitution of size $(M,M)$ having a determining position
and $k\in\mathbb{N}$. For all $n\in\mathbb{N}$, if all the patterns of size
$(n,n)$ of $X^{\sigma}$ have a unique desubstitution by $\sigma^{k}$, then for
all $c\in X^{\sigma}$:

$|\mathcal{L}_{n,n}(c)|\geq
n^{2}\left(1+\frac{1}{\lceil\frac{n}{M^{k}}\rceil^{2}}\right).$

###### Proof.

Let us denote by $\mathcal{P}_{c}(k,l)$ the elements of $\mathcal{L}_{k,l}(c)$
for this proof for better reading of the equations. Using the fact that
$\sigma^{k}$ also has a determining position we build an injective function
$f$ from the disjoint union $\bigsqcup_{0\leq i,j\leq
M^{k}-1}\mathcal{P}_{c}(\lceil\frac{n-i}{M^{k}}\rceil,\lceil\frac{n-j}{M^{k}}\rceil)$
to $\mathcal{P}_{c}(n,n)$. The existence of such an injective function, put
together with the result of Kari and Moutot (Corollary 1), yields the result.

Let $(d_{1},d_{2})\in\\{0,\ldots,M^{k}-1\\}$ be the determining position of
$\sigma^{k}$ obtained in Corollary 6, ${i,j\in\\{0,\ldots,M^{k}-1\\}}$ and
${p\in\mathcal{P}_{c}(\lceil\frac{n-i}{M^{k}}\rceil,\lceil\frac{n-j}{M^{k}}\rceil)}$.
The fact that $p$ appears in $c$ implies that it can be extended to bigger
patterns also appearing in $c$. Then, for each such $p$, we fix one
$p^{ext}\in
A^{\\{-1,\dots,\lceil\frac{n-i}{M^{k}}\rceil+1\\}\times\\{-1,\dots,\lceil\frac{n-j}{M^{k}}\rceil+1\\}}$
such that
$p^{ext}|_{\\{0,\dots,\lceil\frac{n-i}{M^{k}}\rceil-1\\}\times\\{0,\dots,\lceil\frac{n-j}{M^{k}}\rceil-1\\}}=p$
and $p^{ext}\sqsubset c$. We define the following pattern:

$f_{i,j}(p)=\tau_{i-d_{1},j-d_{2}}(\sigma^{k}(p^{ext}))|_{\\{0,\dots,n-1\\}^{2}}$

This is well defined because the fact that
$i,j,d_{1},d_{2}\in\\{0,\dots,M^{k}-1\\}$ implies that

$\\{0,\dots,n-1\\}^{2}\subset(i-d_{1},j-d_{2})+\\{-1,\dots,\left\lceil\frac{n-i}{M^{k}}\right\rceil+1\\}\times\\{-1,\dots,\left\lceil\frac{n-j}{M^{k}}\right\rceil+1\\}.$

By Lemma 4 it holds that $\sigma^{k}(p^{ext})\sqsubset c$ and consequently
that $f_{i,j}(p)\in\mathcal{P}_{c}(n,n)$.

(a) $p$

(b) $p^{ext}$

(c) $\sigma^{k}(p^{ext})$ and
$f_{i,j}(p)=\tau_{i-d_{1},j-d_{2}}(\sigma^{k}(p^{ext}))|_{\\{0,\dots,n-1\\}^{2}}$

Figure 3: Illustration of the definition of $f_{i,j}(p)$ with
$n=6,i=j=1,M=2,k=1$ and $(d_{1},d_{2})=(0,0)$. Figure (a) shows
${p\in\mathcal{P}_{c}(\lceil\frac{n-i}{M^{k}}\rceil,\lceil\frac{n-j}{M^{k}}\rceil)}$,
(b) shows the completion to $p^{ext}$ and (c) the image of $p^{ext}$ by
$\sigma^{k}$ and the pattern $f_{i,j}(p)$ extracted from it.

Then, the following function:

$\begin{array}[]{ccccc}f&:&\bigsqcup_{0\leq i,j\leq
M^{k}-1}\mathcal{P}_{c}(\lceil\frac{n-i}{M^{k}}\rceil,\lceil\frac{n-j}{M^{k}}\rceil)&\to&\mathcal{P}_{c}(n,n)\\\
&&p\in\mathcal{P}_{c}(\lceil\frac{n-i}{M^{k}}\rceil,\lceil\frac{n-j}{M^{k}}\rceil)&\mapsto&f_{i,j}(p)\\\
\end{array}$

is well defined. Now we show that it is injective.

Let
$p\in\mathcal{P}_{c}(\lceil\frac{n-i}{M^{k}}\rceil,\lceil\frac{n-j}{M^{k}}\rceil)$,
$p^{\prime}\in\mathcal{P}_{c}(\lceil\frac{n-i^{\prime}}{M^{k}}\rceil,\lceil\frac{n-j^{\prime}}{M^{k}}\rceil)$
and suppose that $f(p)=f(p^{\prime})$. By definition of $f$, there are some
patterns $p^{ext},p^{\prime ext}\sqsubset c$, such that

$p^{ext}|_{\\{0,\dots,\lceil\frac{n-i}{M^{k}}\rceil-1\\}\times\\{0,\dots,\lceil\frac{n-j}{M^{k}}\rceil-1\\}}=p\text{\leavevmode\nobreak\
and \leavevmode\nobreak\ }p^{\prime
ext}|_{\\{0,\dots,\lceil\frac{n-i^{\prime}}{M^{k}}\rceil-1\\}\times\\{0,\dots,\lceil\frac{n-j^{\prime}}{M^{k}}\rceil-1\\}}=p^{\prime}$

and that:

$f(p)=\tau_{i-d_{1},j-d_{2}}(\sigma^{k}(p^{ext}))|_{\\{0,\dots,n-1\\}^{2}}=\tau_{i^{\prime}-d_{1},j^{\prime}-d_{2}}(\sigma^{k}(p^{\prime
ext}))|_{\\{0,\dots,n-1\\}^{2}}=f(p^{\prime})$ (1)

Moreover as $f(p)=f(p^{\prime})$ is a pattern of size $(n,n)$, it has a unique
desubstitution by $\sigma^{k}$ by assumption. Therefore,
$(i-d_{1},j-d_{2})\equiv(i^{\prime}-d_{1},j^{\prime}-d_{2})\mod(M^{k},M^{k})$.
As $i,j,i^{\prime},j^{\prime}\in\\{0,\dots,M^{k}-1\\}$, it implies that
$(i,j)=(i^{\prime},j^{\prime})$. Consequently, Eq. 1 can be rewritten as

$\tau_{i-d_{1},j-d_{2}}(\sigma^{k}(p^{ext}))|_{\\{0,\dots,n-1\\}^{2}}=\tau_{i-d_{1},j-d_{2}}(\sigma^{k}(p^{\prime
ext}))|_{\\{0,\dots,n-1\\}^{2}}.$

That is,

$\sigma^{k}(p^{ext})|_{\\{0,\dots,n-1\\}^{2}-(i-d_{1},j-d_{2})}=\sigma^{k}(p^{\prime
ext})|_{\\{0,\dots,n-1\\}^{2}-(i-d_{1},j-d_{2})}.$

And as $p^{ext}\sqsubseteq p$,

$\sigma^{k}(p)|_{\\{0,\dots,n-1\\}^{2}-(i-d_{1},j-d_{2})}=\sigma^{k}(p^{\prime})|_{\\{0,\dots,n-1\\}^{2}-(i-d_{1},j-d_{2})},$
(2)

on the domain where $\sigma^{k}(p)$ is defined.

Now we show that $p$ and $p^{\prime}$ are the same on all their domain. Let
${\mathbf{u}\in\\{0,\dots,\lceil\frac{n-i}{M^{k}}\rceil-1\\}\times\\{0,\dots,\lceil\frac{n-j}{M^{k}}\rceil-1\\}}$.
First, we have:

$M^{k}\mathbf{u}+(d_{1},d_{2})\in\\{0,\dots,n-1\\}^{2}-(i-d_{1},j-d_{2}).$

Therefore by using the fact that $(d_{1},d_{2})$ is a determining position,
Eq. 2 gives:

$\displaystyle\sigma^{k}(p)_{M^{k}\mathbf{u}+(d_{1},d_{2})}=\sigma^{k}(p^{\prime})_{M^{k}\mathbf{u}+(d_{1},d_{2})},$
$\displaystyle\Rightarrow\leavevmode\nobreak\ $
$\displaystyle\sigma^{k}(p_{\mathbf{u}})_{d_{1},d_{2}}=\sigma^{k}(p^{\prime}_{\mathbf{u}})_{d_{1},d_{2}},$
$\displaystyle\Rightarrow\leavevmode\nobreak\ $ $\displaystyle
p_{\mathbf{u}}=p^{\prime}_{\mathbf{u}}.$

This being true for all
$\mathbf{u}\in\\{0,\dots,\lceil\frac{n-i}{M^{k}}\rceil-1\\}\times\\{0,\dots,\lceil\frac{n-j}{M^{k}}\rceil-1\\}$,
it eventually holds that $p=p^{\prime}$. This concludes the proof of the
injectivity of $f$.

The injectivity of $f$ yields the following inequality:

$\left|\mathcal{L}_{n,n}(c)\right|\geq\left|\bigsqcup_{0\leq i,j\leq
M^{k}-1}\mathcal{P}_{c}(\lceil\frac{n-i}{M^{k}}\rceil,\lceil\frac{n-j}{M^{k}}\rceil)\right|=\sum_{0\leq
i,j\leq
M^{k}-1}\left|\mathcal{P}_{c}(\lceil\frac{n-i}{M^{k}}\rceil,\lceil\frac{n-j}{M^{k}}\rceil)\right|.$

Since $c$ is aperiodic and by Corollary 1, for all $0\leq i,j\leq M^{k}-1$, it
holds that:

$\displaystyle\left|\mathcal{P}_{c}(\lceil\frac{n-i}{M^{k}}\rceil,\lceil\frac{n-j}{M^{k}}\rceil)\right|$
$\displaystyle\geq\lceil\frac{n-i}{M^{k}}\rceil\times\lceil\frac{n-j}{M^{k}}\rceil+1$
$\displaystyle\geq\lceil\frac{n-i}{M^{k}}\rceil\times\lceil\frac{n-j}{M^{k}}\rceil\left(1+\frac{1}{\lceil\frac{n-i}{M^{k}}\rceil\times\lceil\frac{n-j}{M^{k}}\rceil}\right)$
$\displaystyle\geq\lceil\frac{n-i}{M^{k}}\rceil\times\lceil\frac{n-j}{M^{k}}\rceil\left(1+\frac{1}{\lceil\frac{n}{M^{k}}\rceil^{2}}\right)$

Eventually, taking the sum and using a classical property of the ceiling
function [11, p85] gives:

$\displaystyle\left|\mathcal{L}_{n,n}(c)\right|$ $\displaystyle\geq\sum_{0\leq
i,j\leq
M^{k}-1}\left|\mathcal{P}_{c}(\lceil\frac{n-i}{M^{k}}\rceil,\lceil\frac{n-j}{M^{k}}\rceil)\right|$
$\displaystyle\geq\left(1+\frac{1}{\lceil\frac{n}{M^{k}}\rceil^{2}}\right)\sum_{0\leq
i,j\leq
M^{k}-1}\lceil\frac{n-i}{M^{k}}\rceil\times\lceil\frac{n-j}{M^{k}}\rceil$
$\displaystyle=\left(1+\frac{1}{\lceil\frac{n}{M^{k}}\rceil^{2}}\right)n^{2}.$

∎

Our last lemma gives us a bound on how many times patterns from a substitutive
subshift can be desubstitutued thanks to Solomyak result’s (Theorem 2).

######  Lemma 8.

Let $\sigma$ be an invertible substitution. If all patterns
$p\in\mathcal{L}(X^{\sigma})$ whose support contains
$\\{0,\dots,\rho-1\\}^{2}$ can be uniquely desubstitued by $\sigma$; then for
all $k$, all patterns in $\mathcal{L}(X^{\sigma})$ whose support contains
$\\{0,\dots,(\rho+1)M^{k-1}-2\\}^{2}$ can be uniquely desubstited by
$\sigma^{k}$.

###### Proof.

We show the lemma by induction on $k$. The result is true for $k=1$. Suppose
that the result is true for some $k\geq 1$ fixed. Let $p$ be a pattern in
$\mathcal{L}(X^{\sigma})$ with support $\mathcal{D}$ such that
$\\{0,\dots,(\rho+1)M^{k}-2\\}^{2}\subset\mathcal{D}$. Let
$\mathbf{t},\mathbf{t}^{\prime}\in\mathbb{Z}^{2}$ be such that
$p=\tau_{\mathbf{t}}(\sigma^{k+1}(c))|_{\mathcal{D}}=\tau_{\mathbf{t}^{\prime}}(\sigma^{k+1}(c^{\prime}))|_{\mathcal{D}}$
for some $c,c^{\prime}\in X^{\sigma}$. By supposition, and because
$(\rho+1)M^{k}\geq\rho$, the pattern $p$ can be uniquely desubstitued by
$\sigma$. This, and the fact that
$p=\tau_{\mathbf{t}}(\sigma(\sigma^{k}(c)))|_{\mathcal{D}}=\tau_{\mathbf{t}^{\prime}}(\sigma(\sigma^{k}(c^{\prime})))|_{\mathcal{D}}$
implies that $\mathbf{t}\equiv\mathbf{t}^{\prime}\mod(M,M)$. Therefore
$\mathbf{t}^{\prime}=\mathbf{t}+M\mathbf{u}$ for some
$\mathbf{u}\in\mathbb{Z}^{2}$ and then:

$\tau_{\mathbf{t}}(\sigma(\sigma^{k}(c)))|_{\mathcal{D}}=\tau_{\mathbf{t}+M\mathbf{u}}(\sigma(\sigma^{k}(c^{\prime})))|_{\mathcal{D}},$
$\Rightarrow\sigma(\sigma^{k}(c))|_{\mathcal{D}+\mathbf{t}}=\tau_{au}(\sigma(\sigma^{k}(c^{\prime})))|_{\mathcal{D}+\mathbf{t}}$
$\Rightarrow\sigma(\sigma^{k}(c))|_{\mathcal{D}+\mathbf{t}}=\sigma(\tau_{u}(\sigma^{k}(c^{\prime})))|_{\mathcal{D}+\mathbf{t}}.$

Let $\mathcal{R}=\\{0,\dots,(\rho+1)M^{k-1}-2\\}^{2}$. For some
$\mathbf{v}\in\mathbb{Z}^{2}$, we have the following inclusion (see Fig. 4):

$M(\mathcal{R}+\mathbf{v})+\\{0,\dots,M-1\\}^{2}\subset\\{0,\dots,(\rho+1)M^{k}-2\\}^{2}+\mathbf{t}\subset\mathcal{D}+\mathbf{t}.$
(3)

Figure 4: Inclusion of the sets of Eq. 3.

Then, because $\sigma$ is invertible, we have:

$\displaystyle\sigma(\sigma^{k}(c))|_{\mathcal{D}+\mathbf{t}}=\sigma(\tau_{u}(\sigma^{k}(c^{\prime})))|_{\mathcal{D}+\mathbf{t}}$
$\displaystyle\Rightarrow\leavevmode\nobreak\ $
$\displaystyle\sigma(\sigma^{k}(c))|_{M(\mathcal{R}+\mathbf{v})+\\{0,\dots,M-1\\}^{2}}=\sigma(\tau_{u}(\sigma^{k}(c^{\prime})))|_{M(\mathcal{R}+\mathbf{v})+\\{0,\dots,M-1\\}^{2}}$
$\displaystyle\Rightarrow\leavevmode\nobreak\ $
$\displaystyle\sigma^{k}(c)|_{\mathcal{R}+\mathbf{v}}=\tau_{u}(\sigma^{k}(c^{\prime}))|_{\mathcal{R}+\mathbf{v}}$

By induction hypothesis, this pattern can be uniquely desubstitued by
$\sigma^{k}$, therefore it holds that $\mathbf{u}\equiv 0\mod M^{k}$, which
implies that $\mathbf{t}^{\prime}=\mathbf{t}+M\mathbf{u}\equiv\mathbf{t}\mod
M^{k+1}$ and that $p$ can be uniquely desubstitued by $\sigma^{k+1}$. ∎

One can then desubstitute patterns from a substitutive subshifts as many times
as possible, and bringing the two previous lemmas together gives the lower
bound:

######  Theorem 9.

Let $\sigma$ be a primitive square aperiodic substitution having a determining
position. Then there exists $K>1$ such that

$\forall c\in X^{\sigma},\forall
n\in\mathbb{N},\left|\mathcal{L}_{n,n}(c)\right|\geq Kn^{2}$

###### Proof.

By Theorem 2, there exists $\rho$ such that all patterns $p$ appearing in the
configurations of $X^{\sigma}$ in which a square of size $\rho$ fits can be
uniquely desubstitued by $\sigma$. Let $c\in X^{\sigma}$ and $n\in\mathbb{N}$.
Let $k(n)=\max\\{k\in\mathbb{N}\mid\forall p\in\mathcal{L}_{n,n}(c),p\text{
can be uniquely desubstitued by }\sigma^{k}\\}$. By Lemma 7 it holds that:

$|\mathcal{L}_{n,n}(c)|\geq
n^{2}\left(1+\frac{1}{\lceil\frac{n}{M^{k(n)}}\rceil^{2}}\right)$

Additionnaly, $(\rho+1)M^{k(n)}-1>n$ otherwise Lemma 8 would contradict the
maximality of $k(n)$. Therefore, $\rho+1>\frac{n}{M^{k(n)}}$ and
$\rho+1>\lceil\frac{n}{M^{k(n)}}\rceil$. This turns our bound into:

$|\mathcal{L}_{n,n}(c)|\geq n^{2}\left(1+\frac{1}{(\rho+1)^{2}}\right)$

which gives a constant $K:=1+\frac{1}{(\rho+1)^{2}}>1$ which does depend
neither on $c$ nor on $n$. ∎

######  Remark 2.

A similar bound can be found for non-square uniform substitutions. However, if
one wants to know the pattern complexity $|\mathcal{L}_{n,n}(c)|$ for some
$m,n$, the bound depends on how close the shape of the substitution is with
respect to the rectangle $(m,n)$. More precisely, let $\sigma$ be a
substitution of size $(M,N)$. Then for any $R>0$ such that

$\left|\log_{M}(m)-\log_{N}(n)\right|\leq R,$

there exists $K_{R}>1$ such that for all $c\in X^{\sigma}$,

$|\mathcal{L}_{n,n}(c)|\geq K_{R}mn.$

To conclude this section, let us apply our result to Robinson’s tileset. Let
$X^{R}\subseteq A^{\mathbb{Z}^{2}}$ be the SFT constituted of all valid tiling
by Robinson’s tileset $A$ [24].

A natural way of constructing a substitution on $A$ is to “grow” the pattern
of each tile to a $2\times 2$ pattern. Each tile is composed of a back ground
marking (in light-blue and gray on our drawings) and of a main marking (in
black). We generate the rules of our $2\times 2$ substitution $\sigma$ in the
following way:

  1. 1.

The image of each tile has the same background marking,

  2. 2.

the upper right tile of the image of a tile $t$ has the same main marking than
$t$,

  3. 3.

the lower left tile of the image of $t$ has a main marking with four out going
arrows having its blue parts oriented as the background marking of $t$,

  4. 4.

the two other tiles are chosen to have one main arrows and to incoming arrows
on each side so that the main arrow extends the upper right tile.

Fig. 5 illustrates the construction of $\sigma(t)$ for one tile $t\in A$ and
Fig. 6 its iteration. For space considerations, we do not write the full
substitution obtained by this process.

Figure 5: The four steps of construction of the substitution $\sigma$. Figure
6: Portions of patterns generated by successive applications of $\sigma$.

Even though $X^{R}$ is not a substitutive subshift by itself, as substitutive
patterns of its configurations can be separated by “fracture lines”, one can
show that $X^{\sigma}\subset X^{R}$, and therefore for all $n$,
$|\mathcal{L}_{n,n}(X^{\sigma})|\leq|\mathcal{L}_{n,n}(X^{R})|$.

Now, we partition $A$ into two sets $A_{1}$ and $A_{2}$. $A_{1}$ contains all
the tiles with four outgoing arrows, of the form , and $A_{2}$ all the other
tiles. Then, we define

$f:\begin{cases}A&\rightarrow\\{0,1\\}\\\ t&\mapsto\begin{cases}1&\text{ if
}t\in A_{1}\\\ 0&\text{ if }t\in A_{2}\end{cases}\end{cases}.$

$f$ is compatible with $\sigma$ in the sense that if $f(t_{1})=f(t_{2})$ for
two tiles in $A$, then $f(\sigma(t_{1}))=f(\sigma(t_{2}))$. Therefore the
following substitution is well defined:

$\sigma^{\prime}:\begin{cases}\\{0,1\\}&\rightarrow\\{0,1\\}^{\\{0,1\\}^{2}}\\\
a&\mapsto f((\sigma(t))\text{ with }\sigma(t)=a\end{cases}$

Fig. 7 illustrate $\sigma^{\prime}$ the kind of patters it generates.

Figure 7: The substitution $\sigma^{\prime}$ and an example of pattern
generated by it.

Let $X^{\sigma^{\prime}}$ be the subshift of $\sigma^{\prime}$. We have that
$\mathcal{L}(X^{\sigma^{\prime}})=\\{\sigma^{\prime k}(0),\sigma^{\prime
k}(1)\mid k\in\mathbb{N}\\}=\\{f(\sigma^{k}(t))\mid t\in A,k\in\mathbb{N}\\}$.
Therefore for all $n$, it holds that
$\mathcal{L}_{n,n}(X^{\sigma^{\prime}})=f(\mathcal{L}_{n,n}(X^{\sigma}))$.
This implies that
$|\mathcal{L}_{n,n}(X^{\sigma^{\prime}})|\leq|\mathcal{L}_{n,n}(X^{\sigma})|$.
Then, $\sigma^{\prime}$ is primitive with a determining position and
$X^{\sigma^{\prime}}$ is aperiodic since every $1$ corresponds to one corner
of a square in Robinson’s tiling, therefore Theorem 9 gives that there exists
$K>1$ such that $|\mathcal{L}_{n,n}(X^{\sigma^{\prime}})|\geq Kn^{2}$, and
then:

$|\mathcal{L}_{n,n}(X^{R})|\geq Kn^{2}.$

##  5 Perspectives

We believe that Theorem 9 can be extended to any square uniform primitive
substitution, without the assumption of determining position. In the current
proof, this condition allows us to “recover” information on patterns when de-
substituting in Lemma 7, even when the pattern does not contain a full
$\sigma^{k}(a)$ for some $a\in A$. For now, we do not see how to obtain a
quadratic lower bound without this information, but we think that it is
possible to do. Such a generalization would show that all substitutive
subshifts of square primitive substitutions have complexity
$|\mathcal{L}_{n,n}(X^{\sigma})|\in\Theta(n^{2})$. In addition to a better
understanding of aperiodic subshifts, we hope that this is the first step
towards a classification of two-dimensional substitutions in terms of pattern
complexity, similar to Pansiot’s [21] in the one dimensional case.

It is also unknown if the $mn+1$ complexity bound of Corollary 1 can be
improved for non-substititive aperiodic subhifts.

## Acknowledgements

The authors want to thank Nathalie Aubrun for her guidance and lots of helpful
advice, as well as Guilhem Gamard for his careful proofreading.

## References

  * [1] Pierre Arnoux, Valérie Berthé, Thomas Fernique, and Damien Jamet.  Functional stepped surfaces, flips, and generalized substitutions.  Theoretical Computer Science, 380(3):251–265, 2007.  Combinatorics on Words.  URL: https://www.sciencedirect.com/science/article/pii/S0304397507001995, doi:https://doi.org/10.1016/j.tcs.2007.03.031.
  * [2] Robert Berger.  The Undecidability of the Domino Problem.  PhD thesis, Harvard University, 1964.
  * [3] Robert Berger.  The Undecidability of the Domino Problem.  Memoirs of the American Mathematical Society. The American Mathematical Society, 1966.
  * [4] Valérie Berthé and Laurent Vuillon.  Tilings and rotations on the torus: a two-dimensional generalization of sturmian sequences.  Discrete Mathematics, 223(1):27 – 53, 2000.  doi:10.1016/S0012-365X(00)00039-X.
  * [5] Julien Cassaigne.  Double sequences with complexity mn + 1.  Journal of Automata, Languages and Combinatorics, 4:153–170, 01 1999.
  * [6] Karel Culik.  An aperiodic set of 13 wang tiles.  Discrete Mathematics, 160(1):245–251, 1996.  doi:https://doi.org/10.1016/S0012-365X(96)00118-5.
  * [7] Bruno Durand, Guilhem Gamard, and Anaël Grandjean.  Aperiodic tilings and entropy.  In Arseny M. Shur and Mikhail V. Volkov, editors, Developments in Language Theory, page 166–177. Springer International Publishing, 2014.  doi:10.1007/978-3-319-09698-8_15.
  * [8] Fabien Durand, Bernard Host, and Christian Skau.  Substitutional dynamical systems, bratteli diagrams and dimension groups.  Ergodic Theory and Dynamical Systems, 19:953 – 993, 08 1999.  doi:10.1017/S0143385799133947.
  * [9] Thomas Fernique.  Multidimensional sturmian sequences and generalized substitutions.  International Journal of Foundations of Computer Science, 17(03):575–599, 2006.  arXiv:https://doi.org/10.1142/S0129054106004005, doi:10.1142/S0129054106004005.
  * [10] N P Fogg, V Berthé, S Ferenczi, C Mauduit, and A Siegel.  Substitutions in Dynamics, Arithmetics and Combinatorics.  Lecture Notes in Mathematics. Springer Berlin Heidelberg, 2003.  URL: https://books.google.fr/books?id=UltuCQAAQBAJ.
  * [11] Ronald L. Graham, Donald E. Knuth, and Oren Patashnik.  Concrete Mathematics: A Foundation for Computer Science.  Addison-Wesley, 1989.
  * [12] Branko. Grunbaum and G. C. Shephard.  Tilings and patterns, 1987.
  * [13] Emmanuel Jeandel and Michaël Rao.  An aperiodic set of 11 wang tiles.  Advances in Combinatorics, January 2021.  doi:10.19086/aic.18614.
  * [14] Jarkko Kari.  A small aperiodic set of wang tiles.  Discrete Mathematics, 160(1–3):259–264, 1996.
  * [15] Jarkko Kari and Etienne Moutot.  Decidability and periodicity of low complexity tilings.  In Christophe Paul and Markus Bläser, editors, 37th International Symposium on Theoretical Aspects of Computer Science (STACS 2020), volume 154 of Leibniz International Proceedings in Informatics (LIPIcs), pages 14:1–14:12. Schloss Dagstuhl–Leibniz-Zentrum fuer Informatik, 2020.  URL: https://drops.dagstuhl.de/opus/volltexte/2020/11875, doi:10.4230/LIPIcs.STACS.2020.14.
  * [16] Sébastien Labbé.  Substitutive structure of jeandel–rao aperiodic tilings.  Discrete & Computational Geometry, 65(3):800–855, Apr 2021.  doi:10.1007/s00454-019-00153-3.
  * [17] Douglas A Lind and Brian Marcus.  An Introduction to Symbolic Dynamics and Coding.  Cambridge University Press, 1995.
  * [18] Marston Morse and Gustav A. Hedlund.  Symbolic dynamics.  American Journal of Mathematics, 60(4):815–866, 1938.
  * [19] Shahar Mozes.  Tilings, substitutions systems and dynamical systems generated by them.  Journal d’analyse mathématique, 53:139–186, 1988.  doi:10.1007/BF02793412.
  * [20] Maurice Nivat.  Keynote address at the 25th anniversary of eatcs, during icalp 1997, bologna.  1997\.
  * [21] Jean-Jacques Pansiot.  Complexité des facteurs des mots infinis engendrés par morphismes itérés.  In Jan Paredaens, editor, Automata, Languages and Programming, page 380–389. Springer Berlin Heidelberg, 1984.
  * [22] R. Penrose.  Pentaplexity a class of non-periodic tilings of the plane.  The Mathematical Intelligencer, 2(1):32–37, Mar 1979.  doi:10.1007/BF03024384.
  * [23] E. Arthur Robinson.  Symbolic dynamics and tilings of $\mathbb{R}^{d}$, volume 60.  2004\.  doi:10.1090/psapm/060/2078847.
  * [24] Raphael M. Robinson.  Undecidability and nonperiodicity for tilings of the plane.  Inventiones mathematicae, 12(3):177–209, Sep 1971.  doi:10.1007/BF01418780.
  * [25] Boris Solomyak.  Nonperiodicity implies unique composition for self-similar translationally finite tilings.  Discrete & Computational Geometry, 20(2):265–279, 1998.  Publisher: Springer.
  * [26] Hao Wang.  Proving theorems by pattern recognition – ii.  The Bell System Technical Journal, 40(1):1–41, 1961.  doi:10.1002/j.1538-7305.1961.tb03975.x.

##  Appendix A The Robinson tileset

Figure 8: Robinson’s tileset.


# Uncertainty Quantification in Neural Differential Equations

Olga Graf
Technical University of Munich
<EMAIL_ADDRESS>
&Pablo Flores11footnotemark: 1
Pontificia Universidad Católica de Chile
<EMAIL_ADDRESS>
Pavlos Protopapas
Harvard University
<EMAIL_ADDRESS>
&Karim Pichara
Pontificia Universidad Católica de Chile
<EMAIL_ADDRESS>
These authors have contributed equally to this work.

###### Abstract

Uncertainty quantification (UQ) helps to make trustworthy predictions based on
collected observations and uncertain domain knowledge. With increased usage of
deep learning in various applications, the need for efficient UQ methods that
can make deep models more reliable has increased as well. Among applications
that can benefit from effective handling of uncertainty are the deep learning
based differential equation (DE) solvers. We adapt several state-of-the-art UQ
methods to get the predictive uncertainty for DE solutions and show the
results on four different DE types.

##  1 Introduction

Driven by the growing popularity of deep learning, several areas of research
have obtained state-of-the-art performances with deep neural networks (NNs).
Among other applications, deep NNs have been applied for solving differential
equations (DEs) [10] — a fundamental tool for mathematical modeling in
engineering, finance, and the natural sciences. Deep learning based solutions
of DEs have recently appeared in, e.g., [20], [18], [19], [22], [9], [11],
[17], [5], [21], [8]. Typically, the NN itself approximates the solution of a
DE. Thanks to that, parallelization is natural and, in contrast to classical
numerical methods, the solution at any time can be computed without the burden
of having to compute all previous time steps. Furthermore, NNs are continuous
and differentiable.

Until recently, the focus of deep learning was on achieving better accuracy in
the NN predictions, but now it is increasingly being shifted to measuring the
prediction’s uncertainty, especially if the task at hand is safety critical.
Uncertainty quantification (UQ) has been considered for deep models in
computer vision, medical image analysis, bionformatics, etc [1]. Likewise, UQ
is important for deep models that solve DEs. The uncertainty here stems from
the fact that we cannot train a NN on an infinite time and/or space domain.
Therefore, we seek to estimate the solution’s uncertainty in the regions where
the model was not trained. Moreover, another source of uncertainty comes from
the model’s limitations such as its architecture.

To the best of our knowledge, this is the first work to discuss UQ for deep
models that solve DEs. Contrary to common deep learning setup, we solve DEs
without any observed data, relying only on the samples of time and/or space
and on the mathematical statement that relates functions and their
derivatives. This makes the application of existing UQ methods not so
straightforward. In this paper, we make the following contributions:

  1. 1.

We propose an adaptation of the four state-of-the-art UQ methods in deep
learning — _Bayes By Backprop_ [3], _Flipout_ [24], _Neural Linear Model_ [23,
15], and _Deep Evidential Regression_ [2, 13] — to the case of solving DEs.

  2. 2.

We test the above-mentioned methods on four different DE types: linear
ordinary DE (ODE), non-linear ODE, system of non-linear ODEs, and partial DE
(PDE).

##  2 Preliminaries

###  2.1 Solving differential equations with neural networks

A DE can be expressed as $\mathcal{L}u-f=0$, where $\mathcal{L}$ is the
differential operator, $u(\mathbf{x})$ is the solution that we wish to find on
some (possibly multidimensional) domain $\mathbf{x}$, and $f$ is a known
forcing function. We denote the NN approximation of the true solution by
$u_{N}$. To solve the DE, we minimize the square loss of the residual function
$\mathcal{R}(u_{N}):=\mathcal{L}u_{N}-f$, i.e., the optimization objective is

$\min_{\mathbf{w}}(\mathcal{R}^{2}(u_{N})),$ (1)

where $\mathbf{w}$ are the NN parameters. It is also necessary to inform the
NN about any initial and/or boundary conditions,
$u_{\mathrm{c}}=u(\mathbf{x}_{\mathrm{c}})$. One can achieve that in a
straightforward way by adding a penalizing term to the loss function. However,
the exact satisfaction of initial/boundary conditions is not possible in this
case, causing problems in case of high sensitivity to initial conditions, and
also yielding unnecessary local uncertainty at $\mathbf{x}_{\mathrm{c}}$.
Therefore, we employ an alternative approach and consider a transformation of
$u_{N}$ which enforces the initial/boundary conditions and satisfies them by
construction. E.g., in one-dimensional case, given an initial condition
$u_{0}=u(t_{0})$, we consider a transformation
$\tilde{u}_{N}(t)=u_{0}+(1-e^{-(t-t_{0})})u_{N}(t)$. In general, the
transformation has the form
$\tilde{u}_{N}(\mathbf{x})=A(\mathbf{x},\mathbf{x}_{\mathrm{c}},u_{\mathrm{c}})+B(\mathbf{x},\mathbf{x}_{\mathrm{c}})u_{N}(\mathbf{x})$.
Hereinafter, $\tilde{u}_{N}$ will denote the enforced solution. Accordingly,
we replace $u_{N}$ by $\tilde{u}_{N}$ in the optimization objective (1).
Besides the advantage of satisfying the initial/boundary conditions exactly,
the latter approach can also reduce the effort required during training [12].

###  2.2 Uncertainty quantification under Bayesian framework

While classical learning considers _deterministic model parameters_
$\mathbf{\theta}$, the Bayesian framework introduces uncertainty by
considering a _posterior distribution over the model parameters_ ,
$p(\mathbf{\theta}|\mathcal{D})$, obtained after observing some data
$\mathcal{D}$. The posterior distribution is given by Bayes’ theorem,
$p(\mathbf{\theta}|\mathcal{D})=p(\mathcal{D}|\mathbf{\theta})\cdot
p(\mathbf{\theta})\hskip 2.0pt/\hskip 2.0ptp(\mathcal{D})$, where
$p(\mathcal{D}|\mathbf{\theta})$ is the likelihood, $p(\mathbf{\theta})$ is
the prior distribution over the parameters, and $p(\mathcal{D})$ is the
evidence. The predictions $y$ at a new test point $\mathbf{x}$ are given by
the posterior predictive distribution,

$p(y|\mathbf{x},\mathcal{D})=\textstyle{\int}{p(y|\mathbf{x},\mathbf{\theta})\cdot
p(\mathbf{\theta}|\mathcal{D})d\mathbf{\theta}}.$ (2)

For probabilistic deep models, there are two main strategies of estimating
(2).

Inference through the posterior distribution of model parameters. As stated in
(2), the posterior predictive is obtained by averaging over the posterior
uncertainty in the model parameters. Thus, we can start with estimating the
posterior distribution of the NN weights.

In this case, well-suited is Bayesian NN [14] which places a prior
distribution on all the weights (and biases) $\mathbf{w}$. Since an analytical
solution for the posterior is intractable for Bayesian NNs, we have to use
numerical approximation methods such as MCMC or variational methods. Despite
the need for sampling in both cases, variational methods are computationally
less expensive for high-dimensional parameter spaces and also provide an
analytical approximation.

_Bayes By Backprop_ (_BBB_) is a variational, backpropagation-compatible
method for training a Bayesian NN. Its optimization objective seeks to
minimize the Kullback-Leibler divergence between the true posterior and the
variational posterior which is re-parametrized as
$\mathcal{N}(\mu,\sigma=\log(1+\exp(\rho)))$ to allow for backpropagation. At
each optimization step, weights $\mathbf{w}=\mu+\sigma\circ\epsilon$, where
$\epsilon\sim\mathcal{N}(0,I)$ and $\circ$ is pointwise multiplication, are
obtained by sampling from the variational posterior.

_BBB_ is followed by _Flipout_ which adds a pseudo-independent perturbation to
the weights at each training point $\mathbf{x}_{n}$ in the mini-batch, namely,
$\mathbf{w}_{n}=\mu+(\sigma\circ\epsilon)R_{n}$, where $R_{n}$ is the random
sign matrix. Intuitively, the weights get flipped symmetrically around the
mean with probability $0.5$.

_Neural Linear Model_ (_NLM_) is an alternative to Bayesian NN. It places a
prior distribution only on the last layer’s weights, and learns point
estimates for the remaining layers. One can interpret the output of these
layers as a basis defined by the feature embedding of the data. The last layer
of _NLM_ performs Bayesian linear regression on this feature basis. _NLM_
provides tractable inference under the Gaussian assumption on likelihood; we
get analytical solution for the posterior distribution.

Inference through the higher-order evidential distribution. It is also
possible to infer parameters of the posterior predictive directly, using
Bayesian hierarchical modeling [6, 7]. In _Deep Evidential Regression_
(_DER_), the higher-order, evidential prior is placed over the Gaussian
likelihood function. Choosing Normal Inverse-Gamma (NIG) prior yields an
analytical solution for the model evidence which is maximized by the
optimization objective with respect to the NIG hyperparameters. _DER_ also
proposes an evidence regularizer which minimizes evidence on incorrect
predictions. The posterior predictive mean and variance are computed
analytically using the learned hyperparameters.

##  3 Uncertainty quantification in neural differential equations

Instead of learning a deterministic solution $\tilde{u}_{N}$, we now aim to
learn a probabilistic solution $u_{\theta}$, characterized by a posterior
predictive distribution. We estimate it using some probabilistic model
$g_{\theta}$ parametrized by $\theta$.

###  3.1 Proposed approach

UQ methods described in Section 2.2 rely upon the assumption that the
likelihood function is Gaussian, centered at the model’s prediction and
evaluated at the observed data points. Namely, _DER_ and _NLM_ use it to
derive the analytical form of a loss function and a posterior distribution,
respectively. _BBB_ and _Flipout_ can in principle use any likelihood in the
loss function, but it has to be of known analytical form. In case of DEs, a
natural way of computing likelihood is to evaluate it at the _residual_
$\mathcal{R}$ _on the training domain_ $\mathbf{x}_{\mathrm{t}}$, which can be
seen as a counterpart to the observed data points in the classical setting.
However, we are left with an open problem of choosing the underlying
distribution for the likelihood function. It makes sense to assume that the
probability density is high enough for values close to zero, but no further
assumptions immediately follow. E.g., it may happen that the limitations of NN
architecture do not allow for the perfect fit, i.e., the distribution of
residuals is not centered around zero. To circumvent this problem, we propose
an alternative way of computing likelihood.

In this first work on UQ for neural DE solvers, we will focus on comparing
predictions outside of the training domain given by different UQ methods,
leaving the detailed treatment of the model fit and its associated uncertainty
for future work. Although in Bayesian framework this uncertainty also affects
the uncertainty outside of the training domain, we hypothesize that even a
simplified treatment, i.e., without using residuals’ distribution, gives
reasonable uncertainty estimates. We propose a two-stage training procedure:

1) We first train a classical NN on the training domain
$\mathbf{x}_{\mathrm{t}}$ to find a _deterministic solution_ $\tilde{u}_{N}$,

2) We use $\tilde{u}_{N}$ as _observed data_ for our probabilistic model
$g_{\theta}$ and define the likelihood using a Gaussian assumption,
$p(\tilde{u}_{N}|\mathbf{\theta})=\prod_{\mathbf{x}_{\mathrm{t}}}\mathcal{N}(\tilde{u}_{N};u_{\theta},\varepsilon)$.

Now the optimization objective will be trying to align the probabilistic model
with the given reference $\tilde{u}_{N}$ rather than trying to minimize the
residuals at all costs. We note that despite interpreting $\tilde{u}_{N}$ in
stage two as observed data rather than a function that solves DE, its
associated variance $\varepsilon$ is not of aleatoric nature (i.e. irreducible
variance that comes from the noise inherent to the data), as it would be in
the classical regression problem. It can be still interpreted as a source of
epistemic (reducible) uncertainty coming from the NN model limitations. Here,
we consider a simplified treatment of $\varepsilon$. In case of _BBB_ ,
_Flipout_ , and _NLM_ , we pre-define $\varepsilon$ with some small number. In
case of _DER_ , we learn $\varepsilon$ along with posterior predictive
distribution, but the result is not particularly useful since _DER_ does not
have direct access to residuals during learning.

Eventually, the probabilistic model allows us to find the posterior predictive
distribution of $u_{\theta}$. In case of _BBB_ , _Flipout_ , and _NLM_ , we
have $\tilde{g}_{\theta}(\mathbf{x}_{\mathrm{t}},\tilde{u}_{N})=u_{\theta}$,
i.e., the model outputs a single instance $u_{\theta}$. For _BBB_ and
_Flipout_ , the posterior predictive distribution
$p(u_{\theta}|\mathbf{x}_{\mathrm{t}},\tilde{u}_{N})$ is computed as an
approximation of integral (2) using sampling; for _NLM_ , an analytical form
is available. In case of _DER_ , we have
$g_{\theta}(\mathbf{x}_{\mathrm{t}},\tilde{u}_{N})=(\gamma,\nu,\alpha,\beta)$,
where $(\gamma,\nu,\alpha,\beta)$ are the NIG hyperparameters. The mean of the
posterior predictive distribution is equal to $\tilde{\gamma}$ and the
variance is computed using the remaining hyperparameters. We note that the
predictive uncertainty also requires the initial and/or boundary condition
enforcement; this way we are able to eliminate unnecessary uncertainties at
$\mathbf{x}_{\mathrm{c}}$.

Main drawbacks of the current approach are the double computational burden and
the not so useful uncertainty for NN approximation of the true solution in the
traning region; both of them are subject to further improvement.

##  4 Experiments and discussion

We corroborate our theory with experimental results on four equations: 1.
Linear ODE for squared exponential, $\frac{du}{dt}=-2tu$; 2. Non-linear ODE
for Duffing-type oscillator, $\ddot{u}+\omega^{2}u+\epsilon u^{3}=0$; 3.
Lotka-Volterra equations (system of non-linear ODEs), $\dot{u}=\alpha u-\beta
uv\hskip 4.0pt\wedge\hskip 4.0pt\dot{v}=-\delta u+\gamma uv$; 4. Burgers’
equation (non-linear PDE), $\frac{\partial u}{\partial t}+u\frac{\partial
u}{\partial x}=\nu\frac{\partial^{2}u}{\partial x^{2}}$.

Our implementation is based on a DE solver provided by _neurodiffeq_ [4], a
Python package built with PyTorch [16]. Since we are considering relatively
simple DEs, we use networks with one to three fully-connected hidden layers.
For the prior distribution of the weights in _BBB_ , _Flipout_ , and _NLM_ ,
we use flat Gaussian priors with mean zero. In _BBB_ and _Flipout_ , we
estimate the posterior predictive from 1000 samples. In _DER_ , there is no
need for choosing a weight prior, but an appropriate regularization parameter
has to be chosen instead. Here, we tune the regularization parameter manually.

Figure 1: Uncertainty estimation for neural ODEs. Probabilistic models give
low epistemic uncertainty in the training domain and inflate it outside of the
training domain.

We demonstrate the UQ results for ODEs and for PDE in Figure 1 and Figure 2,
accordingly. For all ODEs, the deterministic solution is able to approximate
the true solution well; we incorporate this fact in our Bayesian inference by
choosing small $\varepsilon$ which yields that there is almost no uncertainty
in the training domain. We observe that the epistemic uncertainty away from
the training domain is high enough for all methods, which is our main desired
result in this paper. For Burgers’ equation, however, we see that the NN is
not able to learn the true solution, and our probabilistic model is
underestimating the epistemic uncertainty in the training domain and outside
of it. In this case, either a better deterministic model or a better UQ
methodology is needed. Nevertheless, even for a non-perfect fit, the
uncertainty starts inflating outside of the training domain which proves our
initial hypothesis.

We have witnessed comparable performance in sampling-dependent (_BBB_ ,
_Flipout_) and sampling-free (_NLM_ , _DER_) methods. Given the computational
expense of sampling during Bayesian NN training, the latter two methods could
be preferable in the case of complex DEs on a multidimentional domain.

We believe that further enhancement in terms of diversifying experiments
(e.g., considering more complex high-dimensional DEs) and developing theory
(e.g., calibrating $\varepsilon$ with residuals at each optimization step)
will help the deep learning based DE solutions to outperform classical ones
and lead to their increased presence in applications.

Figure 2: Uncertainty estimation for NN based solution of Burgers’ equation.
We test _Bayes By Backprop_ on an example of a non-linear PDE. Uncertainty
coming from underfitting is not captured well due to our simplified inference
in terms of this type of uncertainty, nevertheless we see that the uncertainty
inflates outside of the training domain.

## References

  * Abdar et al. [2021] M. Abdar, F. Pourpanah, S. Hussain, D. Rezazadegan, L. Liu, M. Ghavamzadeh, P. Fieguth, X. Cao, A. Khosravi, U. R. Acharya, and et al.  A review of uncertainty quantification in deep learning: Techniques, applications and challenges.  _Information Fusion_ , 76:243–297, Dec 2021.  ISSN 1566-2535.  doi: 10.1016/j.inffus.2021.05.008.  URL http://dx.doi.org/10.1016/j.inffus.2021.05.008.
  * Amini et al. [2019] A. Amini, W. Schwarting, A. Soleimany, and D. Rus.  Deep evidential regression.  _arXiv_ , Oct 2019.  URL http://arxiv.org/abs/1910.02600.
  * Blundell et al. [2015] C. Blundell, J. Cornebise, K. Kavukcuoglu, and D. Wierstra.  Weight uncertainty in neural network.  In _International Conference on Machine Learning_ , pages 1613–1622. PMLR, 2015.
  * Chen et al. [2020] F. Chen, D. Sondak, P. Protopapas, M. Mattheakis, S. Liu, D. Agarwal, and M. Di Giovanni.  Neurodiffeq: A python package for solving differential equations with neural networks.  _Journal of Open Source Software_ , 5:1931, 02 2020.  doi: 10.21105/joss.01931.
  * Flamant et al. [2020] C. Flamant, P. Protopapas, and D. Sondak.  Solving differential equations using neural network solution bundles.  _arXiv preprint arXiv:2006.14372_ , 2020.
  * Gelman and et al. [2006] A. Gelman and et al.  Prior distributions for variance parameters in hierarchical models.  _Bayesian analysis_ , 1(3):515–534, 2006.
  * Gelman et al. [2008] A. Gelman, A. Jakulin, M. G. Pittau, and Y.-S. Su.  A weakly informative default prior distribution for logistic and other regression models.  _The Annals of Applied Statistics_ , 2(4):1360–1383, 2008\.
  * Giovanni et al. [2020] M. D. Giovanni, D. Sondak, P. Protopapas, and M. Brambilla.  Finding multiple solutions of ODEs with neural networks.  In _AAAI Spring Symposium: MLPS_ , 2020.
  * Hagge et al. [2017] T. Hagge, P. Stinis, E. Yeung, and A. M. Tartakovsky.  Solving differential equations with unknown constitutive relations as recurrent neural networks, 2017.
  * Lagaris et al. [1998] I. Lagaris, A. Likas, and D. Fotiadis.  Artificial neural networks for solving ordinary and partial differential equations.  _IEEE Transactions on Neural Networks_ , 9(5):987–1000, 1998.  ISSN 1045-9227.  doi: 10.1109/72.712178.  URL http://dx.doi.org/10.1109/72.712178.
  * Mattheakis et al. [2021] M. Mattheakis, H. Joy, and P. Protopapas.  Unsupervised reservoir computing for solving ordinary differential equations, 2021.
  * McFall and Mahan [2009] K. McFall and J. Mahan.  Artificial neural network method for solution of boundary value problems with exact satisfaction of arbitrary boundary conditions.  _IEEE Transactions on Neural Networks_ , 20(8):1221–1233, 2009.  doi: 10.1109/TNN.2009.2020735.
  * Meinert and Lavin [2021] N. Meinert and A. Lavin.  Multivariate deep evidential regression.  _arXiv preprint arXiv:2104.06135_ , 2021.
  * Neal [2012] R. M. Neal.  Bayesian learning for neural networks.  _Springer Science & Business Media_, 118, 2012.
  * Ober and Rasmussen [2019] S. W. Ober and C. E. Rasmussen.  Benchmarking the neural linear model for regression.  In _Advances in Approximate Bayesian Inference (AABI)_ , pages 2171–2180, 2019.
  * Paszke et al. [2019] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, A. Desmaison, A. Köpf, E. Yang, Z. DeVito, M. Raison, A. Tejani, S. Chilamkurthy, B. Steiner, L. Fang, J. Bai, and S. Chintala.  Pytorch: An imperative style, high-performance deep learning library, 2019\.
  * Paticchio et al. [2020] A. Paticchio, T. Scarlatti, M. Mattheakis, P. Protopapas, and M. Brambilla.  Semi-supervised neural networks solve an inverse problem for modeling Covid-19 spread.  _arXiv preprint arXiv:2010.05074_ , 2020.
  * Piscopo et al. [2019] M. L. Piscopo, M. Spannowsky, and P. Waite.  Solving differential equations with neural networks: Applications to the calculation of cosmological phase transitions.  _Physical Review D_ , 100(1), Jul 2019.  ISSN 2470-0029.  doi: 10.1103/physrevd.100.016002.  URL http://dx.doi.org/10.1103/PhysRevD.100.016002.
  * Raissi [2018] M. Raissi.  Forward-backward stochastic neural networks: Deep learning of high-dimensional partial differential equations.  _arXiv preprint arXiv:1804.07010_ , 2018.
  * Raissi et al. [2019] M. Raissi, P. Perdikaris, and G. Karniadakis.  Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations.  _Journal of Computational Physics_ , 378:686–707, 2019\.  ISSN 0021-9991.  doi: https://doi.org/10.1016/j.jcp.2018.10.045.  URL https://www.sciencedirect.com/science/article/pii/S0021999118307125.
  * Randle et al. [2020] D. Randle, P. Protopapas, and D. Sondak.  Unsupervised learning of solutions to differential equations with generative adversarial networks, 2020.
  * Sirignano and Spiliopoulos [2018] J. Sirignano and K. Spiliopoulos.  DGM: A deep learning algorithm for solving partial differential equations.  _Journal of Computational Physics_ , 375:1339–1364, Dec 2018.  ISSN 0021-9991.  doi: 10.1016/j.jcp.2018.08.029.  URL http://dx.doi.org/10.1016/j.jcp.2018.08.029.
  * Snoek et al. [2015] J. Snoek, O. Rippel, K. Swersky, R. Kiros, N. Satish, N. Sundaram, M. M. A. Patwary, Prabhat, and R. P. Adams.  Scalable Bayesian optimization using deep neural networks.  In _32nd International Conference on Machine Learning (ICML)_ , pages 2171–2180, 2015.
  * Wen et al. [2018] Y. Wen, P. Vicol, J. Ba, D. Tran, and R. Grosse.  Flipout: Efficient pseudo-independent weight perturbations on mini-batches.  _arXiv preprint arXiv:1803.04386_ , 2018.