# New Bounds on Diffsequences

Alexander Clifton111Dept. of Mathematics, Emory University, Atlanta, USA.
Email<EMAIL_ADDRESS>Research supported in part by a George W. Woodruff
Fellowship and NSF Grant DMS-1945200.

###### Abstract

For a set of positive integers $D$, a $k$-term $D$-diffsequence is a sequence
of positive integers $a_{1}<a_{2}<\cdots<a_{k}$ such that $a_{i}-a_{i-1}\in D$
for $i=2,3,\cdots,k$. For $k\in\mathbb{Z}^{+}$ and $D\subset\mathbb{Z}^{+}$,
we define $\Delta(D,k)$, if it exists, to be the smallest integer $n$ such
that every $2$-coloring of $\\{1,2,\cdots,n\\}$ contains a monochromatic
$D$-diffsequence of length $k$. We improve the lower bound on $\Delta(D,k)$
where $D=\\{2^{i}\mid i\in\mathbb{Z}_{\geq{0}}\\}$, proving a conjecture of
Chokshi, Clifton, Landman, and Sawin. We also determine all sets of the form
$D=\\{d_{1},d_{2},\cdots\\}$ with $d_{i}\mid d_{i+1}$ for which $\Delta(D,k)$
exists.

##  1 Introduction

Ramsey theory questions typically ask about whether a graph contains a
monochromatic copy of a particular subgraph when the edges are colored with
$r$ colors. However, more generally, the question at the heart of Ramsey
Theory is whether a set retains a certain property when the elements are
divided into disjoint subsets. For example, the positive integers contain
arbitrarily long arithmetic progressions, but if we split them into $r$
disjoint sets, will there still be a set with arbitrarily long arithmetic
progressions? In 1927, van der Waerden [Wae27] answered this question in the
affirmative by showing that every $r$-coloring of the natural numbers will
contain arbitrarily long monochromatic arithmetic progressions.

It is natural to ask what other types of structures besides arithmetic
progressions exhbit van der Waerden-type theorems. Van der Waerden’s work is
predated by Schur [Sch17] who showed that every $r$-coloring of
$\mathbb{Z}^{+}$ contains some monochromatic triple $(x,y,z)$ satisfying
$x+y=z$. More recent sequences of study include arithmetic progressions modulo
$m$ [LW97, LL94], arithmetic progressions with a fixed common difference
[LW97, LL94, BGL99], and descending waves [BEF90].

Landman and Robertson considered the question of whether an $r$-coloring must
produce a monochromatic sequence whose gaps lie in some fixed set. The answer
to this question will depend on the choice of $r$ as well as the choice of the
set of gaps. In [LR03], Landman and Robertson introduced the notion of a $D$\-
diffsequence:

######  Definition 1.

For a set $D$ of positive integers, a $D$-diffsequence of length $k$ is a
sequence of positive integers $a_{1}<a_{2}\cdots<a_{k}$ such that

$a_{i}-a_{i-1}\in D$

for $i=2,3,\cdots,k$.

######  Definition 2.

A set $D$ is called $r$-accessible if every $r$-coloring of the positive
integers contains arbitrarily long monochromatic $D$-diffsequences.

Note that any $r$-accessible set $D$ is automatically $t$-accessible for
$t<r$. The notion of $r$-accessibility is analogous to van der Waerden’s
theorem with $r$ colors.

The set of Fibonacci numbers is known to be $2$-accessible [LR03] but not
$5$-accessible [AGJ+08]. The set of primes is not $3$-accessible [LR03] but
the question of whether it is $2$-accessible remains unresolved. Results on
accessibility for fixed translates of the set of primes are given in [LR03]
and [LV10]. The set of powers of $n$ is only $2$-accessible when $n=2$ since
otherwise there is a periodic coloring modulo $n-1$ which avoids arbitrarily
long $D$-diffsequences.

Just as the van der Waerden number $W(r,k)$ is the smallest $n$ such that any
$r$-coloring of $[n]$ contains a monochromatic arithmetic progression of
length $k$, we can consider analogues in the diffsequence setting. If $D$ is
$r$-accessible, we let $\Delta(D,k;r)$ denote the smallest positive integer
$n$ such that any $r$-coloring of $[n]:=\\{1,2,\cdots,n\\}$ contains a
monochromatic $D$-diffsequence of length $k$. For simplicity, we let
$\Delta(D,k):=\Delta(D,k;2)$.

Landman and Robertson showed that the set $D=\\{2^{i}\mid
i\in\mathbb{Z}_{\geq{0}}\\}$ is $2$-accessible with
$\Delta(D,k)\leq{2^{k}-1}$. Chokshi, Clifton, Landman, and Sawin [CCLS18]
conjectured the existence of an exponential lower bound. We verify this
conjecture, establishing that $\Delta(D,k)$ is at least exponential in
$\sqrt{k}$.

######  Theorem 1.1.

When $D=\\{2^{i}\mid i\in\mathbb{Z}_{\geq{0}}\\}$,

$\Delta(D,k)\geq{2^{\sqrt{2k}}\left(\frac{(\sqrt{2}-1)k}{8}-\frac{\sqrt{k}}{8}\right)+\frac{\sqrt{k}}{2}}.$

We next consider some sets that are not $2$-accessible.

######  Proposition 1.2.

The set $D=\\{k!\mid k\in\mathbb{Z}^{+}\\}$ is not $2$-accessible.

Generalizing both powers of $2$ and factorials, we consider sets of the form

$D_{\\{a_{n}\\}}:=\\{\displaystyle\prod_{i=1}^{k}a_{i}\mid
k\in\mathbb{Z}^{+}\\}$

for a fixed sequence of positive integers $a_{1},a_{2},a_{3},\cdots$ with
$a_{i}\geq{2}$ when $i\geq{2}$. The elements of the set $D_{\\{a_{n}\\}}$,
arranged in increasing order, form a dividing sequence where each term is a
factor of all subsequent terms.

Note that $\\{a_{n}\\}=(1,2,2,2,2,\cdots)$ recovers powers of $2$ while
$\\{a_{n}\\}=(1,2,3,4,5,\cdots)$ recovers factorials.

We are able to completely characterize when $D_{\\{a_{n}\\}}$ is
$2$-accessible:

######  Theorem 1.3.

$D_{\\{a_{n}\\}}$ is $2$-accessible if and only if $\\{a_{n}\\}$ contains
arbitrarily long strings of consecutive $2$’s.

##  2 Powers of $2$

We will prove Theorem 1.1 by utilizing a series of periodic colorings to
obtain the lower bound.

Let $P_{1}$ be a periodic coloring with a repeating block of “$10$”. This
means that odd numbers are assigned color $1$ and even numbers are assigned
color $0$.

Let $P_{2}$ be a periodic coloring with a repeating block of “$1001$”. This
means numbers that are $1$ or $4\pmod{4}$ are color $1$ while numbers that are
$2$ or $3\pmod{4}$ are color $0$.

Continuing in this manner, we define $P_{t}$ as a periodic coloring with a
repeating block of size $2^{t}$. The first half of the repeating block is the
repeating block of $P_{t-1}$ while the second half is the complement of the
first half. Thus, by construction, any two numbers $x$ and $x+2^{t-1}$ have
different colors with respect to $P_{t}$.

The repeating block of $P_{t}$ can be thought of as the first $2^{t}$ bits of
the Thue-Morse sequence. Note that $P_{3}$ is the periodic coloring that
Landman and Robertson used to obtain their initial linear lower bound in
[LR03].

A monochromatic $D$-diffsequence $a_{1}<a_{2}<\cdots<a_{k}$ with respect to
the coloring $P_{t}$ has no gaps of size $2^{t-1}$. Furthermore, we can upper
bound the number of gaps of each size less than $2^{t-1}$ in a way that does
not depend on $k$:

######  Lemma 2.1.

For $m=0,1,\cdots,t-2$, there are at most $m+1$ gaps of size $2^{m}$.

###### Proof.

Consider a monochromatic $D$-diffsequence with respect to $P_{t}$. For each
$m=0,1,\cdots,t-2$, we can split each repeating block of size $2^{t}$ into
sub-blocks of size $2^{m+1}$. If $a_{i+1}-a_{i}=2^{m}$, then by construction,
$a_{i}$ cannot be in the first half of a sub-block of size $2^{m+1}$. So for
any gap $a_{i+1}-a_{i}$ of size $2^{m}$,
$a_{i}\in\\{2^{m}+1,2^{m}+2,\cdots,2^{m+1}\\}\pmod{2^{m+1}}$. We say that
$a_{i}$ is in the second half mod $2^{m+1}$ while $a_{i+1}$ is in the first
half.

If the gaps $a_{i+2}-a_{i+1},a_{i+3}-a_{i+2},\cdots,a_{j}-a_{j-1}$ are all
larger than $2^{m}$, then $a_{j}$ is in the first half mod $2^{m+1}$ and thus
$a_{j+1}-a_{j}$ cannot be a gap of size $2^{m}$. This means that in between
every two gaps of the same size, there is at least one gap of a smaller size.
Since there is no way to have a gap of size smaller than $1$, there is at most
one gap of size $1$.

We will prove the lemma by induction. Suppose that there exists some
$1\leq{q}\leq{t-2}$ such that there are at most $m+1$ gaps of size $2^{m}$ for
$m=0,1,2,\cdots,q-1$. We will now show that there are at most $q+1$ gaps of
size $2^{q}$.

Assume for the sake of contradiction that there are at least $q+2$ gaps of
size $2^{q}$. The first $q+2$ of these are the gaps between $a_{i_{j}}$ and
$a_{i_{j}+1}$ for some indices $i_{j}$ with $j=1,2,\cdots,q+2$. For
$j=1,2,\cdots,q+1$, the sequence of gaps
$a_{i_{j}+2}-a_{i_{j}+1},a_{i_{j}+3}-a_{i_{j}+2},\cdots,a_{i_{j+1}}-a_{i_{j+1}-1}$,
consists of at least one power of $2$ smaller than $2^{q}$, as well as
possibly some powers of $2$ larger than $2^{q}$, which we can ignore. For
$j=1,2,\cdots,q+1$, let $A_{j}$ denote the sum of the gaps
$a_{i_{j}+2}-a_{i_{j}+1},a_{i_{j}+3}-a_{i_{j}+2},\cdots,a_{i_{j+1}}-a_{i_{j+1}-1}$
which are not divisible by $2^{q}$. For each $j=1,2,\cdots,q+1$, we have that
$(a_{i_{1}}+2^{q})+(A_{1}+2^{q})+(A_{2}+2^{q})+\cdots+(A_{j-1}+2^{q})+A_{j}=a_{i_{1}}+(A_{1}+A_{2}+\cdots+A_{j})+j2^{q}$
is in the second half mod $2^{q+1}$. This means that
$a_{i_{1}}+(A_{1}+A_{2}+\cdots+A_{j})$ is in the second half mod $2^{q+1}$ if
$j$ is even (this includes when $j=0$ and we are just considering $a_{i_{1}}$)
and is in the first half mod $2^{q+1}$ if $j$ is odd.

Beginning with $a_{i_{1}}$, every time we add the next $A_{j}$, we switch
halves mod $2^{q+1}$. This means that $A_{1}+A_{2}+\cdots+A_{r}>(r-1)2^{q}$
for any $r=1,2,\cdots,q+1$. In particular,
$A_{1}+A_{2}+\cdots+A_{q+1}>q2^{q}$. So there exist gaps of size
$1,2,4,\cdots,2^{q-1}$ which collectively sum to at least $q2^{q}+1$. By the
inductive hypothesis, we have at most $m+1$ gaps of size $2^{m}$ for
$m=0,1,\cdots,q-1$. Thus the gaps of size $1,2,4,\cdots,2^{q-1}$ collectively
sum to at most:

$\displaystyle\displaystyle\sum_{m=0}^{q-1}(m+1)2^{m}$
$\displaystyle=\sum_{m=0}^{q-1}(2^{q}-2^{m})$
$\displaystyle=q2^{q}-\sum_{m=0}^{q-1}2^{m}$ $\displaystyle=q2^{q}-(2^{q}-1)$
$\displaystyle=(q-1)2^{q}+1.$

This contradicts that these must sum to at least $q2^{q}+1$, so it’s
impossible for a monochromatic $D$-diffsequence with respect to $P_{t}$ to
have $q+2$ gaps of size $2^{q}$. Thus, there are at most $m+1$ gaps of size
$2^{m}$ for $m=0,1,2,\cdots,q$ and the induction is complete. ∎

Lemma 2.1 gives a constant upper bound on the number of gaps smaller than
$2^{t}$ in a monochromatic $D$-diffsequence with respect to $P_{t}$ and thus a
lower bound in terms of $k$ on the number of gaps of size at least $2^{t}$.
This gives a lower bound for $\Delta(D,k)$ that is linear in terms of $k$ and
whose slope depends on $t$. Picking $t=\lfloor{\sqrt{2k}\rfloor}$ for a given
$k$ is in fact enough to obtain an exponential lower bound of

$\Delta(D,k)\geq{2^{\lfloor\sqrt{2k}\rfloor}(\lfloor\sqrt{2k}\rfloor-2)+2}.$

However, we will instead demonstrate a more refined series of colorings which
gives a slightly better lower bound.

For $u\in\mathbb{Z}_{\geq{0}}$, let $P_{t,u}$ be the periodic coloring with a
repeating block of size $2^{t+u}$ obtained by replacing each bit of $P_{t}$ by
$2^{u}$ copies of itself. We call these $2^{u}$ consecutive $1$’s or $2^{u}$
consecutive $0$’s a sub-block. Note that $P_{t,0}$ is simply ${P_{t}}$.

######  Example 2.2.

$P_{3,1}$ has a repeating block of “$1100001100111100$”, containing $8$ sub-
blocks of length $2$.

######  Example 2.3.

$P_{2,2}$ has a repeating block of “$1111000000001111$”, containing $4$ sub-
blocks of length $4$.

If $a_{i}$ and $a_{i+1}$ are consecutive entries in a monochromatic
$D$-diffsequence with respect to the coloring $P_{t,u}$, there are three
possibilities:

  * i)

$a_{i}$ and $a_{i+1}$ are in different positions within the same sub-block, or

  * ii)

$a_{i}$ and $a_{i+1}$ are in different sub-blocks (possibly in the same block)
but the same position within their respective sub-blocks, or

  * iii)

$a_{i}$ and $a_{i+1}$ are in different sub-blocks and different positions
within their respective sub-blocks.

In the second case, $a_{i+1}-a_{i}$ is a multiple of $2^{u}$, the length of a
sub-block, so $a_{i+1}$ appears $2^{l}$ sub-blocks after $a_{i}$ for some $l$.
In the third case, we note that $a_{i+1}-a_{i}$ cannot be divisible by the
sub-block size $2^{u}$. Therefore, it is at most $2^{u-1}$, so $a_{i}$ and
$a_{i+1}$ are in consecutive sub-blocks.

If we consider what sub-blocks $a_{1},a_{2},\cdots,a_{k}$ lie in and ignore
those that are not the first $a_{i}$ in their respective sub-block, this
corresponds to a monochromatic $D$-diffsequence with respect to the coloring
$P_{t}$.

######  Example 2.4.

$5,6,10,11,12,28$ is a monochromatic $D$-diffsequence in color $0$ with
respect to the coloring $P_{2,2}$. If we look at the sub-blocks these are in,
we have $2,2,3,3,3,7$. Ignoring repeated occurrences of the same sub-block
leaves us with $2,3,7$ which is a monochromatic $D$-diffsequence of color $0$
with respect to $P_{2}$.

Applying Lemma 2.1 to this $D$-diffsequence gives that for $m=0,1,\cdots,t-2$,
there are at most $m+1$ occurrences of $a_{i}$ and $a_{i+1}$ in the original
sequence at a distance of $2^{m}$ sub-blocks apart. In particular, there is at
most one $i$ for which $a_{i}$ and $a_{i+1}$ are in consecutive sub-blocks.
For $m=1,\cdots,t-2$, there at most $m+1$ values of $i$ with
$a_{i+1}-a_{i}=2^{m+u}$.

In the first case, $a_{i+1}$ is in a later position within the sub-block than
$a_{i}$ is. There is at most one time (the third case), where $a_{i+1}$ can be
earlier in its sub-block than $a_{i}$ is. Thus, if we keep track of changes in
the position of $a_{i}$ within its sub-block, it can increase within the range
$1$ to $2^{u}$, decrease once, then increase again, potentially up to $2^{u}$.
This means the gaps accounted for by the first and third cases sum to at most
$2(2^{u})-1$. For $m=1,\cdots,t-2$, we have an upper bound on the number of
gaps of size $2^{m+u}$. We can then obtain a lower bound on the number of gaps
of size at least $2^{t+u}$. (Note that gaps of size $2^{t-1+u}$ cannot occur
in a monochromatic $D$-diffsequence with respect to $P_{t,u}$.)

Although it is possible for gaps of size $2^{u}$ to occur in a monochromatic
$D$-diffsequence with respect to $P_{t,u}$, there is never any reason to use
them if the goal is to minimize $a_{k}$. This is because a gap of size $2^{u}$
can only occur if two sub-blocks of all $1$’s or two sub-blocks of all $0$’s
are adjacent. In that case, it is better to use $2^{u}$ consecutive gaps of
size $1$ in place of one gap of size $2^{u}$. This means there are at most
$\displaystyle\sum_{m=1}^{t-2}(m+1)$ gaps of size less than $2^{t+u}$
accounted for by the second case, in addition to at most $2^{u+1}-1$ gaps of
size less than $2^{t+u}$ accounted for by the first and third cases. There are
a total of $k-1$ gaps $a_{i+1}-a_{i}$, so at least

$(k-1)-\displaystyle\sum_{m=1}^{t-2}(m+1)-(2^{u+1}-1)$

have size at least $2^{t+u}$. Replacing a smaller gap by another gap of size
at least $2^{t+u}$ would increase the sum of the gap sizes, so to minimize the
total size of the gaps, we will take as many as possible of each size less
than $2^{t+u}$.

Thus, the sum of the $k-1$ gaps is at least:

$\displaystyle(2^{u+1}-1)+\displaystyle\sum_{m=1}^{t-2}(m+1)2^{m+u}+[k-1-(2^{u+1}-1)-\displaystyle\sum_{m=1}^{t-2}(m+1)]2^{t+u}$
$\displaystyle=$ $\displaystyle
2^{t+u}(k-2^{u+1}-\frac{(t-2)(t+1)}{2})+(2^{u+1}-1)+(2^{u+t-1}-2^{u+1})+\displaystyle\sum_{m=1}^{t-2}(2^{u+t-1}-2^{u+m})$
$\displaystyle=$ $\displaystyle
2^{t+u}(k-2^{u+1}-\frac{t^{2}}{2}+\frac{t}{2}+1)+(t-2)2^{u+t-1}+2^{u+1}-1$
$\displaystyle=$ $\displaystyle
2^{t+u}(k-2^{u+1}-\frac{t^{2}}{2}+t)+2^{u+1}-1.$

This sum is a lower bound for $a_{k}-a_{1}$ so

$\Delta(D,k)\geq{2^{t+u}(k-2^{u+1}-\frac{t^{2}}{2}+t)+2^{u+1}}.$

We can now strategically choose $t$ and $u$ to maximize the lower bound.
Suppose $t=\left\lfloor\sqrt{2k}\right\rfloor$ and
$u=\left\lfloor\log_{2}{\frac{\sqrt{k}}{2}}\right\rfloor$. Then we have:

$\displaystyle\Delta(D,k)$
$\displaystyle\geq{2^{\left\lfloor\sqrt{2k}\right\rfloor+\left\lfloor\log_{2}{\frac{\sqrt{k}}{2}}\right\rfloor}\left(k-2^{\left\lfloor\log_{2}{\frac{\sqrt{k}}{2}}\right\rfloor+1}-\frac{\left\lfloor\sqrt{2k}\right\rfloor^{2}}{2}+\left\lfloor\sqrt{2k}\right\rfloor\right)+2^{\left\lfloor\log_{2}{\frac{\sqrt{k}}{2}}\right\rfloor+1}}$
$\displaystyle\geq{\frac{2^{\sqrt{2k}}(\frac{\sqrt{k}}{2})}{4}\left(k-\sqrt{k}-k+\sqrt{2k}-1\right)+\frac{\sqrt{k}}{2}}=2^{\sqrt{2k}}\left(\frac{(\sqrt{2}-1)k}{8}-\frac{\sqrt{k}}{8}\right)+\frac{\sqrt{k}}{2}.$

##  3 Factorials

We would like to show that the set, $D$, of factorials is not $2$-accessible.
This means we must demonstrate a coloring of $\mathbb{Z}^{+}$ which avoids
arbitrarily long $D$-diffsequences. We use Beatty sequences to construct the
following $2$-coloring:

The integer $n$ is assigned color $0$ or $1$ depending on the parity of
$\lfloor{n\alpha\rfloor}$ where

$\alpha:=2-\frac{e}{2}-\frac{1}{2e}=1-\displaystyle\sum_{i=1}^{\infty}\frac{1}{(2i)!}.$

Here we prove Proposition 1.2 by showing that this coloring avoids
$D$-diffsequences of length $4$. That is, there are no
$a_{1}<a_{2}<a_{3}<a_{4}$ in $\mathbb{Z}^{+}$ such that
$\lfloor{a_{i}\alpha\rfloor}$ has the same parity for $i=1,2,3,4$ and that
$a_{i+1}-a_{i}\in\\{k!\mid k\in\mathbb{Z}^{+}\\}$ for $i=1,2,3$.

######  Lemma 3.1.

For all $k\geq{1}$, there exists some integer $n$ such that
$2n+\frac{1}{3}\leq{k!\alpha}<2n+1$.

###### Proof.

We split into the cases where $k$ is even and when $k$ is odd:

When $k=2m$ for $m\geq{1}$, we have

$\displaystyle k!\alpha$
$\displaystyle=(2m)!-\displaystyle\sum_{i=1}^{\infty}\frac{(2m)!}{(2i)!}=(2m)!-\displaystyle\sum_{i=1}^{m}\frac{(2m)!}{(2i)!}-\displaystyle\sum_{i=m+1}^{\infty}\frac{(2m)!}{(2i)!}$
$\displaystyle=(2m)!-\displaystyle\sum_{i=1}^{m-1}[(2m)(2m-1)\cdots(2i+1)]-1-\displaystyle\sum_{i=m+1}^{\infty}\frac{1}{(2m+1)(2m+2)\cdots(2i)}$
$\displaystyle=1-\displaystyle\sum_{i=m+1}^{\infty}\frac{1}{(2m+1)(2m+2)\cdots(2i)}\pmod{2}.$

To reach the desired conclusion, it suffices to demonstrate that

$\frac{2}{3}\geq{\displaystyle\sum_{i=m+1}^{\infty}\frac{1}{(2m+1)(2m+2)\cdots(2i)}}>0.$

It is clear that the sum is positive. It can also be bounded above the
following convergent geometric series:

$\displaystyle\displaystyle\sum_{i=m+1}^{\infty}\frac{1}{(2m+1)^{2i-2m}}$
$\displaystyle=\frac{\frac{1}{(2m+1)^{2}}}{1-\frac{1}{(2m+1)^{2}}}=\frac{1}{(2m+1)^{2}-1}$
$\displaystyle\leq{\frac{1}{3^{2}-1}}=\frac{1}{8}.$

When $k=2m+1$ for $m\geq{1}$, we have

$\displaystyle k!\alpha$
$\displaystyle=(2m+1)!-\displaystyle\sum_{i=1}^{\infty}\frac{(2m+1)!}{(2i)!}=(2m+1)!-\displaystyle\sum_{i=1}^{m}\frac{(2m+1)!}{(2i)!}-\displaystyle\sum_{i=m+1}^{\infty}\frac{(2m+1)!}{(2i)!}$
$\displaystyle=(2m+1)!-\displaystyle\sum_{i=1}^{m-1}\frac{(2m+1)!}{(2i)!}-(2m+1)-\displaystyle\sum_{i=m+1}^{\infty}\frac{1}{(2m+2)(2m+3)\cdots(2i)}$
$\displaystyle=1-\displaystyle\sum_{i=m+1}^{\infty}\frac{1}{(2m+2)(2m+3)\cdots(2i)}\pmod{2}.$

To reach the desired conclusion, it suffices to demonstrate that

$\frac{2}{3}\geq{\displaystyle\sum_{i=m+1}^{\infty}\frac{1}{(2m+2)(2m+3)\cdots(2i)}}>0.$

It is clear that the sum is positive. It can also be bounded above by the
following convergent geometric series:

$\displaystyle\displaystyle\sum_{i=m+1}^{\infty}\frac{1}{(2m+2)^{2i-2m-1}}$
$\displaystyle=\frac{\frac{1}{2m+2}}{1-\frac{1}{(2m+2)^{2}}}=\frac{2m+2}{(2m+2)^{2}-1}$
$\displaystyle\leq{\frac{4}{4^{2}-1}}=\frac{4}{15}.$

Only $k=1$ remains to be checked. We have established that
$1-\frac{1}{8}\leq{2\alpha}<1\pmod{2}$. This means that either
$\frac{7}{16}\leq{\alpha}<\frac{1}{2}\pmod{2}$ or
$1+\frac{7}{16}\leq{\alpha}<1+\frac{1}{2}\pmod{2}$. If the latter statement is
true, then $\alpha>1$ or $\alpha<0$ but this is false because
$\alpha<2-\frac{e}{2}<2-\frac{2}{2}=1$ and
$\alpha>2-\frac{3}{2}-\frac{1}{2(2)}=\frac{1}{4}$. Therefore,
$\frac{7}{16}\leq{\alpha}<\frac{1}{2}\pmod{2}$, giving us
$\frac{1}{3}\leq{\alpha}<1\pmod{2}$ as desired. ∎

###### Proof.

of Proposition 1.2
Now suppose that there exist $a_{1},a_{2},a_{3},a_{4}$ all the same color such
that the consecutive gaps lie in $D=\\{i!\mid i\in\mathbb{Z}^{+}\\}$. Since
this sequence is monochromatic,
$\lfloor{a_{i+1}\alpha}\rfloor-\lfloor{a_{i}\alpha}\rfloor$ is even for
$i=1,2,3$. By Lemma 3.1, we know that for $i=1,2,3$,

$2n_{i}+\frac{1}{3}\leq{(a_{i+1}-a_{i})\alpha}<2n_{i}+1$

for some integers $n_{i}$.

If $a_{i}\alpha\in[\frac{2}{3},1)\pmod{2}$, then
$a_{i+1}\alpha\in[1,2)\pmod{2}$. Similarly, if
$a_{i}\alpha\in[\frac{5}{3},2)\pmod{2}$, then $a_{i+1}\alpha\in[0,1)\pmod{2}$.
Thus, in order for $a_{i+1}$ and $a_{i}$ to be the same color, we need
$0\leq{a_{i}\alpha}<\frac{2}{3}\pmod{1}$. In particular, we need
$0\leq{a_{2}\alpha}<\frac{2}{3}\pmod{1}$ and
$0\leq{a_{3}\alpha}<\frac{2}{3}\pmod{1}$.

By similar reasoning, these inequalities force
$0\leq{a_{1}\alpha}<\frac{1}{3}\pmod{1}$ and
$0\leq{a_{2}\alpha}<\frac{1}{3}\pmod{1}$, respectively. However, since
$\lfloor{a_{2}\alpha}\rfloor$ and $\lfloor{a_{1}\alpha}\rfloor$ have the same
parity, this gives $\frac{-1}{3}<a_{2}\alpha-a_{1}\alpha<\frac{1}{3}\pmod{2}$,
contradicting the restriction that
$2n_{1}+\frac{1}{3}\leq{(a_{2}-a_{1})\alpha}<2n_{1}+1$. Therefore, no such
$a_{1},a_{2},a_{3},a_{4}$ exist with respect to this coloring. ∎

##  4 Dividing Sequences

In this section, we prove Theorem 1.3. First we show that $D_{\\{a_{n}\\}}$ is
$2$-accessible whenever $\\{a_{n}\\}$ contains arbitrarily long strings of
consecutive $2$’s. This follows from a similar argument to the one Landman and
Robertson used to show that $\\{2^{i}\mid i\in\mathbb{Z}^{+}\\}$ is
$2$-accessible [LR03].

######  Claim 4.1.

If $\\{a_{n}\\}$ contains arbitrarily long strings of consecutive $2$’s, then
$D_{\\{a_{n}\\}}$ is $2$-accessible.

###### Proof.

Suppose for the sake of contradiction that some coloring of $\mathbb{Z}^{+}$
avoids arbitrarily long monochromatic $D_{\\{a_{n}\\}}$-diffsequences. Suppose
that the longest such diffsequence is of length $k$ and given by
$b_{1}<b_{2}<\cdots<b_{k}$.

Because $\\{a_{n}\\}$ contains arbitrarily long strings of consecutive $2$’s,
there exists some index $j$ such that $a_{j},a_{j+1},\cdots,a_{j+k-1}$ all
equal $2$. Letting $C=\displaystyle\prod_{i=1}^{j-1}a_{i}$, we have that
$2^{t}C\in D_{\\{a_{n}\\}}$ for $t=0,1,\cdots,k$. To avoid extending our
monochromatic diffsequence $b_{1}<b_{2}<\cdots<b_{k}$ to a longer one,
$b_{k}+2^{t}C$ must be the other color for $t=0,1,\cdots,k$. However, this
gives a monochromatic sequence in the other color. The consecutive differences
are

$(b_{k}+2^{t}C)-(b_{k}+2^{t-1}C)=2^{t-1}C$

for $t=1,2,\cdots,k$. As stated above, these are all elements of
$D_{\\{a_{n}\\}}$ meaning we have a monochromatic
$D_{\\{a_{n}\\}}$-diffsequence of length $k+1$, giving a contradiction.

∎

We are left to consider $D_{\\{a_{n}\\}}$ in the case where $\\{a_{n}\\}$ does
not have arbitrarily long strings of consecutive $2$’s. All elements of
$D_{\\{a_{n}\\}}$ are multiples of $a_{1}$, so any monochromatic
$D_{\\{a_{n}\\}}$-diffsequence stays within a residue class modulo $a_{1}$.
Considering each residue class separately, we can avoid the presence of
arbitrarily long monochromatic $D_{\\{a_{n}\\}}$-diffsequences in a residue
class as long as it possible to avoid arbitrarily long monochromatic
$T_{\\{a_{n}\\}}$-diffsequences in $\mathbb{Z}^{+}$ where
$T_{\\{a_{n}\\}}=\\{\frac{d}{a_{1}}\mid d\in D_{\\{a_{n}\\}}\\}$. Thus, it
suffices to only consider when $a_{1}=1$.

If there is no string of $k$ consecutive $2$’s in $\\{a_{n}\\}$, we will show
that there exists a coloring of $\mathbb{Z}^{+}$ which avoids monochromatic
$D_{\\{a_{n}\\}}$-diffsequences of length $2^{k}+1$. As in Section 3, we will
use Beatty sequences to make such a coloring. Every $n\in\mathbb{Z}^{+}$ will
be colored according to the parity of $\lfloor{n\alpha\rfloor}$ where our
choice of $\alpha$ depends on $\\{a_{n}\\}$. For convenience, let
$d_{t}:=\displaystyle\prod_{i=1}^{t}a_{i}$ so that
$D_{\\{a_{n}\\}}=\\{d_{t}\mid t\in\mathbb{Z}^{+}\\}$.

######  Lemma 4.2.

There exists some $\alpha>0$ such that for all $t\in\mathbb{Z}^{+}$, there
exists some integer $n_{t}$ such that

$2n_{t}+\frac{1}{2^{k}}\leq{d_{t}\alpha}\leq{2n_{t}+1}.$

Using this lemma, we can finish the proof of Theorem 1.3.

###### Proof.

Suppose for the sake of contradiction that we have a monochromatic
$D_{\\{a_{n}\\}}$-diffsequence $b_{1}<b_{2}<\cdots<b_{2^{k}+1}$. Since this
sequence is monochromatic,
$\lfloor{b_{i+1}}\alpha\rfloor-\lfloor{b_{i}\alpha\rfloor}$ is even for
$i=1,2,\cdots,2^{k}$. By Lemma 4.2, we know that for $i=1,2,\cdots,2^{k}$,

$2m_{i}+\frac{1}{2^{k}}\leq{(b_{i+1}-b_{i})\alpha}\leq{2m_{i}+1}$ (1)

for some integers $m_{i}$. In order for $b_{i+1}$ and $b_{i}$ to be the same
color, we need $0\leq{b_{i}\alpha}<1-\frac{1}{2^{k}}\pmod{1}$, so
$0\leq{b_{2^{k}}\alpha}<1-\frac{1}{2^{k}}\pmod{1}$. Since
$\lfloor{b_{2^{k}}\alpha}\rfloor$ and $\lfloor{b_{2^{k}-1}\alpha}\rfloor$ have
the same parity, this in turn necessitates

$0\leq{(b_{2^{k}-1})\alpha}<1-\frac{2}{2^{k}}\pmod{1},$

which necessitates

$0\leq{(b_{2^{k}-2})\alpha}<1-\frac{3}{2^{k}}\pmod{1}.$

Continuing in this manner, we get that

$0\leq{b_{2}\alpha}<1-\frac{2^{k}-1}{2^{k}}=\frac{1}{2^{k}}\pmod{1}.$

However, $\lfloor{b_{1}\alpha}\rfloor$ and $\lfloor{b_{2}\alpha}\rfloor$ have
the same parity so either $0\leq{b_{2}\alpha-
b_{1}\alpha}<\frac{1}{2^{k}}\pmod{2}$ or $1<b_{2}\alpha-
b_{1}\alpha<2\pmod{2}$. This directly contradicts Equation 1. ∎

###### Proof.

of Lemma 4.2

For each $t\in\mathbb{Z}^{+}$, we will construct closed intervals
$I_{b}^{t}:=[C_{b}^{t},D_{b}^{t}]$ for $b=1,\cdots,t$ such that if
$d_{b}\gamma\in I_{b}^{t}\pmod{2}$, then $d_{h}\gamma\in I_{h}^{t}\pmod{2}$
for $h=b+1,b+2,\cdots,t$. These will be constructed in such a way that
$\frac{1}{2^{k}}\leq{C_{b}^{t}}<D_{b}^{t}\leq{1}$ for $b=1,\cdots,t$. Then,
for any $\gamma\in I_{1}^{t}$, we have $d_{1}\gamma\in I_{1}^{t}\pmod{2}$ and
thus $d_{h}\gamma\in I_{h}^{t}\pmod{2}\subset[\frac{1}{2^{k}},1]\pmod{2}$ for
$h=1,2,\cdots,t$.

We then define
$J_{t}:=[C_{t},D_{t}]:=\displaystyle\bigcap_{b=1}^{t}I_{1}^{b}$. Thus, for
$\gamma\in J_{t}$, we have $d_{b}\gamma\in[\frac{1}{2^{k}},1]\pmod{2}$ for
$b=1,2,\cdots,t$. Also we have $J_{t+1}\subset J_{t}$ for all
$t\in\mathbb{Z}$, so $C_{t}\leq{C_{t+1}}<{D_{t+1}}\leq{D_{t}}$ for
$t\in\mathbb{Z}$.

Consider the sequence $\\{C_{t}\\}$ for $t=1,2,\cdots$. This is a non-
decreasing sequence and it is bounded above by $1$. Therefore, it converges to
some limit $\alpha$. Similarly, the sequence $\\{D_{t}\\}$ for $t=1,2,\cdots$
is non-increasing and bounded below by $0$ so it also converges to some limit
$\beta\geq{\alpha}$. For each $t\in\mathbb{Z}^{+}$, we have
$C_{t}\leq{\alpha}\leq{\beta}\leq{D_{t}}$. Thus, this $\alpha$ lies in $J_{t}$
for all $t\in\mathbb{Z}^{+}$, meaning that
$d_{t}\alpha\in[\frac{1}{2^{k}},1]\pmod{2}$ for all $t\in\mathbb{Z}^{+}$.

It now suffices to construct each $I_{b}^{t}$. We begin with
$I_{t}^{t}:=[\frac{1}{2},1]$. Then for $b\geq{2}$, we recursively construct
$I_{b-1}^{t}$ from $I_{b}^{t}$ as follows:

  * •

If $a_{b}$ is odd:

$C_{b-1}^{t}=\frac{a_{b}-1+C_{b}^{t}}{a_{b}}\hskip
72.26999ptD_{b-1}^{t}=\frac{a_{b}-1+D_{b}^{t}}{a_{b}}.$

  * •

If $a_{b}$ is even:

$C_{b-1}^{t}=\frac{a_{b}-2+C_{b}^{t}}{a_{b}}\hskip
72.26999ptD_{b-1}^{t}=\frac{a_{b}-2+D_{b}^{t}}{a_{b}}.$

We claim this construction satisfies the requisite properties:

######  Claim 4.3.

For $b=1,2,\cdots,t$,

$\frac{1}{2^{k}}\leq{C_{b}^{t}}<D_{b}^{t}\leq{1}.$

######  Claim 4.4.

If $d_{b}\gamma\in I_{b}^{t}\pmod{2}$, then $d_{h}\gamma\in I_{h}^{t}\pmod{2}$
for $h=b+1,b+2,\cdots,t$.

###### Proof.

of Claim 4.3 We set $C_{t}^{t}:=\frac{1}{2}$ and $D_{t}^{t}:=1$. We have
$D_{b}^{t}>C_{b}^{t}$ for $t=b$ and whenever $D_{b}^{t}>C_{b}^{t}$, we get
that

$a_{b}-2+D_{b}^{t}>a_{b}-2+C_{b}^{t},\hskip
72.26999pta_{b}-1+D_{b}^{t}>a_{b}-1+C_{b}^{t},$

so by induction, $D_{b}^{t}>C_{b}^{t}$ for $b=1,\cdots,t$.

Now we need to show that $C_{b}^{t}\geq{\frac{1}{2^{k}}}$ and
$D_{b}^{t}\leq{1}$ for $b=1,\cdots,t-1$. We will proceed by induction by
showing that for $b=2,\cdots,t$, $C_{h}^{t}\geq{\frac{1}{2^{k}}}$ for
$h=b,b+1,\cdots,t$ implies $C_{b-1}^{t}\geq{\frac{1}{2^{k}}}$ and that
$D_{b}^{t}\leq{1}$ implies $D_{b-1}^{t}\leq{1}$. There are two cases to
consider:

  * i)

$a_{b}>2$ or

  * ii)

$a_{b}=a_{b+1}=a_{b+2}=\cdots=a_{b+q-1}=2$ with $b+q-1\leq{t}$ but either
$a_{b+q}>2$ or $b+q=t+1$. Note that that the sequence $\\{a_{n}\\}$ has no
string of $k$ consecutive $2$’s so $q<k$.

In either case,

$D_{b-1}^{t}\leq{\frac{a_{b}-1+D_{b}^{t}}{a_{b}}}\leq{\frac{a_{b}-1+1}{a_{b}}}=1.$

Since $C_{b}^{t}<D_{b}^{t}\leq{1}$, we know that $1\geq{1-C_{b}^{t}}>0$. In
the first case,

$C_{b-1}^{t}=1-\frac{1-C_{b}^{t}}{a_{b}}.$

for odd $a_{b}$ is minimized when $a_{b}=3$, in which case,

$C_{b-1}^{t}=1-\frac{1-C_{b}^{t}}{3}\geq{\frac{2}{3}}.$

For the first case with even $a_{b}$,

$C_{b-1}^{t}=1-\frac{2-C_{b}^{t}}{a_{b}},$

which is minimized when $a_{b}=4$, yielding,

$C_{b-1}^{t}=1-\frac{2-C_{b}^{t}}{4}\geq{\frac{1}{2}}.$

This means that when $a_{b}>2$, we not only get
$C_{b-1}^{t}\geq{\frac{1}{2^{k}}}$, but we get the potentially stronger claim
of $C_{b-1}^{t}\geq{\frac{1}{2}}$.

Now consider the second case. If $b+q\leq{t}$ and $a_{b+q}>2$, then a
consequence of the first case is that $C_{b+q-1}^{t}\geq{\frac{1}{2}}$. If
$b+q=t+1$, our initial definition of $C_{t}^{t}$ also gives
$C_{b+q-1}^{t}\geq{\frac{1}{2}}$.

Now since $a_{b}=a_{b+1}=\cdots=a_{b+q-1}=2$, we have for
$h=b-1,b,\cdots,b+q-2$,

$\displaystyle C_{h}^{t}$
$\displaystyle=\frac{a_{h+1}-2+C_{h+1}^{t}}{a_{h+1}}$
$\displaystyle=\frac{2-2+C_{h+1}^{t}}{2}=\frac{C_{h+1}^{t}}{2}.$

This means

$\displaystyle
C_{b-1}^{t}=\frac{C_{b+q-1}^{t}}{2^{q}}\geq{\frac{1/2}{2^{q}}}=\frac{1}{2^{q+1}}.$

Recall that $q<k$, so this quantity is at least $\frac{1}{2^{k}}$ and the
proof is complete.

∎

###### Proof.

of Claim 4.4

It suffices to show that if $d_{b}\gamma\in I_{b}^{t}\pmod{2}$, then
$d_{b+1}\gamma\in I_{b+1}^{t}\pmod{2}$.

If $d_{b}\gamma\in I_{b}^{t}\pmod{2}$ and $a_{b+1}$ is odd,

$\displaystyle C_{b}^{t}$
$\displaystyle\leq{d_{b}\gamma}\leq{D_{b}^{t}}\pmod{2}$
$\displaystyle\frac{a_{b+1}-1+C_{b+1}^{t}}{a_{b+1}}$
$\displaystyle\leq{d_{b}\gamma}\leq\frac{a_{b+1}-1+D_{b+1}^{t}}{a_{b+1}}\pmod{2}$
$\displaystyle a_{b+1}-1+C_{b+1}^{t}$
$\displaystyle\leq{a_{b+1}d_{b}\gamma}\leq{a_{b+1}-1+D_{b+1}^{t}}\pmod{2}$
$\displaystyle C_{b+1}^{t}$
$\displaystyle\leq{d_{b+1}\gamma}\leq{D_{b+1}^{t}}\pmod{2}.$

Going from the second line to the third line uses that
$0\leq{C_{b+1}^{t}}\leq{D_{b+1}^{t}}\leq{1}$ and that $a_{b+1}-1$ is even.

Similarly if $d_{b}\gamma\in I_{b}^{t}\pmod{2}$ and $a_{b+1}$ is even,

$\displaystyle C_{b}^{t}$
$\displaystyle\leq{d_{b}\gamma}\leq{D_{b}^{t}}\pmod{2}$
$\displaystyle\frac{a_{b+1}-2+C_{b+1}^{t}}{a_{b+1}}$
$\displaystyle\leq{d_{b}\gamma}\leq\frac{a_{b+1}-2+D_{b+1}^{t}}{a_{b+1}}\pmod{2}$
$\displaystyle a_{b+1}-2+C_{b+1}^{t}$
$\displaystyle\leq{a_{b+1}d_{b}\gamma}\leq{a_{b+1}-2+D_{b+1}^{t}}\pmod{2}$
$\displaystyle C_{b+1}^{t}$
$\displaystyle\leq{d_{b+1}\gamma}\leq{D_{b+1}^{t}}\pmod{2}.$

Going from the second line to the third line uses that
$0\leq{C_{b+1}^{t}}\leq{D_{b+1}^{t}}\leq{1}$ and that $a_{b+1}-2$ is even.

∎

∎

##  5 Concluding Remarks

In this paper, we improved the lower bound for $\Delta(D,k)$ when
$D=\\{2^{i}\mid i\in\mathbb{Z}_{\geq{0}}\\}$, establishing that $\Delta(D,k)$
grows faster than any polynomial. The exponent of the lower bound is
$\Theta(\sqrt{k})$ while the exponent of the upper bound is $\Theta(k)$.
Values computed up to $k=12$ align closely with the lower bound so we
conjecture that it is asymptotically tight:

######  Conjecture 5.1.

When $D=\\{2^{i}\mid i\in\mathbb{Z}_{\geq{0}}\\}$,

$\Delta(D,k)=2^{\Theta(\sqrt{k})}.$

It does not make sense to consider $D=\\{n^{i}\mid
i\in\mathbb{Z}_{\geq{0}}\\}$ for any other integer $n$ because $n^{i}\equiv
1\pmod{n-1}$, allowing us to use a periodic coloring modulo $n-1$ to show $D$
is not $2$-accessible. However, it may be interesting to consider other sets
that exhibit exponential growth such as
$D_{\alpha,f}:=\\{\lfloor{\alpha^{i}}\rfloor\mid i\in\mathbb{Z}_{\geq{0}}\\}$
or $D_{\alpha,c}=\\{\lceil{\alpha^{i}}\rceil\mid i\in\mathbb{Z}_{\geq{0}}\\}$
for $\alpha>1$.

When $\alpha=2^{\frac{1}{n}}$ for some $n\in{\mathbb{Z}_{\geq{0}}}$, then both
$D_{\alpha,f}$ and $D_{\alpha,c}$ contain all powers of $2$ and therefore are
$2$-accessible. While not quite in this form, the set of Fibonacci numbers,
which is $2$-accessible[AGJ+08], exhibits similar growth to $D_{\alpha,f}$ and
$D_{\alpha,c}$ for $\alpha=\frac{1+\sqrt{5}}{2}$.

######  Question 5.2.

Is there any $\alpha\in(1,2)$ for which $D_{\alpha,f}$ or $D_{\alpha,c}$ is
not $2$-accessible?

######  Question 5.3.

Is there any $\alpha>2$ for which $D_{\alpha,f}$ or $D_{\alpha,c}$ is
$2$-accessible?

If the answer to either of these questions is yes, we can also consider
whether the set of such $\alpha$ has nonzero measure.

For sets $D$ which are $r$-accessible, there are examples where
$\Delta(D,k;r)$ is polynomial in $k$ such as $D=\\{a\in\mathbb{Z}^{+}\mid
m\nmid a\\}$ for any integer $m\geq{3}$ and $r=2$ [CCLS18],
$D=\\{2\\}\cup\\{2a-1\mid a\in\mathbb{Z}^{+}\\}$ for $r\leq{3}$ [Lan97], and
$D=\\{ma\mid a\in\mathbb{Z}^{+}\\}$ for any $m\in\mathbb{Z}^{+}$ and any
$r\in\mathbb{Z}^{+}$ . There are also cases like $D=\\{2^{i}\mid
i\in\mathbb{Z}_{\geq{0}}\\}$ with $r=2$ where $\Delta(D,k;r)$ is at most
exponential in $k$, but still grows faster than any polynomial. We consider
what other types of behavior are possible for the function $\Delta(D,k;r)$ for
sets $D$ which are $r$-accessible.

######  Question 5.4.

For which $t\in\mathbb{R}^{+}$ does there exist an $r$-accessible set
$D\subset\mathbb{Z}^{+}$ for which

$\Delta(D,k;r)=e^{\Omega(k^{t})}\hskip 7.22743pt?$

One possible mode of attack is as follows. For any infinite set
$T\in\mathbb{Z}^{+}$, the set $D=\\{i-j\mid i,j\in T,i>j\\}$ is $r$-accessible
for any $r$ [LR03]. This means we can construct arbitrarily sparse sets which
are still $r$-accessible and we generally expect a sparser set to have higher
values for $\Delta(D,k;r)$.

Acknowledgments. The author would like to thank Bruce Landman for introducing
him to diffsequences at the 2015 University of West Georgia REU.

## References

  * [AGJ+08] Hayri Ardal, David S. Gunderson, Veselin Jungić, Bruce M. Landman, and Kevin Williamson.  Ramsey results involving the fibonacci numbers.  Fibonacci Quarterly, 46-47:10–17, 2008.
  * [BEF90] T. C. Brown, P. Erdős, and A. R. Freedman.  Quasi-progressionsand descending waves.  Journal of Combinatorial Theory, Series A, 53:81–95, 1990.
  * [BGL99] Tom C. Brown, Ronald L. Graham, and Bruce M. Landman.  On the set of common differences in van der waerden’s theorem on arithmetic progressions.  Canadian Mathematical Bulletin, 42(1):25–36, 1999.
  * [CCLS18] Kajal Chokshi, Alexander Clifton, Bruce Landman, and Oliver Sawin.  Ramsey functions for sequences with restricted gaps.  Journal of Combinatorial Mathematics and Combinatorial Computing, 105:35–44, 2018.
  * [Lan97] Bruce M. Landman.  Avoiding arithmetic progressions (mod $m$) and arithmetic progressions.  Utilitas Mathematica, 52:173–182, 1997.
  * [LL94] B. M. Landman and A. F. Long.  Ramsey functions for sequences with adjacent differences in a specified congruence class.  Proceedings of the Twenty-fifth Southeastern International Conference on Combinatorics, Graph Theory and Computing (Boca Raton, FL, 1994), 103:3–20, 1994.
  * [LR03] Bruce Landman and Aaron Robertson.  Avoiding monochromatic sequences with special gaps.  SIAM Journal on Discrete Mathematics, 21(3):794–801, 2003.
  * [LV10] Bruce M. Landman and Kevin Ventullo.  Avoiding monochromatic sequences with gaps in a fixed translation of the primes.  Utilitas Mathematica, 82:207–214, 2010.
  * [LW97] Bruce M. Landman and Beata Wysocka.  Collections of sequences having the ramsey property only for few colours.  Bulletin of the Australian Mathematical Society, 55:19–28, 1997\.
  * [Sch17] Issai Schur.  Über die kongruenz $x^{m}+y^{m}\equiv z^{m}\pmod{p}$.  Jahresbericht der Deutschen Mathematiker-Vereinigung, 25:114–117, 1917.
  * [Wae27] B. L. van der Waerden.  Beweis einer baudetschen vermutung.  Nieuw Archief voor Wiskunde, 15:212–216, 1927.


# Model-Free Forecasting of Partially Observable Spatiotemporally Chaotic
Systems

Vikrant Guptaa, Larry K. B. Lib, Shiyi Chena, Minping Wana
aGuangdong Provincial Key Laboratory of Turbulence Research and Applications,
Department of Mechanics and Aerospace Engineering, Southern University of
Science and Technology, Shenzhen 518055, PR China
bDepartment of Mechanical and Aerospace Engineering, Hong Kong University of
Science and Technology, Hong Kong, PR China

###### Abstract

Reservoir computing is a powerful tool for forecasting turbulence because its
simple architecture has the computational efficiency to handle large systems.
Its implementation, however, often requires full state-vector measurements and
knowledge of the system nonlinearities. We use nonlinear projector functions
to expand the system measurements to a high dimensional space and then feed
them to a reservoir to obtain forecasts. We demonstrate the application of
such reservoir computing networks on spatiotemporally chaotic systems, which
model several features of turbulence. We show that using radial basis
functions as nonlinear projectors enables complex system nonlinearities to be
captured robustly even with only partial observations and without knowing the
governing equations. Finally, we show that when measurements are sparse or
incomplete and noisy, such that even the governing equations become
inaccurate, our networks can still produce reasonably accurate forecasts, thus
paving the way towards model-free forecasting of practical turbulent systems.

###### keywords:

Time-series forecasting , Machine learning , Chaos , Turbulence , Reservoir
computing , Echo-state networks

††journal: ArXiv

##  1 Introduction

The availability of large datasets and advances in computational hardware have
led to revolutionary applications of machine learning (ML), such as in speech
recognition [1] and medical diagnosis [2]. This has led to renewed interest in
the application of ML to predict the behaviour of turbulent flow systems, such
as for weather and climate forecasting [3, 4, 5, 6, 7, 8, 9]. Turbulent
systems are defined by processes spanning a wide range of spatial and temporal
scales. These processes are coupled nonlinearly such that even small changes
at the microscopic scale can lead to large changes at the macroscopic scale
[10], a phenomenon poetically termed the “butterfly effect” by Lorenz [11].
Modelling such nonlinearly coupled dynamical systems with recurrent neural
networks (RNNs) would require the network dimension to be exceedingly large,
limiting their applicability to either coarse predictions [3, 4] or
necessitating dimensionality reduction of the input states [12, 13, 14].

Reservoir-computing-based RNNs (RC-RNNs) have been shown to be more
computationally efficient at learning large chaotic systems as compared to
backpropagation RNNs [15, 16]. Their implementation, however, often requires
knowledge of the type of nonlinearities in the governing equations [17, 15].
This hinders the application of RC-RNNs to natural and engineering systems
where full state-vector measurements are often not available and/or the
governing equations are not known. In this paper, we introduce a modification
to RC-RNNs whereby nonlinear projectors are used to first expand the
observables to a high dimensional space before they are fed to the reservoir.
We demonstrate the application of such networks for forecasting
spatiotemporally chaotic systems, which are nonlinear dynamical systems with a
large number of degrees of freedom that can model several important features
of turbulence [18, 19]. In particular, we explore the use of radial basis
functions (RBFs) as nonlinear projectors to enable the networks to capture
complex system nonlinearities without knowing the governing equations. We then
explore scenarios in which only sparse or noisy and incomplete measurements
are available, and compare the predictions from RC-RNNs with those from
conventional time-marching.

##  2 Network details

We choose the echo-state network (ESN) of Jaeger [20] for its ability to
forecast chaotic systems. However, we emphasize that our work is applicable to
other RNNs as well. In ESNs the input state
$\mathbf{u}(t)\in\mathbb{R}^{d_{u}}$, which consists of $d_{u}$ concurrent
measurements in space, is connected to a high-dimensional dynamical system,
say $\mathcal{R}$, known as the reservoir. The reservoir has the activation
state $\mathbf{r}(t)\in\mathbb{R}^{d_{r}}$, which is effectively a nonlinear
function of the input history $\mathbf{u}(t),\mathbf{u}(t-\Delta t),\dots$,
where $\Delta t$ is a discrete time step (Eq. 1a). The output
$\mathbf{v}(t)\in\mathbb{R}^{d_{v}}$ is then obtained as a linear readout from
$\mathbf{r}(t)$ using Eq. 1b.

$\mathbf{r}(t)=\mathcal{N}\left(A\mathbf{r}(t-\Delta
t)+W_{in}\mathbf{u}(t)\right),$ (1a) $\mathbf{v}(t)=W_{out}\mathbf{r}(t).$
(1b)

Here, $W_{in}\in\mathbb{R}^{d_{r}\times d_{u}}$ connects the input state
$\mathbf{u}(t)$ to the reservoir nodes, $A\in\mathbb{R}^{d_{r}\times d_{r}}$
is an adjacency matrix that interconnects the reservoir nodes, $\mathcal{N}$
is a known nonlinear function ($\tanh$ here) that acts on each element
separately (i.e.
$\mathcal{N}(\left[a_{1},a_{2},\dots\right])=\left[\mathcal{N}(a_{1}),\mathcal{N}(a_{2}),\dots\right]$)
and $W_{out}\in\mathbb{R}^{d_{v}\times d_{r}}$ connects the reservoir nodes to
the output (see the schematic in Fig. 1 (a)).

Figure 1: (a) Schematic of a standard ESN. The input nodes are connected to a
high-dimensional reservoir, which is the only hidden layer in the network. A
large enough reservoir can capture the system nonlinearities from which the
output (future predictions here) is obtained via linear readout. (b) The
network in the prediction phase is closed-loop, i.e. the output of the network
becomes its input at the next step. (c) Schematic of the ESN from [17]; it has
two improvements: (i) it is parallelized and (ii) the reservoir state is
expanded to a higher-dimension (such as from $r$ to $[r,r^{2}]$). (d)
Schematic of the proposed modified ESN; the observed vector $\mathbf{u}(t)$ is
expanded to $\mathbf{u}_{nl}(t)$ via nonlinear projectors $\mathscr{A}_{1}$,
$\mathscr{A}_{2}$, etc (i.e.
$u_{nl}(t)=[u(t),\mathscr{A}_{1}(u(t)),\mathscr{A}_{2}(u(t)),...]$. This means
that the nonlinearity in the modified network is introduced before the
reservoir instead of after it as in the network in (c).

The elements of $W_{in}$ and $A$ are fixed by drawing them from random
distributions, such as from a uniform distribution in $[-1,1]$. The elements
of $A$ are scaled to keep its spectral radius ($\rho$) usually below one so
that the network satisfies the fading memory property, i.e. the future state
depends only on a finite sequence of past inputs. Finally, the elements of
$W_{out}$ are optimized via a one-step linear regression that minimizes the
mean-square error between $\mathbf{v}(t)$ and the target signal
$\widehat{\mathbf{v}}(t)$ during the training time $-T\leq t\leq 0$. In
forecasting problems, which are of interest here, the target signal
$\widehat{\mathbf{v}}(t)=g(\mathbf{u}(t+1))$, where $g$ is a bijective
function. Consequently, the network becomes closed-loop in the prediction
phase, i.e. the model output becomes its input at the next step as shown in
Fig. 1 (b). See A for further details and [21] for practical guidance on
designing ESNs.

If the complex nonlinearities of turbulent systems are to be captured, $d_{r}$
may need to be several orders of magnitude larger than $d_{u}$. This can be
computationally costly for turbulent systems for which $d_{u}$ is already
large. Recognizing this, Pathak et al. (2018) [17] proposed two improvements
(see Fig. 1 (c)). First, they parallelized the ESNs based on the spatial
locality of turbulent interactions, which facilitates parallel computing. The
parallelized network has two additional hyperparameters: the number of
parallel reservoirs $N_{g}$ and the interaction length $N_{i}$. Second, they
explicitly expanded the reservoir activation state to a higher-dimensional
space such that the output is obtained as
$\mathbf{v}(t)=W_{out}\left[\mathbf{r},\mathbf{r}^{2}\right]$. The main
drawback of this approach is that in order to capture the system
nonlinearities, the nonlinearity added explicitly should closely match that
present in the system equations themselves [15].

We will show that the network from Pathak et al. (2018) [17] can be improved
if the expansion to the high-dimensional space happens before the reservoir.
In other words, we expand the observables (i.e. $u\in\mathbb{R}^{d_{u}}$) to a
high-dimensional space using nonlinear projector functions $\mathscr{A}_{l}$
($l=1,2,...,L$). The expanded vector $u_{nl}\in\mathbb{R}^{(L+1)d_{u}}$ is
then the input to the reservoir from which the output is trained via linear
readout (see Fig. 1 (d)).

##  3 Predictions from existing and modified networks

Figure 2: (a) The exact signal for the KS system and its predictions from
different RC-RNNs: KS$r$ (an instance of a standard network shown in Fig. 1
(a)), KS$r^{2}$ (an instance of the network from [17] shown in Fig. 1 (c)) and
KS$y^{2}$ (an instance of the proposed network shown in Fig. 1 (d)). All
networks here are of size $d_{r}=1000$ and $N_{g}=1$ (i.e. the reservoir size
is only an order of magnitude larger than $d_{u}$). The KS$r$ network, which
has no explicitly added nonlinearity, cannot capture the system nonlinearities
(unless maybe for much larger $d_{r}$). The networks with added square
nonlinearities (KS$r^{2}$ and KS$y^{2}$) are able to capture the system
nonlinearities well. The mean prediction errors (b) and the valid prediction
time (c) from the different networks are also shown. The modified network
KS$y^{2}$ has relatively slower growth of error as compared to the network of
[17] (KS$r^{2}$). The improvement is significant considering that the
KS$r^{2}$ network would require four times more computational power to achieve
similar performance to the KS$y^{2}$ network. The KS$G_{rbf}$ network is also
an instance of the proposed network but with input to reservoir as
$[y,\exp(-0.03y^{2})]$, this network is discussed with Fig. 4 (d-f).

We consider the Kuramoto–Sivashinsky (KS) equation as our first example:

$y_{t}=-yy_{x}-y_{xx}-\nu y_{xxxx},\;\;\;\;\;\;x\in[0,L].$ (2)

The scalar field $y(x,t)$ is periodic in space. The subscripts $x$ and $t$
denote the partial derivatives in space and time, respectively. The parameter
$\nu$ acts like the viscosity in the Navier–Stokes equations [22] and $L$ is
the system size with which the attractor dimension increases linearly [23]. We
numerically solve Eq. 2 for $(\nu,L)=(1,22)$ with $N=128$ gridpoints and
timestep $dt=0.25$ [24]. In this study $\Delta t$ of networks is set as $dt$.

Fig. 2 (a) shows the exact signal from the simulations and its prediction from
different ESNs. All the networks feature a single reservoir of size
$d_{r}=1000$, which is only an order of magnitude larger than $d_{u}=128$. The
networks are trained for $-10000\leq t\leq 0$, training with much smaller time
history is possible but the predictive accuracy will then be lower as shown in
[17]. We note that the KS$r$ network, which is an instance of the standard
network (Fig. 1 (a)), is unable to learn the system dynamics, with the
corresponding forecast soon becoming uninterpretable. Nevertheless, the
results drastically improve when the networks with added nonlinearities are
used. Here KS$r^{2}$ and KS$y^{2}$ are instances of the network from [17]
(Fig. 1 (c)) and the proposed network (Fig. 1 (d)), respectively. Both
networks predict the exact signal in the short-term after which the patterns
in the predictions appear to be spatiotemporally shifted. This shift is an
inherent property of chaotic systems: their trajectories diverge exponentially
in time even when they start from infinitesimally close initial conditions.
The rate of separation is given by $\exp(\Lambda_{max}t)$, where
$\Lambda_{max}$ is the maximal Lyapunov exponent and is positive for chaotic
systems [25].

We quantify the prediction error as
$\mathcal{E}(t)=\sqrt{\langle(\mathbf{v}(t)-\widehat{\mathbf{v}}(t))^{2}\rangle}/\sigma_{y}$,
where $\langle.\rangle$ denotes spatial averaging and $\sigma_{y}$ is the
standard deviation of $y(x,t)$. To account for variability when comparing the
performance of different ESNs, we take the mean of the error
($\overline{\mathcal{E}(t)}$) for a series of predictions (see B for details).
Fig. refKSpreds (b) shows $\overline{\mathcal{E}(t)}$ for predictions from
various ESNs. The mean error increases exponentially from any infinitesimal
differences in the predictions at $t=0$. A slower increase in
$\overline{\mathcal{E}(t)}$ would reduce the effect of both measurement and
model uncertainties on future predictions. However, we note that the mean-
square-error contains only limited information [26]. For the predictions to be
physically interpretable, their probability distribution should be close to
that of the exact signal [27]. (This is the reason for the popularity of
generative-adversarial networks [28], which can give better interpretable
results often at the cost of increased mean-square-error [29].) If the
predictions are physically interpretable in the long-term,
$\overline{\mathcal{E}(t)}$ would saturate to a value close to $\sqrt{2}$.

As expected from the predictions in Fig. 2 (a), $\overline{\mathcal{E}(t)}$
from the standard network KS$r$ not only increases rapidly but it also does
not saturate to $\sqrt{2}$, i.e. the predictions are not interpretable. For
all other networks used here, $\overline{\mathcal{E}(t)}$ saturates to
$\sqrt{2}$ thus indicating statistically consistent predictions. We note that
the modified network KS$y^{2}$ has smaller growth-rate of error as compared to
that for the network from [17] (KS$r^{2}$). In fact, the KS$r^{2}$ network
requires almost four times more computational power (i.e. $N_{g}=4$) to
achieve a similar level of accuracy as the KS$y^{2}$ network. Although both
networks have explicitly added square nonlinearities, the KS$y^{2}$ network
outperforms the KS$r^{2}$ network. This is because in the KS$y^{2}$ network
the reservoir receives more information about the system. In other words, the
reservoir in the KS$y^{2}$ network is preconditioned with square nonlinearity
and can thus better approximate the system nonlinearity.

Another metric recently introduced for evaluating the predictive performance
of models is the valid prediction time (VPT) [16]. It refers to the time at
which $\overline{\mathcal{E}(t)}$ crosses a threshold; the threshold is set as
0.5 following [16]. The advantage of using VPT is that it is a scalar measure
that quantifies the predictive accuracy of different methods. Networks with
slower growth-rate of the mean-square error also have larger VPT as shown in
Fig. 2 (c) for the KS system; this trend is independent of the threshold
value.

##  4 Partially observable systems: selection of nonlinear projectors

Figure 3: (a) The observed signal, $y_{1}$, of the GS system and its
predictions from the GS$P^{3}$, GS$P^{5}$, GS$G_{rbf}$ and GS$2G_{rbf}$
networks of size $d_{r}=1500$ and $N_{g}=4$ (i.e. the total reservoir size is
6000, which is only an order of magnitude larger than $d_{u}$). Although all
four networks can predict the short-term dynamics, the networks with monomial
projectors (GS$P^{3}$ and GS$P^{5}$) fail to give statically consistent
results while the networks with Gaussian RBF projectors (GS$G_{rbf}$ and
GS$2G_{rbf}$) give statistically consistent results. Also shown are the mean
prediction errors (b) and the valid prediction time (c) from the different
ESNs. In addition to statistically consistent results, the networks with
Gaussian RBF projectors also show slower growth-rates of the error and hence
higher VPTs.

For the KS system, we showed that the networks with explicitly added square
nonlinearities (KS$r^{2}$ and KS$y^{2}$) work well. As noted in [15], this is
because the added nonlinearity in the networks matches the nonlinearity in the
system equations (i.e. square nonlinearity in the KS system here). The form of
system nonlinearity, however, is not always known; for example the governing
equations are not known or full state-vector measurements are not available.
To illustrate, we consider the Gray–Scott (GS) system, which models the
reaction-diffusion of two species, as our second example,

$y_{1t}=\left(F-y_{1}y_{2}^{2}\right)-Fy_{1}+D_{1}y_{1xx},$ (3a)
$y_{2t}=y_{1}y_{2}^{2}-(F+k)y_{2}+D_{2}y_{2xx},$ (3b)

where $y_{1}(x,t)$ and $y_{2}(x,t)$ are the concentration of two chemical
species. We choose $D_{1}=2.0e^{-5}$, $D_{2}=1.0e^{-5}$, $F=0.015$, $k=0.0415$
and $x\in[0,1]$ with periodic boundary conditions for which the system
exhibits annihilation type spatiotemporal chaos [30]. We numerically solve the
GS system with 200 chebyshev points and $dt=0.5$ and save the results on an
equispaced grid containing 256 nodes [31]. The time evolution of the
concentration of species one ($y_{1}$) is shown in Fig. 3 (a).

We now assume that only $y_{1}$ measurements are available
($\mathbf{u}(t)=y_{1}(x,t)$), i.e. we do not know the full state-vector. One
may argue that we can still form the full state-vector in terms of
$(y_{1},y_{1t})$ or calculate $y_{2}$ using $(y_{1},y_{1t},y_{1xx})$. This
will, however, require knowledge of Eq. 3 and noise free measurements of
$y_{1}$ with high spatiotemporal resolution. Neither are usually available in
practical situations, nor are they needed by our method (see the results in
Figs. 5 and 6).

With only $y_{1}$ measurements available, the form of system nonlinearity is
not known. The challenge is therefore to select nonlinear projectors
$\mathscr{A}_{l}$ that can enable the network to capture the nonlinear
dynamics of the GS system from $y_{1}$ measurements alone. First we try
monomial projectors, such as $y_{1}^{2}$, $y_{1}^{3}$, etc., as used in new-
generation-ESNs [32] and often used to approximate nonlinearity in turbulent
systems [33]. (GS$P^{n}$ refers to an ESN with nonlinear projectors
$\mathscr{A}_{l}=y_{1}^{l+1},\;l=1,2,...,n-1$.) We note from Fig. 3 (b) that
the growth of $\overline{\mathcal{E}}$ slows down slightly with increasing
number of monomials but $\overline{\mathcal{E}}$ does not saturate to
$\sqrt{2}$ for any of the GS$P^{n}$ networks. Fig. 3 (a) shows poor quality of
the predictions from the GS$P^{n}$ networks. The statistics here is calculated
based on the closed-loop predictions from $0\leq t\leq 25000$. We find that
adding more monomials or reservoir units does not significantly improve the
quality of predictions from the GS$P^{n}$ networks (results not shown).

Next, we try Gaussian RBFs as nonlinear projectors. This choice is inspired by
the robust approximation properties of RBFs [34] on which RBF networks are
based [35, 36]. For each Gaussian RBF projector
$\mathscr{A}=\exp\left(-\gamma\left(y_{1}(x,t)-y_{c}(x)\right)^{2}\right)$,
there are two parameters: (i) a shape factor $\gamma$ and (ii) a space-
dependent centre $y_{c}$. As in RBF networks, we choose the centre in an
unsupervised manner and optimize for $\gamma$ via a trial-and-correction
procedure. In the GS system, $y_{1}(x,t)$ varies from 0 to 1, and the degree
of nonlinearity increases from 1 to 0. We therefore choose $y_{c}(x)=0$ as the
centre ($y_{c}(x)=1$ is an equally good choice). Fig. 3 shows that the
GS$G_{rbf}$ network, which uses a Gaussian RBF projector
$\exp\left(-2y_{1}^{2}\right)$, performs noticeably better, both in terms of
short-term accuracy as well as the statistics of the predictions, as compared
to any of the GS$P^{n}$ networks. The results are further improved,
particularly in terms of the statistics, when two Gaussian RBF projectors
($\exp\left(-2y_{1}^{2}\right)$ and $\exp\left(-3y_{1}^{2}\right)$) are used
in the GS$2G_{rbf}$ network.

Fig. 3 (c) shows the VPT from the GS$P^{5}$ and GS$2G_{rbf}$ networks, which
highlights the significant improvement in the predictive accuracy when
Gaussian RBFs are used instead of monomials as nonlinear projectors. The small
VPT values from all the networks indicate the difficulty in obtaining accurate
predictions from partial observations. This difficulty can be further
appreciated from Fig. 6 in which the governing equations are used for
obtaining the predictions.

##  5 Relevance of nonlinear approximation properties of nonlinear projectors

We explore the connection between the nonlinear approximation properties of
the projector functions and the performance of the corresponding ESNs. Fig. 4
(a) shows the variation of the nonlinear term $y_{1}y_{2}^{2}(y_{1})$, its
mean and mid-range values ($m$) with the observed variable $y_{1}$ in the GS
system. Here $m$ is the mean of the maximum and minimum values of
$y_{1}y_{2}^{2}$ at a given $y_{1}$. The error in approximating $m$ by a
nonlinear projector function ($\mathscr{A}$) is given as

$E_{mid}=\min_{a_{1},a_{2}}\sqrt{\langle\left(a_{1}+a_{2}\mathscr{A}-m\right)^{2}\rangle}/\sqrt{\langle
m^{2}\rangle},$ (4)

where $\mathscr{A}=\exp(-\gamma y_{1}^{2})$ for Gaussian RBFs (the
corresponding $E_{mid}$ is shown in subfigure (b)) and $\mathscr{A}=y_{1}^{p}$
for monomials (the corresponding $E_{mid}$ is shown in subfigure (c)). There
are three points to note here. First, the value of $\gamma$ at which
GS$G_{rbf}$ performs best coincides with that at which $E_{mid}$ is minimum
(i.e. $\gamma\approx 2$), demonstrating that a good nonlinear projector does
indeed approximate the system nonlinearity. Second, minimization of the error
with respect to the mid-range value (rather than the mean value) indicates
that a good nonlinear projector minimizes the chances of maximum error (i.e.
$L_{\infty}$-error) [37]. Third, $E_{mid}$ is lower for the Gaussian RBFs than
for the monomials, thus explaining the superior performance of the GS$G_{rbf}$
and GS$2G_{rbf}$ networks over the GS$P^{n}$ networks in Fig. 3. Such
statistical interpretation indicates that RBFs can be useful in other RNNs as
well.

Figure 4: Nonlinear terms (a, d) and errors in approximating their mid-range
values via Gaussian RBFs (b, e) and monomials (c, f) for the GS system (a-c)
and the KS system (d-f). For the GS system, the Gaussian RBFs have relatively
lower error, thus explaining the superior performance of the networks with the
Gaussian RBF projectors. For the KS system, even when the nonlinearity is a
monomial, the Gaussian RBF can approximate the nonlinear term well.

Furthermore, the robustness of RBFs is evidenced by their ability to
satisfactorily approximate most types of nonlinearities. For example,
$\mathscr{A}=y^{2}$ is the best projector for the KS system for which the
nonlinear term itself is $y^{2}$. However, Gaussian RBFs can still approximate
$y^{2}$ nonlinearity well (see Fig 4 (e)) and the corresponding ESN
(KS$G_{rbf}$) predicts the KS system relatively well (see Fig. refKSpreds
(b-c)). Here we have only shown the results for Gaussian RBFs; the results are
similarly good when multiquadric RBFs are used (see C).

##  6 Comparison between model-free and governing-equations-based predictions

###  6.1 Application to Kuramoto–Sivashinsky system with sparse measurements

Figure 5: (a) Sparsely observed signals of the KS system and their predictions
from the governing equations and the KS$y^{2}$ network ($d_{r}=2000$ and
$N_{g}=1$). (b) The VPT of predictions at different levels of sparsity. The
observed signal becomes increasingly sparse as $d_{u}$ decreases (or
$dX\approx 128/d_{u}$ increases). (c) The VPT for the governing-equations-
based predictions decreases rapidly with increasing $dX$, varying with the
inverse of the square-root of the root-mean-square-error in the initial
condition ($e(dX)$). (d) The VPT for the KS$y^{2}$ predictions are relatively
unaffected as long as the autocorrelation between observations at neighbouring
spatial locations ($C_{x}(dX)$) remains close to one.

We now consider a situation where only sparse measurements of the KS system
are available. We compare the predictive performance of the governing-
equations-based time-marching with that of the KS$y^{2}$ network. We take the
ground truth as the original simulation with $N=128$ gridpoints. The observed
signal ($y_{obs}(x,t)$), shown in the middle column of Fig. 5 (a), consists of
measurements at $d_{u}$ equispaced points at distance $dXL/N$. This implies
that the observations become increasingly sparse with decreasing $d_{u}$ (or
increasing $dX$).

The governing-equations-based predictions are obtained by time-marching from
the initial condition at $t=0$. The initial condition, say $y_{int}(x)$, is
obtained by linearly interpolating $y_{obs}(x,t=0)$ from the observed space
(consisting of $d_{u}$ points) to the numerical grid (consisting of $N=128$
points). The predictions and the VPT for different levels of sparsity in the
observations are shown in Fig. 5 (a) and (b), respectively. The accuracy of
the predictions decreases with increasing sparsity as indicated by the
reduction in VPT. Fig. 5 (c) shows that the reduction in VPT is proportional
to $\sqrt{e(dX)^{-1}}$, where
$e(dX)=\sqrt{\langle\left(y(x,t=0)-y_{int}(x,t=0)\right)^{2}\rangle}$ is the
root-mean-square error in the initial condition. We note that the error in the
initial conditions could be reduced by using other interpolation schemes or
data assimilation [38, 39]. However, the trend would remain the same because
an interpolation scheme cannot entirely compensate for the missing information
at smaller scales. We also note that the system parameters and boundary
conditions are known perfectly here, which often contribute to error in the
governing-equations-based predictions.

The predictions from the KS$y^{2}$ network are obtained by training the
network with the sparsely observed data, i.e. the network never sees the
original data. The predictions and VPT from the KS$y^{2}$ network for
different levels of sparsity in the observations are shown in Fig. 5 (a) and
(b), respectively. We find that the predictions from the RC-RNN are
insensitive to the discretization or interpolation error; a similar level of
accuracy is maintained even when $d_{u}$ is reduced from $N$ to $N/4$. Fig. 5
(d) shows that VPT from KS$y^{2}$ is proportional to $\sqrt{C_{x}(dX)}$, where
$C_{x}(dX)$ is the mean autocorrelation between $y(x)$ and $y(x+dXL/N)$
signals. The results show that information about the well correlated elements
of the state-vector can be learned from the time history of partial
measurements. This suggests that RC-RNNs can unravel the space-time dynamics
in agreement with Takens’ embedding theorem [40], which has recently been
suggested to be the case also for data-driven discovery of the governing
equations from partial measurements [41].

The results show that RC-RNNs can provide excellent predictions up to moderate
levels of sparsity even when the governing equations may suffer from high
discretization errors. The accuracy of both methods can be improved, such as
by using better interpolation schemes for the governing equations or by
increasing the reservoir dimension for the RC-RNN. The key point, however, is
that RC-RNNs are less affected by the sparsity as compared to governing-
equations-based time-marching. When the sparsity is high such that the
spatially adjacent measurements are beyond the first zero in the
autocorrelation, we find that the RC-RNNs perform poorly. We expect the same
limitations for other data-driven methods as well

###  6.2 Application to Gray–Scott system with partial and noisy measurements

Figure 6: (a) Partially observed noisy signals for the GS system and their
predictions from the governing equations and the GS$2G_{rbf}$ network (of size
$d_{r}=1500$ and $N_{g}=4$). (b) VPT for the predictions from the governing
equations and the GS$2G_{rbf}$ network. A reduction in predictive accuracy
occurs with increasing noise level. The GS$2G_{rbf}$ network, however, is
relatively less affected by the noise. The small values of VPT from both
methods demonstrate the challenge of forecasting chaotic systems with partial
observability.

We now consider another practical scenario in which only partial and noisy
observables are available. We take the ground truth as the original simulation
of the GS system. The concentration of species one with measurement noise, say
$y_{1obs}$, is the observable (shown in the middle column of Fig. 6 (a)). The
measurement noise is additive, white in space and time, with uniform
distribution between $[-\xi\sigma_{y_{1}},\xi\sigma_{y_{1}}]$, where
$\sigma_{y_{1}}$ is the standard deviation of $y_{1}$ and $\xi$ is the noise
level. Similar to the sparse measurement cases, the governing-equations-based
predictions are obtained by time-marching from the initial conditions that are
calculated from the noisy measurements of $y_{1}$ at $t=0$, while RC-RNN-based
predictions are obtained by training the GS$2G_{rbf}$ network with the noisy
measurements.

The governing equations, their parameters and the boundary conditions are all
known perfectly, but the initial condition for $y_{2}$ is now estimated from
the noisy measurements of $y_{1}$. This introduces high levels of error in the
initial condition and, consequently, the accuracy of the governing-equations-
based predictions drops drastically with increasing noise-level as shown in
Fig. 6. The RC-RNNs are fundamentally different in that they predict the
future by learning patterns in the time history. By design, therefore, they
are robust to small levels of measurement noise. The results in Fig. 6 show
that although the accuracy of the RC-RNN-based predictions also decrease with
increasing noise level, they do not suffer as much as the governing-equations-
based predictions.

##  7 Conclusions

The main conclusions of this paper are threefold. The first two concern the
explicitly added nonlinearity in RC-RNNs. Fig. 2 and Refs. [17, 15] show it to
be essential to explicitly add nonlinearities in RC-RNNs to enable them to
capture complex system nonlinearities. Firstly, we find it more
computationally efficient to add the nonlinearity before the reservoir, such
as in the KS$y^{2}$ network, than after the reservoir, such as in the
KS$r^{2}$ network. Adding the nonlinearity before the reservoir preconditions
the network towards the system dynamics and hence approximate the system
nonlinearities better. Secondly, we find that the use of RBFs as nonlinear
projectors (to explicitly add nonlinearities in RC-RNNs) enables the networks
to capture the system nonlinearities even with partial measurements and/or
without knowing the governing equations. We show that the robust approximation
properties of RBFs, which makes them suitable for nonlinear approximation in
numerical schemes and RBF networks [35, 34, 36], also makes them a good choice
as nonlinear projectors in the RC-RNNs proposed here. The best nonlinear
projectors approximate the mid-range value of the system nonlinearities, thus
minimize the $L_{\infty}$-error. This statistical interpretation enhances the
viability of using RBFs in other RNNs for turbulence prediction.

The third conclusion is on viability of RC-RNNs as an alternative to
conventional methods in practical situations. Observations in natural and
engineering systems can often have missing data, coarse resolution and noise;
thus several ML methods have been developed to extract information from
limited and noisy measurements [6, 42]. Prediction of turbulent systems,
however, is extremely sensitive to any estimation errors. This is evident from
the results in Figs. 5 and 6, which show that the governing equations, even
with perfectly known boundary conditions and parameters, quickly lose
predictive accuracy as the measurements become sparse, incomplete and noisy.
We find that the predictive performance of RC-RNNs is relatively less affected
by sparsity and noise in the measurements, thus can potentially outperform the
governing equations in such cases.

In summary, we have demonstrated the capabilities of our proposed RC-RNN
framework to forecast (i) partially observable spatiotemporally chaotic
systems with no knowledge of the form of system nonlinearity and (ii) sparsely
observable systems or partially observable systems with measurement noise for
which even the governing equations become inaccurate. No recourse to the
governing equations is required, thus making our RC-RNN framework truly model-
free. We emphasize that the partial and sparse observability in our work is
fundamentally different from the model-reduction or partial observability
discussed in recent studies on turbulence forecasting [12, 14]. In those
studies, the extraction of dominant system modes may require noise-free
measurements of the full state-vector. Finally, we also highlight the
difficulty in predicting noisy and partially observable systems, indicating
that much larger reservoirs than those used in this study may be required for
useful predictions in real-world applications, such as weather forecasting.
Understanding the relationship between the system dynamics and hyperparameters
of RC-RNNs and the uncertainty quantification of the predictions would be
vital to unlocking the full potential of applying RC-RNNs to practical
turbulent systems.

##  8 Acknowledgments

We wish to thank Prof. Luca Biferale for his feedback during the preparation
of this manuscript. This work was funded by the National Natural Science
Foundation of China (Grant Nos. 12002147, 12050410247 and 11988102), the
Shenzhen Science and Technology Program (Grant No.
KQTD20180411143441009 ) and the Department of Science and Technology of
Guangdong Province (2020B1212030001). We acknowledge support from the Centers
for Mechanical Engineering Research and Education at MIT and SUSTech, as well
as from the Center for Computational Science and Engineering at SUSTech.. LLKB
gratefully acknowledges the funding from the Research Grants Council of Hong
Kong (Projects 16210419, 16200220, 16215521).

## References

  * [1] G. Hinton, L. Deng, D. Yu, G. Dahl, A. R. Mohamed, N. Jaitly, A. Senior, V. Vanhoucke, P. Nguyen, T. Sainath, B. Kingsbury, Deep neural networks for acoustic modeling in speech recognition: The shared views of four research groups, IEEE Signal Processing Magazine 29 (2012) 82–97.  doi:10.1109/MSP.2012.2205597.
  * [2] M. Alber, A. B. Tepole, W. R. Cannon, S. De, S. Dura-Bernal, K. Garikipati, G. Karniadakis, W. W. Lytton, P. Perdikaris, L. Petzold, E. Kuhl, Integrating machine learning and multiscale modeling—perspectives, challenges, and opportunities in the biological, biomedical, and behavioral sciences, Digital Medicine 2 (12 2019).  doi:10.1038/s41746-019-0193-y.
  * [3] X. Shi, Z. Gao, L. Lausen, H. Wang, D.-Y. Yeung, W.-K. Wong, W.-C. Woo, Deep learning for precipitation nowcasting: A benchmark and a new model, 2017, pp. 5622–5632.
URL http://arxiv.org/abs/1706.03458

  * [4] Z. Y. Wan, T. P. Sapsis, Reduced-space gaussian process regression for data-driven probabilistic forecast of chaotic dynamical systems, Physica D 345 (2017) 40–55.  doi:10.1016/j.physd.2016.12.005.
URL
http://arxiv.org/abs/1611.01583http://dx.doi.org/10.1016/j.physd.2016.12.005

  * [5] Z. Y. Wan, P. Vlachas, P. Koumoutsakos, T. Sapsis, Data-assisted reduced-order modeling of extreme events in complex dynamical systems, PLoS ONE 13 (5 2018).  doi:10.1371/journal.pone.0197704.
  * [6] M. Reichstein, G. Camps-Valls, B. Stevens, M. Jung, J. Denzler, N. Carvalhais, Prabhat, Deep learning and process understanding for data-driven earth system science, Nature 566 (2019) 195–204.  doi:10.1038/s41586-019-0912-1.
  * [7] K. Kashinath, M. Mustafa, A. Albert, J. L. Wu, C. Jiang, S. Esmaeilzadeh, K. Azizzadenesheli, R. Wang, A. Chattopadhyay, A. Singh, A. Manepalli, D. Chirila, R. Yu, R. Walters, B. White, H. Xiao, H. A. Tchelepi, P. Marcus, A. Anandkumar, P. Hassanzadeh, Prabhat, Physics-informed machine learning: Case studies for weather and climate modelling, Phil. Trans. R. Soc. A 379 (2021) 20200093.  doi:10.1098/rsta.2020.0093.
  * [8] O. Watt-Meyer, N. D. Brenowitz, S. K. Clark, B. Henn, A. Kwa, J. McGibbon, W. A. Perkins, C. S. Bretherton, Correcting weather and climate models by machine learning nudged historical simulations, Geophysical Research Letters 48 (8 2021).  doi:10.1029/2021GL092555.
  * [9] J. Pathak, S. Subramanian, P. Harrington, S. Raja, A. Chattopadhyay, M. Mardani, T. Kurth, D. Hall, Z. Li, K. Azizzadenesheli, P. Hassanzadeh, K. Kashinath, A. Anandkumar, Fourcastnet: A global data-driven high-resolution weather model using adaptive fourier neural operators, arXiv.2202.11214 (2 2022).
URL http://arxiv.org/abs/2202.11214

  * [10] C. E. Leith, R. H. Kraichnan, Predictability of turbulent flows, J. Atmos. Sci. 29 (1972) 1041–1058.
  * [11] E. N. Lorenz, The essence of chaos (1993).
  * [12] P. R. Vlachas, W. Byeon, Z. Y. Wan, T. P. Sapsis, P. Koumoutsakos, Data-driven forecasting of high-dimensional chaotic systems with long short-term memory networks, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 474 (5 2018).  doi:10.1098/rspa.2017.0844.
  * [13] Z. Li, N. Kovachki, K. Azizzadenesheli, B. Liu, K. Bhattacharya, A. Stuart, A. Anandkumar, Fourier neural operator for parametric partial differential equations, arXiv.2010.08895 (10 2020).
URL http://arxiv.org/abs/2010.08895

  * [14] P. R. Vlachas, G. Arampatzis, C. Uhler, P. Koumoutsakos, Multiscale simulations of complex systems by learning their effective dynamics, Nature Machine Intelligence 4 (2022) 359–366.  doi:10.1038/s42256-022-00464-w.
  * [15] A. Chattopadhyay, P. Hassanzadeh, D. Subramanian, Data-driven predictions of a multiscale lorenz 96 chaotic system using machine-learning methods: Reservoir computing, artificial neural network, and long short-term memory network, Nonlinear Processes in Geophysics 27 (2020) 373–389.  doi:10.5194/npg-27-373-2020.
  * [16] P. R. Vlachas, J. Pathak, B. R. Hunt, S. T. P., M. Girvan, E. Ott, P. Koumoutsakos, Backpropagation algorithms and reservoir computing in recurrent neural networks for the forecasting of complex spatiotemporal dynamics, Neural Networks 126 (2020) 191–217.
URL https://www.elsevier.com/open-access/userlicense/1.0/

  * [17] J. Pathak, B. Hunt, M. Girvan, Z. Lu, E. Ott, Model-free prediction of large spatiotemporally chaotic systems from data: A reservoir computing approach, Physical Review Letters 120 (1 2018).  doi:10.1103/PhysRevLett.120.024102.
  * [18] V. Yakhot, Large-scale properties of unstable systems governed by the Kuramoto-Sivashinksi equation, Physical Review A 24(1) (1981) 642–643.
  * [19] M. C. Cross, P. C. Hohenberg, Pattern formation outside of equilibrium, Reviews of Modern Physics 65 (1993) 851–1112.  doi:10.1103/RevModPhys.65.851.
  * [20] H. Jaeger, H. Haas, Harnessing nonlinearity: Predicting chaotic systems and saving energy in wireless communication, Science 304 (2004) 78–81.
  * [21] M. Lukoševičius, A practical guide to applying echo state networks (2012).
  * [22] V. Bratanov, F. Jenko, D. R. Hatch, M. Wilczek, Nonuniversal power-law spectra in turbulent systems, Phys. Rev. Lett. 111 (2013) 075001.
  * [23] P. Manneville, Liapounov exponents for the kuramoto-sivashinsky model (1985).
  * [24] A. K. Kassam, L. N. Trefethen, Fourth-order time-stepping for stiff pdes, SIAM Journal on Scientific Computing 26 (2005) 1214–1233.  doi:10.1137/S1064827502410633.
  * [25] G. Boffetta, M. Cencini, M. Falcioni, A. Vulpiani, Predictability: a way to characterize complexity, Physics Reports 356 (2002) 367–474.
  * [26] A. J. Majda, N. Chen, Model error, information barriers, state estimation and prediction in complex multiscale systems, Entropy 20 (8 2018).  doi:10.3390/e20090644.
  * [27] T. DelSole, Predictability and information theory. part i: Measures of predictability, J. Atmos. Sci. 61 (2004) 2425–2440.
  * [28] I. J. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozair, A. Courville, Y. Bengio, Generative adversarial nets, Curran Associates, Inc., 2014.
URL http://www.github.com/goodfeli/adversarial

  * [29] H. Kim, J. Kim, S. Won, C. Lee, Unsupervised deep learning for super-resolution reconstruction of turbulence, J. Fluid Mech. 910 (2021) A29.  doi:10.1017/jfm.2020.1028.
URL http://arxiv.org/abs/2007.15324http://dx.doi.org/10.1017/jfm.2020.1028

  * [30] Y. Nishiura, D. Ueyama, Spatio-temporal chaos for the gray-scott model, Physica D 150 (2001) 137–162.
  * [31] T. A. Driscoll, N. Hale, L. N. Trefethen, Chebfun Guide, Pafnuty Publications, 2014\.
  * [32] D. J. Gauthier, E. Bollt, A. Griffith, W. A. Barbosa, Next generation reservoir computing, Nature Communications 12 (12 2021).  doi:10.1038/s41467-021-25801-2.
  * [33] S. L. Brunton, J. L. Proctor, J. N. Kutz, W. Bialek, Discovering governing equations from data by sparse identification of nonlinear dynamical systems, Proceedings of the National Academy of Sciences of the United States of America 113 (2016) 3932–3937.  doi:10.1073/pnas.1517384113.
  * [34] M. D. Buhmann, Radial basis functions, Acta Numerica 9 (2000) 1–38.
  * [35] F. Girosi, T. Poggio, Networks and the best approximation property, Biol. Cybern. 63 (1990) 169–176.
  * [36] R. Zemouri, D. Racoceanu, N. Zerhouni, Recurrent radial basis function network for time-series prediction, Engineering Applications of Artificial Intelligence 16 (2003) 453–463.  doi:10.1016/S0952-1976(03)00063-0.
  * [37] A. Giloni, M. Padberg, Alternative methods of linear regression, Mathematical and Computer Modelling 35 (2002) 361–374.
  * [38] J. S. Whitaker, G. P. Compo, J. N. Thépaut, A comparison of variational and ensemble-based data assimilation systems for reanalysis of sparse observations, Monthly Weather Review 137 (2009) 1991–1999.  doi:10.1175/2008MWR2781.1.
  * [39] J. Brajard, A. Carrassi, M. Bocquet, L. Bertino, Combining data assimilation and machine learning to emulate a dynamical model from sparse and noisy observations: A case study with the lorenz 96 model, Journal of Computational Science 44 (7 2020).  doi:10.1016/j.jocs.2020.101171.
  * [40] F. Takens, Detecting strange attractors in turbulence (1981).
  * [41] J. Bakarji, K. Champion, J. N. Kutz, S. L. Brunton, Discovering governing equations from partial measurements with deep delay autoencoders, arXiv:2201.05136 (1 2022).
URL http://arxiv.org/abs/2201.05136

  * [42] M. Raissi, A. Yazdani, G. E. Karniadakis, Hidden fluid mechanics: Learning velocity and pressure fields from flow visualizations, Science 367 (2020) 1026–1030.
URL http://science.sciencemag.org/

##  Appendix A Details of the reservoir networks

In an echo-state network (ESN), $W_{in}\in\mathbb{R}^{d_{r}\times d_{u}}$
connects the input vector $\mathbf{u}(t)\in\mathbb{R}^{d_{u}}$ to the
reservoir nodes. Following [17], each reservoir node is connected to a single
element from the input vector, while each element of the input vector is
connected to $d_{r}/d_{u}$ reservoir nodes. Connecting each reservoir node by
a single input element ensures that uncorrelated input elements, such as from
two spatially uncorrelated locations, do not get mixed in the reservoir. The
non-zero elements of $W_{in}$ are randomly chosen from a normal distribution
in [-1, 1]. These elements are then scaled by $\sigma$, which can be a scalar
or a vector. In most of the networks used in the present study, $\sigma$ is a
scalar. In KS$G_{rbf}$, GS$P^{n}$, GS$G_{rbf}$ and GS$2G_{rbf}$ networks,
$\sigma\in\mathbb{R}^{L+1}$ (where $L$ is the number of nonlinear projectors)
such that input from each nonlinear projector is scaled independently.

The adjacency matrix $A\in\mathbb{R}^{d_{r}\times d_{r}}$ determines the
internal connections between the reservoir nodes. The elements of $A$ are also
randomly chosen from a normal distribution in [-1, 1], they are then rescaled
such that the largest eigenvalue of $A$ (also known as the spectral radius) is
equal to $\rho$. It is sometimes recommended to make $A$ sparse by randomly
fixing most of its elements to be equal to 0. However, we do not find any
advantage in doing that in our case.

The elements of $W_{out}\in\mathbb{R}^{d_{v}\times d_{r}}$ are optimized to
minimize the squared-error between the output $\mathbf{v}(t)$ and the target
signal $\widehat{\mathbf{v}}(t)$. The optimized elements are obtained by
performing the linear regression with Tikhonov regularization:

$W_{out}=\widehat{\mathbf{V}}\mathbf{R}^{T}\left(\mathbf{R}\mathbf{R}^{T}+\beta\mathbf{I}\right)^{-1},$
(5)

where $\widehat{\mathbf{V}}=\left[\widehat{\mathbf{v}}(-T+N_{tran}\Delta
t),\widehat{\mathbf{v}}(-T+N_{tran}\Delta t+\Delta
t),...,\widehat{\mathbf{v}}(-\Delta t)\right]$,
$\mathbf{R}=\left[\mathbf{r}(-T+N_{tran}\Delta t),\mathbf{r}(-T+N_{tran}\Delta
t+\Delta t),...,\mathbf{r}(-\Delta t)\right]$, $\beta=1e^{-5}$ is the
regularization parameter and $\mathbf{I}$ is an identity matrix. The vectors
$\widehat{\mathbf{v}}(-T)$, $\widehat{\mathbf{v}}(-T+\Delta t)$, …,
$\widehat{\mathbf{v}}(-\Delta t)$ and $\mathbf{r}(-T)$, $\mathbf{r}(-T+\Delta
t)$, …, $\mathbf{r}(-\Delta t)$ are the target signal and reservoir states
during the training period, i.e. from $t=-T$ to 0, $N_{tran}=50$ is the number
of transient steps discarded during the training and $\Delta t$ is the
discrete time step of the network.

##  Appendix B Calculations of the mean prediction error

There are two sources of randomness in the calculations of the prediction
error (${\mathcal{E}(t)}$): (i) the initial condition of the system and (ii)
the random seed used when drawing the elements of $W_{in}$ and $A$ matrices.
Once the network is trained using data from $t=-T$ to 0, we calculate the
predictions by reinitializing the network for 200 initial conditions. Given
that ESNs satisfy the fading memory property, we only need a finite sequence
of past inputs, say from $t_{r}-N_{ReInit}\Delta t$ to $t_{r}-\Delta t$, to
obtain the predictions starting from $t_{r}$. We find that a finite sequence
of $N_{ReInit}=10$ past inputs is sufficient to reinitialize the ESNs in our
study. During this reinitialization, the activation state is set to zero, i.e.
$\mathbf{r}(t_{r}-N_{ReInit}\Delta t)=0$, and is then evolved using the inputs
$\mathbf{u}(t-N_{ReInit}\Delta t)$, $\mathbf{u}(t_{r}-N_{ReInit}\Delta
t+\Delta t)$, …, $\mathbf{u}(t_{r}-\Delta t)$. Such calculations are repeated
for 10 random seeds. The mean prediction error ($\overline{\mathcal{E}(t)}$)
is obtained by averaging the prediction error (${\mathcal{E}(t)}$) over 200
initial conditions and 10 random seeds, i.e. for 2000 calculations.

The mean prediction error ($\overline{\mathcal{E}(t)}$) presented in the paper
is approximately optimized over all the hyperparameters for each network. The
main hyperparameters in a network of given size $d_{r}$ and the number of
parallel units $N_{g}$ are: (i) the interaction length $N_{i}$, (ii) the
spectral radius $\rho$ of $A$ and (iii) the wieght $\sigma$ of input elements.
When nonlinear projector functions are used, there are additional
hyperparameters associated with the projector functions.

##  Appendix C Different radial basis functions

In the main text, we only used Gaussian radial basis functions (RBFs) as
nonlinear projector for which we fixed the center at $y_{c}(x)=0$ to predict
the Gray-Scott system for which only concentration of one species,
$y_{1}(x,t)$, is observable. Figure 2 shows that using Gaussian RBF with
center at $y_{c}=1$ and multiquadric RBF are equally good options, thus
showing the robustness of RBFs as nonlinear projectors. A multiquadric RBF is
given as $\sqrt{1+\gamma\left(y_{1}-y_{c}\right)^{2}}$, where $\gamma$ is the
shape factor and $y_{c}$ is the center.

Figure 7: (a) The observed signal $y_{1}$ of the Gray-Scott system and (b-d)
its predictions from different RC-RNNs. The networks GS$G_{rbf}$ ($y_{c}=0$)
uses $\exp(-2y_{1}^{2})$, GS$G_{rbf}$ ($y_{c}=1$) uses $\exp(-2(y_{1}-1)^{2})$
and GS$M_{rbf}$ ($y_{c}=0$) uses $\sqrt{1+3y_{1}^{2}}$ as nonlinear
projectors.