# A note on the space of evolutionary operators in population genetics and
folding dynamics

Sergio Lukić
Corresponding address: Simons Center for Systems Biology, School of Natural
Sciences, Institute for Advanced Study, Princeton NJ 08540, USA.
<EMAIL_ADDRESS>
###### Abstract.

Discrete dynamical systems defined by the iteration of a polynomial map of the
unit simplex to itself appear in the context of population genetic systems
evolving under mutation, recombination and weak selection. Although
exceptional progress has been made in finding particular solutions to these
systems, our knowledge of the general properties of the space of all possible
dynamical systems of this kind is still limited. We prove that the space of
bounded-degree polynomial maps of the unit simplex to itself is a compact and
convex subset of a Euclidean space. We provide an explicit characterization of
such a space and of its boundary. A special class of maps in the boundary, the
folding maps, which generalize the logistic map for any dimension and degree
are defined and constructed. Finally, we use numerical methods to study the
ergodic and mixing properties of maps in the neighborhood of several of these
folding maps.

###### Key words and phrases:

Mathematical Population Genetics, Discrete Dynamical Systems, Real Algebraic
Geometry

##  1\. Introduction

Growing access to genome data generated by next-generation sequencing
technologies is revitalizing the field of mathematical population genetics. In
particular, the problem of determining the simplest historical scenarios that
have likely shaped patterns of genetic variation is motivating the development
of new mathematical and computational tools that can get the most out of these
large datasets. For instance, evolutionary biologists are often interested in
ascertaining the genealogical relationships that underlie a sample of
homologous DNA sequences. In this case, one can model the background
demographic history of the populations from which the samples were taken (e.g.
divergence between populations, migratory flows between populations, changes
in the population sizes, etc.), together with the effects of mutation,
recombination and natural selection. Using these kinds of models,
investigators can infer the simplest scenarios that best fit the observed
patterns of genetic variation and the most probable genealogies that underlie
a given set of DNA sequences [SH13].

Some of the mathematical methods that are used in these applications consist
of:

  * i)

Computations of the likelihood function by means of integrals defined over the
space of all possible genealogies associated with the sample [Fel88].

  * ii)

Sequentially Markovian approximations. These methods capture the complexities
of the genealogical relationships by means of Hidden Markov Models defined as
stochastic processes along the DNA sequences [MC05].

  * iii)

Computations of Allele Frequency Spectra (AFS) by means of random and/or
deterministic dynamical systems defined on the space of population frequencies
[LH12]. In this case, rather than working directly with DNA sequence data, one
uses summary statistics of the data defined by the frequencies of the
alternative forms of the genes in the sample.

Recent progress in the development of AFS-based techniques in addition to
their computational efficiency and mathematical elegance make this approach
particularly appealing.

Allele frequency spectra summarize the distribution of counts of alternative
forms of genes (alleles) that are observed in finite samples of DNA sequences.
This distribution is a function of the population frequencies of the
corresponding gene variants, or alleles. It is on this space of population
allele frequencies in which the dynamical systems that model the evolutionary
process is defined. The space can be represented as the direct product of unit
simplices. The dimensionality of the simplices depends on the number of
alleles associated with the genes, while the number of simplices in the direct
product corresponds to the number of populations. For simplicity, here we
consider one single population. If the number of variants of a gene in the
population is finite and equal to $n+1$, we can enumerate each allele by an
integer number $i$ and denote by $x_{i}$ the proportion of individuals that
share this particular form of the gene. The genetic composition of the
population $\\{x_{i}\\}_{i=1}^{n+1}$ is then described by a point in the
$n$-dimensional unit simplex. In the absence of genetic drift111Genetic drift
denotes the sampling noise that arises in populations with finite number of
individuals after sampling the DNA sequences that constitute each new
generation of individuals., the dynamics of the allele frequencies is
deterministic and very often can be described by means of a polynomial map
from this simplex into itself. For instance, the dynamics of a population
evolving under recurrent mutations are described by means of finite Markov
chains. Assuming that the population mates at random, the corresponding time
evolution is simply given by a linear map

(1.1) $x_{i}(t+1)=\sum_{j=1}^{n+1}M_{i\leftarrow j}x_{j}(t),$

where $M_{i\leftarrow j}$ denotes the proportion of individuals of type $j$
that mutate to type $i$ in one generation. The existence of non-zero
recombination between genes would require the addition of degree two terms in
the polynomial map [HS03], while the effects due to genetic drift can be
incorporated by means of diffusion approximations [LH12].

An understanding of the properties of maps from the simplex into itself is
motivated by the study of these dynamical systems. On the one hand, the
deterministic dynamical systems mentioned above are specified by one of these
maps. On the other hand, solving the diffusion partial differential equations
(PDEs) that arise in the presence of genetic drift, requires the use of
numerical methods that rely on the choice of a particular map from the simplex
to itself222 For instance, in the case of the Finite Element Method, often one
is interested in constructing an adaptive mesh on the simplex on which the PDE
is defined. This is because the resolution needed to numerically solve the PDE
is higher in certain parts of the simplex than in others. A natural way to
specify such an adaptive mesh consists of defining a uniform mesh on the
simplex, and an appropriate map from the simplex into itself. The choice of
the map will be such that the preimages of the subdomains where higher
resolution is needed will span higher volumes than the preimages of the
subdomains where lower resolution is needed..

A particularly interesting class of maps are the polynomial maps of bounded-
degree, which generalize the linear maps associated with finite Markov chains
to maps of arbitrary degree. In the linear case, as represented in Eq. (1.1),
the general structure of the space of linear maps can be determined easily. In
particular, as $\sum_{i=1}^{n+1}M_{i\leftarrow j}=1$ for each $j$, and every
matrix element $M_{i\leftarrow j}$ is non-negative, the space of Markov
matrices is the direct product of $n+1$ simplices of dimension $n$. We denote
such a $n(n+1)$-polytope by $\mathsf{PMaps}(\Delta^{n},1)$, and by
$\Delta^{n}$ we denote the $n$-dimensional unit simplex. As the largest
eigenvalue of $M$ is $1$ and because of the Perron-Frobenius theorem (see
[Sal97]), the dynamics of any Markov chain in the interior of
$\mathsf{PMaps}(\Delta^{n},1)$ always converges to a unique fixed point in the
$n$-simplex. This is because while one eigenvalue is $1$, the remaining
eigenvalues are in the interior of the unit disk in the complex plane. This
compares with the Markov chains that are located in the boundary of the
polytope $\mathsf{PMaps}(\Delta^{n},1)$, where the particular arrangements of
zeroes in the stochastic matrix span different faces. In particular,
stochastic matrices located in faces of decreasing dimensionality correspond
to matrices with an increasing number of zeroes. The spectrum of these
matrices consists of one eigenvalue equal to $1$, a set of eigenvalues located
in the interior of the unit disk in the complex plane, and another set of
eigenvalues located in the boundary of the unit disk. Indeed, when
$M_{i\leftarrow j}$ has the maximum number of zero entries and its determinant
is not zero, i.e. the matrix is a permutation matrix (see Appendix A for an
elementary proof of this statement), all the eigenvalues are located in the
boundary of the unit disk. The spectrum of any of these stochastic matrices
allows us to determine some dynamical properties of the corresponding Markov
chains; e.g. the existence of cycles is associated with more than one
eigenvalue in the unit circle of the complex plane (see for instance [Sal97]).

In this paper, we are able to extend some of these known results in the linear
case to the general case of the space of bounded-degree stochastic polynomial
maps from the simplex into itself. First, we demonstrate that these spaces are
compact and convex subsets of Euclidean spaces. Second, we characterize a
class of maps in the boundary of these spaces that generalize the notion of
stochastic matrix with maximum number of zeroes and nonzero determinant.
Geometrically, these polynomial maps are defined by maximizing the number of
preimages of the boundary of the simplex. Indeed, this is satisfied by a
family of maps that we introduce here, the _minimal polynomial folding maps_ ,
which, as we show below, are related to the theory of simplicial tessellations
of the Euclidean space [HW88, Ves91]. Finally, we show how these maps can be
constructed by solving certain systems of polynomial equations, we compute
several examples in one and two dimensions, and we study the dynamical
properties of these maps and of their deformations in the space of polynomial
maps.

###  1.1. Basic notation and definitions.

We denote by $\Delta^{n}\subset\mathbb{R}^{n+1}$ the standard $n$-simplex:

(1.2)
$\Delta^{n}=\left\\{x=(x_{1},x_{2},\ldots,x_{n+1})\in\mathbb{R}^{n+1}\colon\,\,\sum_{i=1}^{n+1}x_{i}=1,\,\,x_{j}\geq
0,\,1\leq j\leq n+1\right\\}.$

For simplicity, the vertices, edges, faces and other boundary components of
the $n$-simplex are denoted as its $0$-faces, $1$-faces, $2$-faces, …, and
$(n-1)$-faces. The interior of the $n$-simplex,
$\Delta^{n}\backslash\partial\Delta^{n}$, is also referred to as the $n$-face
of the simplex. All the $i$-faces except for the $0$-faces, which are closed
sets in $\mathbb{R}^{n+1}$, are open sets. Hence, any simplex can be expressed
as the disjoint union of its faces. This implies that if $\mathcal{S}$ is the
set of all the faces in $\Delta^{n}$, there exists an onto map
$F_{c}\colon\Delta^{n}\to\mathcal{S}$, the _face map_ , which assigns to each
point in the simplex its corresponding face, i.e., $x\in F_{c}(x)$.

Any polynomial map from the $n$-simplex to itself is defined by means of $n+1$
polynomials $\\{P_{i}(x_{1},x_{2},\ldots,x_{n+1})\\}_{i=1}^{n+1}$ satisfying

$\sum_{i=1}^{n+1}P_{i}(x)=1\quad\mathrm{and}\quad P_{j}(x)\geq 0,\,1\leq j\leq
n+1,$

for all $x\in\Delta^{n}$. We define the degree of a given map as the degree of
the polynomial $P_{i}(x)$ of the largest degree. This definition is ambiguous
because any pair of polynomials, $P_{i}(x)$ and
$P_{i}(x)(x_{1}+x_{2}+\cdots+x_{n+1})^{N}$ restricted to $\Delta^{n}$, are
identical for any non-negative integer $N$. To avoid this ambiguity, we assume
throughout that the degree of a polynomial $P(x)$ is the smallest degree
possible after factoring out $(x_{1}+x_{2}+\cdots+x_{n+1})$ terms. We denote
by $\mathsf{PMaps}(\Delta^{n},k)$ the space of polynomial maps of degree less
than or equal to $k$ from $\Delta^{n}$ to itself. By $\Pi_{k}^{n}$ we denote
the vector space of polynomials in $n$ variables of degree at most $k$.

A $d$-folding map from the simplex to itself is a $d$-to-$1$ and onto
continuous mapping in which the number of preimages is $d$ everywhere in the
interior of $\Delta^{n}$. A polynomial fold of $\Delta^{n}$ is a particular
instance of a folding map that is also a polynomial map.

##  2\. The geometry of the space of stochastic polynomial maps

For a given positive integer $d$, the corresponding $d$-folding polynomial
maps located in the space of polynomial maps $\mathsf{PMaps}(\Delta^{n},k)$
with the lowest degree $k$ possible, are a very special set of points in
$\mathsf{PMaps}(\Delta^{n},k)$. This will become clearer after we describe in
more detail the geometry of $\mathsf{PMaps}(\Delta^{n},k)$.

######  Theorem 1.

The space of polynomial maps of degree $k$ from the $n$-simplex to itself,
$\mathsf{PMaps}(\Delta^{n},k)$, is a compact and convex subset of the vector
space of polynomial maps of degree $k$ from $\mathbb{R}^{n}$ to itself.

_Proof_. The proof consists of an explicit construction of
$\mathsf{PMaps}(\Delta^{n},k)$, such that the properties stated in the theorem
follow from such construction. For simplicity, we work with the projection of
the simplex $\Delta^{n}\subset\mathbb{R}^{n+1}$ down to $\mathbb{R}^{n}$,
where $\Delta^{n}$ is defined by the set of inequalities

$x_{i}\geq 0,\,\,1\leq i\leq n,$ $\sum_{i=1}^{n}x_{i}\leq 1.$

Any polynomial map $f$ in $\mathsf{PMaps}(\Delta^{n},k)$ can be defined by
means of $n$ polynomials $\\{P_{i}(x)\\}_{i=1}^{n}$ of degree less than or
equal to $k$, that obey the following conditions:

(2.1) $\displaystyle P_{i}(x)\geq 0,\,\,\forall i\in[1,n],\,\forall
x\in\Delta^{n},$ (2.2) $\displaystyle 1-\sum_{i=1}^{n}P_{i}(x)\geq
0,\,\,\forall x\in\Delta^{n},$

and

(2.3) $x^{\prime}=f(x),\quad x^{\prime}_{i}=P_{i}(x),$

which states how preimages $x$ are mapped to images $x^{\prime}$ and
therefore, it completes the definition of $f$. The condition in Eq. (2.1) can
be simplified thanks to an important theorem of Pólya, which provides a
systematic process for deciding whether a given polynomial is strictly
positive on the simplex [HLP52]. Using Pólya’s result, one can show how the
polynomials that are non-negative on $\Delta^{n}$ span a convex cone in
$\Pi_{k}^{n}$, the space of polynomials in $n$ variables of degree at most
$k$. To show this, we first need to construct homogeneous representatives of
the polynomials in $\Pi_{k}^{n}$. In particular, if $P(x)$ is a polynomial in
$\Pi_{k}^{n}$, we construct its homogenous polynomial representative
$P_{H}(x)$ in $\mathbb{R}^{n+1}$ as follows. We write $P(x)$ as a sum of
monomials

$P(x)=\sum_{|\alpha|\leq
k}c_{\alpha_{1}\alpha_{2}\cdots\alpha_{n}}x^{\alpha_{1}}_{1}x^{\alpha_{2}}_{2}\cdots
x^{\alpha_{n}}_{n}$

with $|\alpha|=\sum_{i=1}^{n}\alpha_{i}$, and multiply the terms with
$|\alpha|<k$ by $(x_{1}+x_{2}+\cdots+x_{n+1})^{k-|\alpha|}$. In this way, we
obtain a homogenous polynomial of degree $k$ which is identical to $P(x)$ when
restricted to $\Delta^{n}\subset\mathbb{R}^{n+1}$: $P_{H}(x)=$

$\sum_{|\alpha|\leq
k}c_{\alpha_{1}\alpha_{2}\cdots\alpha_{n}}\left(\delta_{|\alpha|,k}+(1-\delta_{|\alpha|,k})(x_{1}+x_{2}\cdots+x_{n+1})^{k-|\alpha|}\right)x^{\alpha_{1}}_{1}x^{\alpha_{2}}_{2}\cdots
x^{\alpha_{n}}_{n},$

with $\delta_{|\alpha|,k}=1$ when $|\alpha|=k$ and $\delta_{|\alpha|,k}=0$
otherwise. Using this homogenous polynomial $P_{H}(x)$, Pólya’s theorem simply
says that if $P(x)$ is positive on $\Delta^{n}$, i.e. $P(x)>0$, $\forall
x\in\Delta^{n}$, then there exists a non-negative integer $N$ such that all
the coefficients of $(x_{1}+x_{2}+\cdots+x_{n+1})^{N}P_{H}(x)$ are positive.
If $P(x)$ is positive on $\Delta^{n}$ then $\lambda P(x)$, with
$\lambda\in\mathbb{R}^{+}$, is also positive. This implies that the
polynomials of degree $k$ with $P(x)>0$, $\forall x\in\Delta^{n}$, span an
open convex cone $\mathcal{K}_{+}(\Delta^{n},k)\subset\Pi_{k}^{n}$. It also
follows that $\dim\,\mathcal{K}_{+}(\Delta^{n},k)=\dim\,\Pi_{k}^{n}$.
Furthermore, one can use this construction to describe the set of non-negative
polynomials on the simplex. More precisely, if $P(x)$ is strictly positive on
the $n$-simplex and $Q(x)$ is non-negative for all $x\in\Delta^{n}$, then any
deformation of $Q(x)$ by $\epsilon P(x)$ makes $Q(x)+\epsilon P(x)$ an element
of $\mathcal{K}_{+}(\Delta^{n},k)$ for any positive real value of $\epsilon$.
Therefore, the closure of $\mathcal{K}_{+}(\Delta^{n},k)$ in $\Pi_{k}^{n}$,
which we denote by
$\mathcal{K}(\Delta^{n},k)=\mathrm{cl}\left(\mathcal{K}_{+}(\Delta^{n},k)\right)$,
is the set of non-negative polynomials on $\Delta^{n}$. In other words,

$\\{P(x)\in\Pi_{k}^{n}\colon\,\,P(x)\geq 0,\,\,\forall
x\in\Delta^{n}\\}=\mathcal{K}(\Delta^{n},k).$

Thus, the set of polynomials $\\{P_{i}(x)\\}_{i=1}^{n}$ that satisfy Eq. (2.1)
are elements of the cartesian product of $n$ closed convex cones
$\mathcal{K}(\Delta^{n},k)$. This product of cones is itself a closed convex
cone in $\oplus_{i=i}^{n}\Pi_{k}^{n}$, the direct sum of vector spaces of
polynomials on the simplex. We denote such a convex cone by
$\times_{i=1}^{n}\mathcal{K}(\Delta^{n},k)$.

The space $\mathsf{PMaps}(\Delta^{n},k)\subset\oplus_{i=i}^{n}\Pi_{k}^{n}$ and
its compactness property arise after we impose Eq. (2.2) on the convex cone
$\times_{i=1}^{n}\mathcal{K}(\Delta^{n},k)$. By the Heine-Borel theorem it is
sufficient to show that $\mathsf{PMaps}(\Delta^{n},k)$ is a bounded subset of
the cone
$\times_{i=1}^{n}\mathcal{K}(\Delta^{n},k)\subset\oplus_{i=i}^{n}\Pi_{k}^{n}$,
in order to demonstrate that $\mathsf{PMaps}(\Delta^{n},k)$ is compact. First,
let

$\mathrm{span}_{\lambda}\\{\lambda P_{i}(x)\\}_{i=1}^{n}$

for $\lambda\in\mathbb{R}^{+}$ positive and real, denote the ray in
$\times_{i=1}^{n}\mathcal{K}(\Delta^{n},k)$ generated by
$\\{P_{i}(x)\\}_{i=1}^{n}$. The maximum of the function
$\sum_{i=1}^{n}P_{i}(x)$ on the simplex is strictly positive for any non-zero
set of polynomials
$\\{P_{i}(x)\\}_{i=1}^{n}\in\times_{i=1}^{n}\mathcal{K}(\Delta^{n},k)$. Let us
denote such a maximum value by $S_{max}$. Now, Eq. (2.2) restricted on this
ray is equivalent to requiring the non-negativity of
$1-\lambda\sum_{i=1}^{n}P_{i}(x)$ for all $x$ in the simplex. This is only
satisfied by the segment $\\{\lambda P_{i}(x)\\}_{i=1}^{n}$ parametrized by
$\lambda\in[0,1/S_{max}]$. Similarly, for every ray in
$\times_{i=1}^{n}\mathcal{K}(\Delta^{n},k)$ only a bounded and finite
subinterval of the ray satisfies Eq. (2.2). It follows that
$\mathsf{PMaps}(\Delta^{n},k)$ is a bounded and closed subset of the convex
cone $\times_{i=1}^{n}\mathcal{K}(\Delta^{n},k)$.

The convexity of $\mathsf{PMaps}(\Delta^{n},k)$ follows from the convexity of
the simplex $\Delta^{n}$. In particular, if the image of $x$ defined by
$\\{P_{i}(x)\\}_{i=1}^{n}\in\mathsf{PMaps}(\Delta^{n},k)$ is in $\Delta^{n}$
for each $x\in\Delta^{n}$, and the same is true for a second map
$\\{Q_{i}(x)\\}_{i=1}^{n}\in\mathsf{PMaps}(\Delta^{n},k)$, then
$\\{tP_{i}(x)+(1-t)Q_{i}(x)\\}_{i=1}^{n}$ is in $\Delta^{n}$ and yields a
well-defined map in $\mathsf{PMaps}(\Delta^{n},k)$ for any
$t\in[0,1]$.$\square$

---
Figure 1.  Convex set of quadratic maps of the type $f(x)=x(A+Bx)$ on the
$1$-simplex. The horizontal axis parametrizes the $A$-coefficient and the
vertical axis parametrizes the $B$-coefficient. For pairs $\\{A,\,B\\}$
outside the convex set, the image of the quadratic map is not in the simplex.
The red dot denotes the identity map, while the black dot denotes the
polynomial two-folding map
$f_{2}(x)=\frac{1-\cos(2\arccos(1-2x))}{2}=4x(1-x).$

Our construction of $\mathsf{PMaps}(\Delta^{n},k)$ relies on Pólya’s criterion
to decide whether a given polynomial vanishes somewhere in the simplex
$\Delta^{n}$ or not. In particular, we have showed how the interior of
$\mathsf{PMaps}(\Delta^{n},k)$ consists of maps for which every defining
polynomial, $\\{P_{i}(x)\\}_{i=1}^{n}$ plus the polynomial
$1-\sum_{j}P_{j}(x)$, are strictly positive on $\Delta^{n}$, and therefore,
their corresponding images are in the interior of $\Delta^{n}$. This compares
with the boundary of $\mathsf{PMaps}(\Delta^{n},k)$, consisting of maps for
which some of the defining polynomials do vanish on some algebraic subsets of
the standard simplex. For these maps the preimages of the boundary of
$\Delta^{n}$ are not the empty set. The case of finite Markov chains is
particularly illuminating in this regard. When $k=1$, the defining polynomials
are non-negative linear forms on the simplex, while the boundary of
$\mathsf{PMaps}(\Delta^{n},1)$ consists of maps that have a non-empty preimage
of the boundary $\partial\Delta^{n}$. In an extreme limiting case, there exist
maps with non-singular Jacobian matrices that maximize the size of the
corresponding preimages of $\partial\Delta^{n}$. These maps consist of
invertible stochastic matrices that have the largest number possible of zero
entries, i.e. the permutation matrices. Geometrically, these matrices are
special vertices in the boundary of the $n(n+1)$-dimensional polytope
$\mathsf{PMaps}(\Delta^{n},1)$, among other things because they are the only
bijective maps therein (see Lemma 1 in the Appendix).

One can define nonlinear analogs of these extremal maps in
$\mathsf{PMaps}(\Delta^{n},k)$ for $k>1$. They consist of polynomial maps
$f\colon\,\Delta^{n}\to\Delta^{n}$ with a non-singular Jacobian almost
everywhere in $\Delta^{n}$, that maximize the complexity of the algebraic
subset in $\Delta^{n}$ that is the preimage of the boundary under $f$, i.e.
they maximize the complexity of
$f^{-1}\left(\partial\Delta^{n}\right)\subset\Delta^{n}$. These maps are easy
to construct in the case of the $1$-simplex. For instance, when $n=1$ and
$k=2$ the largest algebraic subset of the interval that can be mapped to its
boundary by a quadratic map (with a non-singular Jacobian almost everywhere in
$[0,1]$) consists of three points (see Fig. 1). In general, the largest
algebraic subset of $\Delta^{1}$ that can be mapped to the boundary of
$\Delta^{1}$ by a map in $\mathsf{PMaps}(\Delta^{1},k)$ consists of $k+1$
points. These maps can be constructed explicitly [HW88, Ves91]. They coincide
with the Chebyshev maps and their permutations:

(2.4) $f_{k}\colon\,[0,1]\to[0,1],\quad f_{k}(x)=\frac{1-T_{k}(1-2x)}{2}.$

Here, $T_{k}$ denotes the $k^{th}$ Chebyshev polynomial of the first kind,
which can be evaluated as $T_{k}(x)=\cos\left(k\arccos x\right)$. One can show
from the properties of the Chebyshev maps that every map $f_{k}$ is a
$k$-folding map of the $1$-simplex.

As far as we know, there do not exist previous attempts in the literature to
construct the higher dimensional analogs of these extremal maps on the unit
simplex. Although these maps are indeed harder to construct, we show below how
some particular families of extremal maps admit a simple geometrical
description. In particular, the concept of a polynomial map $f$ of fixed
degree that maximizes the complexity of
$f^{-1}\left(\partial\Delta^{n}\right)$ finds an instance in any polynomial
$d$-folding map of minimal degree. As every polynomial $d$-folding map $f$ is
continuous, onto and $d$-to-$1$ everywhere in the interior of the image
simplex, one can decompose the preimage simplex $\Delta^{n}$ as the union of
$d$ sets, $\Delta^{n}=\cup_{a=1}^{d}\delta_{a}$, such that the restriction of
the map on each subset, $f|_{\delta_{a}}\colon\,\delta_{a}\to\Delta^{n}$, is
one-to-one and onto. In this case, the preimage of the boundary is simply the
union of the boundaries of each subset $\delta_{a}$,
$f^{-1}\left(\partial\Delta^{n}\right)=\cup_{a=1}^{d}\partial\delta_{a}$.
Therefore, those $d$-folding maps $f$ located in spaces
$\mathsf{PMaps}(\Delta^{n},k)$ that have the smallest $k$ possible maximize
the complexity of $f^{-1}\left(\partial\Delta^{n}\right)\subset\Delta^{n}$
among all maps in $\mathsf{PMaps}(\Delta^{n},k)$. In the next section we
define and construct these minimal polynomial folding maps. First, we propose
a restricted notion of homotopy that allows us to group together equivalent
folding maps. For instance, we consider $4x(1-x)$ and
$2x\Theta(1/2-x)+(2-2x)\Theta(x-1/2)$ to be topologically equivalent
$2$-folding maps of the $1$-simplex (here, $\Theta$ denotes the Heaviside
function); however, we do not consider the polynomial $2$-fold $4x(1-x)$ and
the piecewise linear $3$-fold
$3x\Theta(1/3-x)+(2-3x)\Theta(x-1/3)\Theta(2/3-x)+(3x-2)\Theta(x-2/3)$ to be
equivalent. Then, we use this notion of homotopy to define the _minimal
polynomial folding maps_. Finally, we show how to construct these maps as the
zero loci of certain polynomial systems equations and compute several
examples.

##  3\. Polynomial Folding Maps

In order to be able to compare polynomial folds with non-polynomial folds,
e.g. piecewise linear folding maps, we need an appropriate criterion to
establish a topological equivalence relation between them. To that end, we
introduce a restricted notion of homotopy equivalence between folding maps as
follows.

|
---|---
Figure 2.  A piecewise linear $2$-fold of the $1$-simplex. Here, the
$1$-simplex spans the interval $[0,1]$, the horizontal axis denotes the
preimage of the map, and the vertical axis denotes its image. The standard
notion of homotopy equivalence between maps is not valid to compare folding
maps. The $2$-folding map, given by the function
$f(x)=10x\Theta(\frac{1}{10}-x)+(1-\frac{10}{9}(x-\frac{1}{10}))\Theta(x-\frac{1}{10})$,
is homotopic to the polynomial fold $g(x)=4x(1-x)$ (see left plot). The same
$2$-folding map is also homotopic to the identity map $I(x)=x$, (see right
plot).

We say that two maps, $f\colon\,\Delta^{n}\to\Delta^{n}$ and
$g\colon\,\Delta^{n}\to\Delta^{n}$, are _face-constrained homotopic_ (or _fc-
homotopic_) if and only if there exist two continuous functions,
$H\colon\,\Delta^{n}\times[0,1]\to\Delta^{n}$ and
$K\colon\,\Delta^{n}\times[0,1]\to\Delta^{n}$, from the product of the
$n$-simplex with the unit interval $[0,1]$ to the $n$-simplex such that:

  * i)

$K(x,t)\colon\,\Delta^{n}\to\Delta^{n}$ is one-to-one and onto for all
$t\in[0,1]$.

  * ii)

$\\{K(x,0),\,H(x,0)\\}=\\{x,\,f(x)\\}$ and
$\\{K(x,1),\,H(x,1)\\}=\\{y,\,g(y)\\}$ for all $x\in\Delta^{n}$.

  * iii)

The _face-constrained_ condition requires that if $F_{c}(x)$ is the face where
$x$ belongs, $x\in F_{c}(x)\subset\Delta^{n}$, and $F_{c}\left(f(x)\right)$
denotes the face where the image $f(x)$ is located, $f(x)\in
F_{c}\left(f(x)\right)\subset\Delta^{n}$, then the deformation defined by $K$
and $H$ preserves such constraints, i.e. $K(x,t)\in F_{c}(x)$ and $H(x,t)\in
F_{c}\left(f(x)\right)$ for all $t\in[0,1]$.

Note that if the condition iii) is not satisfied, one recovers the standard
notion of homotopy by choosing $K$ to be the identity map.

---
Figure 3.  A _face-constrained homotopy_ of the piecewise linear fold defined
by
$f(x)=10x\Theta(\frac{1}{10}-x)+(1-\frac{10}{9}(x-\frac{1}{10}))\Theta(x-\frac{1}{10})$.
The preimage points $x=0$ and $x=1$ are located on $0$-faces, while both of
their images are on the same $0$-face ($y=0$). The images of all the interior
points, except $x=1/10$, are also interior points. The image of the interior
point $x=1/10$ is located on the $0$-face $y=1$. These constraints are
preserved by any _fc-homotopy_. Hence, while $f(x)$ is _fc-homotopic_ to
$g(x)=4x(1-x)$, it is not _fc-homotopic_ to the identity map.

######  Definition 1.

_Minimal polynomial folding map_. Given an fc-homotopy class of folding maps
on the $n$-simplex, we say that a particular polynomial representative is
minimal if and only if the corresponding polynomial map exists in
$\mathsf{PMaps}(\Delta^{n},k)$ and there do not exist fc-homotopically
equivalent maps in $\mathsf{PMaps}(\Delta^{n},l)$, with $l<k$.

Let us assume that a given minimal polynomial folding map
$f\colon\,\Delta^{n}\to\Delta^{n}$ is defined by the polynomials
$\\{P_{i}(x)\\}_{i=1}^{n}$ and $P_{n+1}(x)=1-\sum_{j=1}^{n}P_{j}(x)$. Then,
the preimage of each of the $n+1$ facets of the simplex is an algebraic subset
of the unit simplex defined either by

$\\{x\in\Delta^{n}\colon\,P_{i}(x)=0\\},$

for the facet $y_{i}=0$, or by

$\\{x\in\Delta^{n}\colon\,1-\sum_{j=1}^{n}P_{j}(x)=P_{n+1}(x)=0\\},$

for the facet $\sum_{i=1}^{n}y_{i}=1$. We denote by
$\\{z_{i}\subset\Delta^{n}\backslash\partial\Delta^{n},\,b_{i}\subset\partial\Delta^{n}\\}$
the codimension-one preimage of the $i^{th}$ facet, decomposed in its interior
and boundary parts. In this case, $z_{i}$ is the subset of the preimage
located in the interior of the simplex and $b_{i}$ is the complementary subset
located in the boundary of the simplex. It is important to note that as
$b_{i}\subset\partial\Delta^{n}$ is the zero locus of a polynomial defined on
$\Delta^{n}$, $b_{i}$ is either the empty set, or a finite set of facets.
Using other words, $b_{i}$ cannot be composed of finite subsets within facets
because such subsets are not algebraic; $b_{i}$ can only be composed of whole
facets, because they are the only algebraic subsets of codimension one in the
boundary of the simplex, e.g. $x_{i}=0$, or $1-\sum_{j=1}^{n}x_{j}=0$. Indeed,
there exist folding maps that do not admit fc-homotopic polynomial
representatives; this is because whenever different subsets within a single
facet are mapped to distinct facets, the folding map cannot be polynomial.
This will become clearer after we show the following.

######  Theorem 2.

Let $f\in\mathsf{PMaps}(\Delta^{n},k)$ be a minimal polynomial folding map of
degree $k$. Then, each defining polynomial $P_{i}(x)$, with $1\leq i\leq n+1$,
can be factored as

(3.1) $P_{i}(x)=l_{i}(x)(q_{i}(x))^{2}m_{i}(x),$

where $l_{i}(x)$, $q_{i}(x)$ and $m_{i}(x)$ obey the following:

  * •

The zero locus of $l_{i}(x)$ in $\Delta^{n}$ is exactly
$b_{i}\subset\partial\Delta^{n}$.

  * •

The degree of $l_{i}(x)$ equals the number of facets contained in $b_{i}$.
Furthermore, the product of all the $l_{i}(x)$’s is

$\prod_{i=1}^{n+1}l_{i}(x)=(1-\sum_{j}x_{j})\prod_{i=1}^{n}x_{i}.$

  * •

The zero locus of $q_{i}(x)$ in the interior of the $n$-simplex is exactly
$z_{i}\subset\Delta^{n}\backslash\partial\Delta^{n}$.

  * •

$l_{i}(x)$ and $m_{i}(x)$ are strictly positive in the interior of the
$n$-simplex. The zero set of $m_{i}(x)$ is either empty or it is contained in
the zero set of $l_{i}(x)$.

  * •

The degree of $P_{i}(x)$ is bounded by $k$, $\deg l_{i}+2\deg q_{i}+\deg
m_{i}\leq k$.

_Proof_. The preimage of $\Delta^{n}\backslash\partial\Delta^{n}$ under $f$
consists of $d$ disconnected open subsets in the interior of $\Delta^{n}$. The
restriction of $f$ to any of these subsets gives rise to a one-to-one mapping
that has the whole interior of the simplex as codomain. The complementary set
of the union of these $d$ open subsets in $\Delta^{n}$ is
$f^{-1}(\partial\Delta^{n})$. These facts follow from the properties of the
folding maps: they are continuous and onto, for every point in the interior of
the simplex there exist $d$ different preimages, and they map the boundary
$\partial\Delta^{n}$ to a subset in the boundary of the image simplex (see
Lemma 2 in the Appendix). In addition to this, the fact that a minimal
polynomial folding map is also polynomial imposes further constraints. In
particular, the properties of the defining polynomials stated in Theorem $2$
are derived from the algebraic properties of the preimage set
$f^{-1}(\partial\Delta^{n})$. Algebraically,
$f^{-1}\left(\partial\Delta^{n}\right)$ is the union of irreducible varieties
of codimension one in $\Delta^{n}$. The set defined by the zero locus
$P_{i}(x)=0$ denotes the algebraic variety in
$f^{-1}\left(\partial\Delta^{n}\right)$ that is preimage of the $i^{th}$
facet. As such a preimage has to be the union of a finite number of facets and
a codimension one variety in the interior of the simplex, $P_{i}(x)$ can be
factored as

$P_{i}(x)=l_{i}(x)r_{i}(x)m_{i}(x).$

Here, the zero locus of $l_{i}(x)$ is the union of facets in the boundary, the
zero locus of $r_{i}(x)$ is a codimension one variety in the interior of the
simplex, and $m_{i}(x)$ is a non-negative polynomial of degree equal or less
than $k-\deg l_{i}(x)-\deg r_{i}(x)$, whose zero locus is either empty or it
is contained in the zero locus of $l_{i}(x)$. The existence of
$m_{i}(x)\in\mathcal{K}(\Delta^{n},$ $k-\deg l_{i}(x)-\deg r_{i}(x))$ follows
from the facts that $P_{i}(x)$ has a degree equal or less than $k$, that
$P_{i}(x)\geq 0$ for all $x\in\Delta^{n}$ and that the zero set of $P_{i}(x)$
is already specified by $r_{i}(x)$ and $l_{i}(x)$. Because of these
constraints, the existence of an additional polynomial factor of degree less
than or equal to $k-\deg l_{i}(x)-\deg r_{i}(x)$ is possible.

As the $i^{th}$-facet of the $n$-simplex is defined by the zero locus of the
degree-one monomial $x_{i}$ (when $i=n+1$, $x_{n+1}=1-\sum_{j=1}^{n}x_{j}$),
the lowest degree polynomial whose zero locus consists of a predefined set of
facets $b_{i}\subset\Delta^{n}$ is

$l_{i}(x)=\prod_{j\in b_{i}}x_{j}.$

If $b_{i}$ is empty then $l_{i}(x)=1$. It follows that $l_{i}(x)>0$ for all
$x$ in the interior of the simplex. Furthermore, as $f(\partial\Delta^{n})$ is
contained in $\partial\Delta^{n}$, the images of every facet in $\Delta^{n}$
are in the union of all the image facets, and therefore

$\prod_{i=1}^{n+1}l_{i}(x)$

has to be a polynomial divisible by
$(1-\sum_{j=1}^{n}x_{j})\prod_{i=1}^{n}x_{i}$. As the preimage of each facet
is the union of irreducible varieties of codimension one, each facet has to be
mapped to a whole facet. This implies that the monomial $x_{j}$ whose zero
locus defines the facet cannot appear as a factor in more than one polynomial
$l_{i}(x)$. Otherwise, the image of such $j^{th}$-facet would be a face of
codimension larger than one. It follows that $\prod_{i=1}^{n+1}l_{i}(x)$ is
not only divisible by $(1-\sum_{j=1}^{n}x_{j})\prod_{i=1}^{n}x_{i}$ but it is
a constant multiple of $(1-\sum_{j=1}^{n}x_{j})\prod_{i=1}^{n}x_{i}$. By
convention, such a constant factor is absorbed in the $m_{i}(x)$ terms, and

$\prod_{i=1}^{n+1}l_{i}(x)=(1-\sum_{j=1}^{n}x_{j})\prod_{l=1}^{n}x_{l}.$

Finally, it remains to show that the factor $r_{i}(x)$ is the square of a
polynomial: $r_{i}(x)=\left(q_{i}(x)\right)^{2}$. As the zero set $z_{i}=$
$\\{x\in\Delta^{n}\colon\,r_{i}(x)=0\\}$ is a codimension-one variety in the
interior of the simplex, and $r_{i}(x)$ is strictly positive for every
$x\in\Delta^{n}\backslash z_{i}$, it follows from the Taylor expansion of
$r_{i}(x)$ around any $x\in z_{i}$ that the gradient of $r_{i}(x)$ restricted
to $z_{i}$ has to vanish. Therefore, $z_{i}$ has to be the zero locus of a
polynomial $q_{i}(x)$ which might have negative values on a subset of the
simplex and satisfies $r_{i}(x)=(q_{i}(x))^{2}$.$\square$

A priori, the factorization in Eq. (3.1) can be used to construct any minimal
polynomial folding map. First, one fixes an fc-homotopy class of folding maps
that admits polynomial representatives. Then, the assignment of facets to
image facets defined by the fc-homotopy class is used to determine the
polynomials $\\{l_{i}(x)\\}_{i=1}^{n+1}$. Unfortunately, specifying the
remaining factors in Eq. (3.1) is a harder problem. The approach that we
follow consists of parametrizing the set of all the polynomials,
$\\{m_{i}(x,\phi^{i})\\}_{i=1}^{n+1}$ and $\\{q_{i}(x,\theta)\\}_{i=1}^{n+1}$,
that are compatible with the given class of folding maps. Then, we determine
the particular parameters $\\{\hat{\theta},\,\hat{\phi}^{i}\\}_{i=1}^{n+1}$
that define the folding map by solving a system of polynomial equations. More
precisely, we use the defining equation

(3.2)
$l_{n+1}(x)\left(q_{n+1}(x,\theta)\right)^{2}m_{n+1}(x,\phi^{n+1})=1-\sum_{i=1}^{n}l_{i}(x)\left(q_{i}(x,\theta)\right)^{2}m_{i}(x,\phi^{i}),$

to determine the parameters $\\{\hat{\theta},\,\hat{\phi}^{i}\\}_{i=1}^{n+1}$
that satisfy Eq. (3.2) for all $x\in\Delta^{n}$. In practice, this is a
tedious approach because the set of compatible polynomials usually has several
independent components and each component is a complicated semialgebraic set
in a high-dimensional vector space. Nevertheless, it is possible to identify
each component and construct the corresponding spaces of polynomials. In
particular, every component is associated with a different partition
$\\{(a_{i},b_{i})\\}_{i}$ of the degrees of $(q_{i}(x))^{2}m_{i}(x)$; i.e.
partitions of $k-\deg l_{i}(x)=2a_{i}+b_{i}$, with $a_{i}=\deg q_{i}(x)$,
$b_{i}=\deg m_{i}(x)$ and $k$ the degree of the polynomial map. Given a
particular partition $\\{(a_{1},b_{1}),$ $(a_{2},b_{2}),$
$\ldots(a_{n+1},b_{n+1})\\}$ of the degree of every polynomial factor, it is
possible to construct the space of compatible polynomial factors with these
prescribed degrees and find the solutions of Eq. (3.2).

This construction becomes more transparent in the next subsection, where we
show several examples in one and two dimensions. However, before we construct
particular examples, it is useful to state some general properties of these
families of polynomials. For instance, the family of polynomials
$\\{q_{i}(x,\theta)\\}_{i=1}^{n+1}$ can be determined by requiring the
corresponding zero loci $\\{z_{i}\subset\Delta^{n}\\}_{i=1}^{n+1}$ to be
compatible with the given class of folding maps. This allows for the
possibility of the existence of maps in the fc-homotopy class for which the
preimage of $\partial\Delta^{n}$ restricted to the interior of the simplex is
exactly defined by the set of hypersurfaces $\\{z_{i}\\}_{i=1}^{n+1}$. The
parameter space $\\{\phi^{i}\\}$ admits a simpler description. The set of
compatible polynomials $m_{i}(x,\phi^{i})$ consists of those polynomials of
degree $b_{i}$

$m_{i}(x,\phi^{i})=\sum_{|\alpha|\leq
b_{i}}\phi^{i}_{\alpha_{1}\alpha_{2}\cdots\alpha_{n}}x_{1}^{\alpha_{1}}x_{2}^{\alpha_{2}}\cdots
x_{n}^{\alpha_{n}}$

whose zero loci are either non-existent or are contained in the zero locus of
$l_{i}(x)$.

###  3.1. Examples of minimal polynomial folding maps

In the following examples we fix the fc-homotopy class where the minimal
polynomial folding map sits by defining a non-polynomial folding map in the
same class of maps. In particular, a well studied family of folding maps are
the piecewise linear folding maps derived from periodic tessellations of the
$n$-dimensional Euclidean space [HW88]. In this construction one defines an
embedding of the simplex into a Euclidean space of the same dimension, such
that the action of the group of reflections through the facets of the simplex
generates a tessellation of the space (e.g. see Fig. 4). One then uses this
tessellation to determine the associated folding maps. More precisely, we can
construct these piecewise linear folding maps by means of (see [HW88] for more
details on this construction):

  * i)

An affine embedding of $\Delta^{n}$ in $\mathbb{R}^{n}$ such that the embedded
simplex tessellates $\mathbb{R}^{n}$ by means of reflections through the
hyperplanes associated with the $(n-1)$-faces of the simplex. We denote this
affine map by $A\colon\,\Delta^{n}\to\mathbb{R}^{n}$, and we assume that the
image of one of the vertices of $\Delta^{n}$ is the zero vector in
$\mathbb{R}^{n}$.

  * ii)

An integer dilation of the image simplex in $\mathbb{R}^{n}$. In particular,
we require that for any integer dilation, the dilated simplex is tessellated
by a finite number of smaller replicas of itself (e.g. see Fig. 4.).

  * iii)

A map that assigns to each point in every tile of the tessellation of
$\mathbb{R}^{n}$ the equivalent point in the tile that is the image of
$\Delta^{n}$ under $A$. This map is unique; it is constructed by means of
reflections through hyperplanes; and we denote it by
$h\colon\,\mathbb{R}^{n}\to A(\Delta^{n})$.

Now, we define a folding map by first dilating the affine simplex
$A(\Delta^{n})$ $m$ units, i.e. we multiply every vector in
$A(\Delta^{n})\subset\mathbb{R}^{n}$ by the positive integer $m$. We denote
such dilation by $m\circ A(\Delta^{n})$, and by $d$ we denote the number of
smaller replicas of the simplex contained in $m\circ A(\Delta^{n})$.
Therefore, the piecewise linear $d$-folding map is defined as the following
composition of mappings

$A^{-1}\circ h\circ m\circ A\colon\,\Delta^{n}\to\Delta^{n}.$

In what follows we define several piecewise linear folding maps in one and two
dimensions, and then apply the polynomial factorization of Theorem 2 to
determine their associated minimal polynomial folding maps.

####  3.1.1. Folding maps of the one-simplex

The piecewise linear $d$-folding maps in dimension one can be derived from
partitions of the unit interval in $d$ subsegments. In particular, if we
assume that $x=0$ is mapped to $0$ then a partition determines uniquely a
piecewise linear folding map. The simplest non-trivial case, the two-folding
map, can be written as $f_{2,a}(x)=(x/a)\Theta(a-x)+((1-x)/(1-a))\Theta(x-a)$
with $a\in(0,1)$; although for simplicity we consider only $a=1/2$. It follows
from the definition of $f_{2,a}$ that there does not exist a strictly linear
map that is fc-homotopic to this two-folding map. However, if there exists a
polynomial map of degree two that is fc-homotopic to the two-fold, then by
Theorem 2 its corresponding defining polynomials will factor as

(3.3) $P_{1}(x)=Ax(1-x),\quad P_{2}(x)=1-P_{1}(x)=(B+Cx)^{2},$

with $A>0$, $B$ and $C$ being unknowns to be determined. The fc-homotopy
equivalence to $f_{2,a}(x)$ requires that the preimages of $y=0$ are $x=0$ and
$x=1$, and that there exists an interior point $x_{i}=-B/C\in(0,1)$ that is
mapped to $y=1$. We can determine the coefficients $A$, $B$ and $C$ by
expanding the defining equation $P_{2}(x)=1-P_{1}(x)$ in powers of $x$

$(B+Cx)^{2}=1-Ax(1-x),\quad B^{2}-1+(2BC+A)x+(C^{2}-A)x^{2}=0.$

If this is true for every $x\in[0,1]$, then $B=1$, $C=-A/2$, $A^{2}/4=A$, and
$A=4$. Therefore, the corresponding minimal polynomial folding map is defined
by $P_{1}(x)=4x(1-x)$ and $P_{2}(x)=(1-2x)^{2}$, i.e. the logistic map
[May76].

$k$ | $P_{1}(x)$ | $P_{2}(x)=1-P_{1}(x)$
---|---|---
$1$ | $x$ | $1-x$
$2$ | $4x(1-x)$ | $(1-2x)^{2}$
$3$ | $x(3-4x)^{2}$ | $(1-x)(1-4x)^{2}$
$4$ | $16x(1-x)(1-2x)^{2}$ | $\left(8x^{2}-8x+1\right)^{2}$
$5$ | $x(5-20x+16x^{2})^{2}$ | $(1-x)\left(16x^{2}-12x+1\right)^{2}$
$6$ | $4x(1-x)\left((1-4x)(3-4x)\right)^{2}$ | $\left((1-2x)(1-16x+16x^{2})\right)^{2}$
Table 1.  Polynomial factorization of the first six Chebyshev maps from the
$1$-simplex to itself.

One can repeat this construction for every fc-homotopy class of $d$-folding
maps on the one-simplex. As an example, we computed the first six minimal
polynomial folding maps (see Table 1). It is easy to notice that these
polynomials are the normalized Chebyshev polynomials of the first kind
$P_{1,d}(x)=\frac{1-\cos(d\arccos(1-2x))}{2}$. This is indeed expected from
the definition of the Chebyshev polynomials, which employs the piecewise
linear folding function [AR64].

####  3.1.2. Folding maps of the two-simplex

From the theory of affine Weyl groups, we know that every tessellation of the
two-dimensional plane that is derived from the reflections of an isosceles
triangle is isomorphic to one of two possibilities (see Fig. 4.). We denote
each possibility by the name of its corresponding affine Weyl group ($A_{2}$
and $B_{2}$). Of these, only the folding maps derived from the $B_{2}$
tessellation admit polynomial representatives. This is because the folding
maps derived from the $A_{2}$ tessellation map every edge in the boundary of
the triangle to more than one different edge, and therefore the pre-image of
an edge cannot be an algebraic set.

|
---|---
Figure 4.  Two tessellations of the real plane by $2$-simplices. The
tessellation on the left is associated with the Weyl group $A_{2}$ and the
equilateral triangle. The right picture represents a tessellation by isosceles
right triangles, which is associated with the Weyl group $B_{2}$. Both
tessellations give rise to corresponding piecewise linear folding maps on the
$2$-simplex. The dashed green triangles associated with the $g=1$ dilation,
give rise to the $1$-fold; the dashed blue triangles correspond to $g=2$
dilations associated with $4$-folds; the red $g=3$ dilations give rise to
$9$-folds. Of these folds only the ones associated with the isosceles right
triangles admit fc-homotopic polynomial maps.

Given the fc-homotopy class associated with the $B_{2}$ two-fold, using
Theorem 2, and assuming that the minimal polynomial folding map has degree
$k=2$, the defining polynomials have to be of the form

(3.4) $\left(\begin{array}[]{c}P_{1}(x,y)\\\ P_{2}(x,y)\\\
1-P_{1}(x,y)-P_{2}(x,y)\end{array}\right)=\left(\begin{array}[]{c}A(y-Bx)^{2}\\\
(1-x-y)(C+Dx+Ey)\\\ Fxy\end{array}\right).$

Here, $A,B,\ldots,F$, are unknown parameters that satisfy $A>0$, $B>0$, $F>0$
and $C+Dx+Ey>0$ for all $(x,y)\in\Delta^{2}$. Now, we impose Eq. (3.2) on this
family of polynomials, which gives rise to the equation

$Fxy=1-A(y-Bx)^{2}-(1-x-y)(C+Dx+Ey).$

Expanding this equation in the standard basis of polynomials of degree two
with two variables, we get $0=$

$(1-C)+(C-D)x+(C-E)y+(E+D+2AB-F)xy+(D-AB^{2})x^{2}+(E-A)y^{2}.$

Therefore, finding the parameters $(A,B,\dots,F)$ that satisfy this equation
for all $(x,y)\in\Delta^{2}$ amounts to solving the following system of
polynomial equations:

(3.5) $\begin{array}[]{c}1-C=0,\\\ C-D=0,\\\ C-E=0,\\\ E+D+2AB-F=0,\\\
D-AB^{2}=0,\\\ E-A=0.\end{array}$

From the first, second and third equations it follows that $C=D=E=1$, which
yields the strictly positive linear form
$1+x+y\in\mathcal{K}_{+}(\Delta^{n},1)$. From the last equation we find that
$A=1$, and therefore $F$ and $B$ have to obey the fourth and fifth equations
with $A=C=D=E=1$. This implies that $2=F-2B$ and $1-B^{2}=0$. As $B$ has to be
positive in order to yield a folding map compatible with the fc-homotopy
class, $B$ has to be $1$ and $F$ has to be $4$. This shows that there is a
unique polynomial two-fold of degree two in this fc-homotopy class. We can
compactly write this map as

(3.6) $\left(\begin{array}[]{c}P_{1}(x,y)\\\ P_{2}(x,y)\\\
1-P_{1}(x,y)-P_{2}(x,y)\end{array}\right)=\left(\begin{array}[]{c}(y-x)^{2}\\\
(1-x-y)(1+x+y)\\\ 4xy\end{array}\right).$

$f_{d}=$ $\\{P_{1}(x,y)$, $P_{2}(x,y)$,
$P_{3}(x,y)=1-P_{1}(x,y)-P_{2}(x,y)\\}$
---
$f_{1}=$ $\\{x$, $y$, $1-x-y\\}$
$f_{2}=$ $\\{(x-y)^{2}$, $(1-x-y)(1+x+y)$, $4xy\\}$
$f_{4}=$ $\\{(1-2x^{2}-2y^{2})^{2}$, $8xy(1-2xy)$,
$4(1-x-y)(1+x+y)(x-y)^{2}\\}$
$f_{8}=$ $\\{\left(1-4x^{2}+4x^{4}-8xy-4y^{2}+24x^{2}y^{2}+4y^{2}\right)^{2}$,
$8(1-x-y)(1+x+y)(1-2x^{2}+2x^{4}+4xy-2y^{2}-4x^{2}y^{2}+2y^{4})(x-y)^{2}$,
$32xy(1-2xy)(1-2x^{2}-2y^{2})^{2}\\}$
$f_{9}=$
$\\{x(1-y)(2-9x+24x^{2}-16x^{3}+9y-24y^{2}+16y^{3})(3-4x)^{2}(1-4y)^{2}$,
$y(1-x)(2-9y+24y^{2}-16y^{3}+9x-24x^{2}+16x^{3})(3-4y)^{2}(1-4x)^{2}$,
$(1-x-y)^{2}\left(1-8x+16x^{2}-8y-16xy+16y^{2}\right)^{2}\\}$
Table 2.  Polynomial factorization of five different folding maps from the
$2$-simplex to itself.

One can repeat this construction for higher degree folding maps (see Table 2
and Fig. 5). The minimal polynomial folding maps of degree four and eight are
unique and they can be expressed as a composition of the two-folding map; i.e.
$f_{4}=f_{2}\circ f_{2}$ and $f_{8}=f_{2}\circ f_{2}\circ f_{2}$. In the case
of the nine-folding map we only explored the space of maps corresponding to a
single partition of the degrees: $\\{(2,4),\,(2,4),\,(2,4)\\}$. In this
partition the interior pre-images of the boundary correspond to quadratic
curves (see last row in Table 2 and Fig. 5). Although we found a unique
minimal folding map, other different partitions of the degrees in the defining
polynomials could give rise to additional minimal maps. The fact that there
may exist more than one equivalent minimal folding map can be exemplified by
two folding maps of degree eighteen: $f_{9}\circ f_{2}$ and $f_{2}\circ
f_{9}$, which are at the same time distinct and fc-homotopically equivalent
maps.

|
---|---
|
|
|
Figure 5.  Several folding maps on the $2$-simplex. The lines within each
triangle denote the different preimages of the boundary of the simplex under a
given folding map. Each row consists of an fc-homotopy class of folding maps
represented by a piecewise linear map (see left) and a minimal polynomial
folding map (see right and Table 2). The orange dots denote the fixed points
of the minimal polynomial folding maps.

####  3.1.3. Dynamical properties of some polynomial folding maps

The set of Chebyshev maps on the one-simplex has very rich dynamical
properties [AR64]. First, one can show using simple trigonometry that every
Chebyshev map of degree $d$ preserves the measure

$\mu(B)=\frac{1}{\pi}\int_{B}\frac{dx}{\sqrt{x(1-x)}},\quad\,B\subseteq[0,1].$

In addition to this, the corresponding dynamical systems are strongly mixing
with respect to $\mu$; they form a commutative semigroup under the composition
of maps ($f_{d}\circ f_{e}=f_{e}\circ f_{d}=f_{d\times e}$); and they are
orthogonal as polynomial functions under the inner product $L^{2}([0,1],\mu)$
defined by the invariant measure $\mu$ [AR64]. Some of these properties are
also present in deformations of the Chebyshev maps that are not folding maps.
For instance, it is well known that many ergodic properties of the logistic
map $f_{2}$ are preserved by deformations of the type
$f_{2}^{(\lambda)}=\lambda x(1-x)$, where $\lambda$ belongs to a subset of
full measure in the interval $\sim(3.57,4]$ (see [May76]).

Despite these results in dimension one, our knowledge regarding polynomial
folding maps on higher dimensional simplices is non-existent. Because of this,
we performed several numerical experiments with the polynomial folding maps of
degree two and nine that we derived above (Table 2). In the case of the
polynomial two-fold, we determined the fixed points of the first three
iterations of the map (i.e., $f_{2}$, $f_{4}$ and $f_{8}$). We found that the
spectrum of the Jacobian matrix associated with these fixed points lies
outside the unit disk in the complex plane, which is suggestive of chaotic
behavior. We iterated $10^{4}$ times the two-folding map for different random
initial conditions, and found that the dynamics never converged to an
attracting periodic orbit. We repeated this numerical analysis for thousands
of maps in the neighborhood of $f_{2}$ in $\mathsf{PMaps}(\Delta^{2},2)$, and
found that the closest maps to $f_{2}$ exhibited either ergodic behavior or
had an attracting periodic orbit of size larger than $10^{4}$ iterations (see
Fig. 6). In order to perform these experiments, we endowed
$\mathsf{PMaps}(\Delta^{2},2)$ with an $L^{2}$ metric, and computed
deformations of the two-folding map by constructing explicitly a polytopic
approximation of the interior of $\mathsf{PMaps}(\Delta^{2},2)$ (see Appendix
B for details).

---
Figure 6.  Dynamical properties of deformations of the polynomial two-folding
map on the $2$-simplex (see Appendix B). The horizontal axis denotes the
$L^{2}$ distance to the two-folding map. The vertical axis denotes the
smallest absolute value of the eigenvalues of the Jacobian matrix associated
with the fixed points of the map. The green dots denote maps that are either
ergodic or their dynamics converge to a periodic orbit of size larger than
$10^{4}$ iterations. The red dots denote maps whose dynamics converge either
to a finite set of absorbing fixed points or to a periodic orbit of size
smaller than $10^{4}$ iterations.

In the case of the nine-folding polynomial map we found thirteen fixed points.
The Jacobian matrices associated with eleven of these fixed points have their
spectrum located in the exterior of the unit disk, while in the remaining
fixed points, $(1,0)$ and $(0,1)$, the Jacobian matrix has eigenvalues in the
interior of the unit disk. These vertices are attracting fixed points, and
each one is located next to two repelling fixed points (within a Euclidean
distance of less than $1/100$). We computed the dynamics of $10,000$ random
initial conditions close to the origin $(0,0)$ and computed their fixation
times (see Fig. 7). As expected, every scenario converged either to the vertex
$(0,1)$ or $(1,0)$, and the distribution of fixation times closely follows a
log-normal distribution with mean $2.557$ and standard deviation $0.482$.

---
Figure 7.  Distribution of fixation times associated with the dynamics of a
minimal polynomial nine-fold (see Table 2). Here, time is measured as the
number of iterations. The vertex points $(1,0)$ and $(0,1)$ of the triangle
(see Fig. 5) are the only stable fixed points of the map. The distribution
represents the fixation times of $10^{4}$ uniformly distributed random points
in the square $0<x<1/100$, $0<y<1/100$. The shape of the resulting
distribution did not change appreciably when smaller squares were chosen. A
fit of the log-normal distribution ($\hat{\mu}=2.557,\hat{\sigma}=0.482$) to
the simulated data (green curve) is shown superposed.

##  4\. Conclusion

Discrete-generation multi-genic models with non-trivial linkage maps are
difficult to study because of the nonlinearity associated with the dynamical
equations. For instance, simple neutral mutation-recombination models with no
genetic drift require quadratic maps from the simplex to itself to describe
the corresponding population dynamics. More generally, the addition of weak
selection effects combined with time-dependent intensities of selection,
mutation and recombination can give rise to very complex dynamical systems
also described by the iteration of polynomial maps from the simplex to itself
(e.g. see [Lyub92, Gro07, Baa01, NHB99]). Although these maps introduce a
plethora of technical problems that are not present in the case of finite
Markov chains, some particular solutions can be found in biologically
interesting regimes, such as the loose-linkage limit [NHB99] and the tight-
linkage low-mutation limit [Baa01].

In this paper, rather than stressing particular models, we have described some
general properties of the space of stochastic polynomial maps on the simplex.
Our motivation here is twofold. On the one hand we are interested in
developing efficient parametrizations of the space of polynomial maps of
bounded degree on the simplex, with the goal of providing a large class of
diffeomorphisms on the simplex that are useful in practical applications (e.g.
defining adaptive meshes on the simplex). On the other hand, a global
description of the space of polynomial maps allows us to systematically study
properties of the associated dynamical systems. For instance, we have
characterized an important class of maps in the boundary of
$\mathsf{PMaps}(\Delta^{n},k)$, the minimal polynomial folding maps. We have
shown how the dynamics of some of these folding maps and their deformations in
$\mathsf{PMaps}(\Delta^{n},k)$ can exhibit ergodic and mixing properties,
similar to the logistic map in dimension one [May76]. Still, it would be
interesting to explore whether there exist more efficient algorithms to
construct these folding maps. The construction that we have provided in this
paper requires significant computational resources. However, given that the
one-dimensional minimal polynomial folding maps can be derived easily using
classical methods, it remains an open question whether the higher dimensional
folding maps can be derived using simpler methods than the ones that we use
here.

## Acknowledgements

We thank Ethan Akin, Ben Greenbaum, Julien Keller, Stan Leibler, Arnold
Levine, Rami Pugatch, Tiberiu Tesileanu and Tsvi Tlusty for helpful comments.
S.L. is the Addie and Harold Broitman Member in the Simons Center for Systems
Biology at the Institute for Advanced Study, Princeton, NJ.

##  Appendix A Some Elementary Lemmas

######  Lemma 1.

Every one-to-one and onto stochastic matrix in $\mathsf{PMaps}(\Delta^{n},1)$
is an $n+1$ by $n+1$ permutation matrix.

_Proof._ Let $M$ be a one-to-one and onto stochastic matrix, represented by an
$n+1$ by $n+1$ invertible matrix. Let the standard simplex $\Delta^{n}$ be
defined as the convex hull of the standard basis elements in
$\mathbb{R}^{n+1}$, i.e. $e_{1}=(1,0,0,\cdots)$, $e_{2}=(0,1,0,\cdots)$,
$\ldots\,e_{n+1}=(0,0,0,\cdots,1)$. Therefore, the image of $\Delta^{n}$ under
$M$ in $\mathbb{R}^{n+1}$ is the convex hull of the vertices
$\\{Me_{i}\\}_{i=1}^{n+1}$. As $M$ is a stochastic matrix, the image of the
standard simplex is a subset of the standard simplex. Furthermore, as $M$ is
onto, the image simplex is indeed the standard simplex, i.e. the convex hull
of the standard basis elements. This implies that $Me_{i}=e_{j}$ for every
$i\in[1,n+1]$ and $j\in[1,n+1]$. As $M$ is invertible, the image of the
standard basis elements under $M$ spans all the standard basis elements, and
therefore $M$ has to be a permutation matrix.$\square$

######  Lemma 2.

Let $f\colon\,\Delta^{n}\to\Delta^{n}$ be a folding map. The image of the
boundary of the simplex under $f$ is contained in the boundary,
$f(\partial\Delta^{n})\subset\partial\Delta^{n}$.

_Proof._ It follows from the definition of folding map that $f$ is continuous,
onto and $d$-to-$1$, such that the number of preimages is $d$ everywhere in
the interior of $\Delta^{n}$. Let us assume that there exists a point
$x\in\partial\Delta^{n}$ whose image under $f$ is in the interior of the
simplex. Then, one can construct a small open ball $B_{f(x)}$ also in the
interior of the simplex that contains $f(x)$. As $f(x)$ is the image of a
point in the boundary of $\Delta^{n}$ we can divide $B_{f(x)}$ in three
subsets. First, one subset is an open disk of dimension $n-1$ that is the
image under $f$ of a subset in the boundary of $\Delta^{n}$ that contains $x$.
This $n-1$ dimensional disk, that we denote by $D_{f(x)}$, divides the ball in
two open half-balls $B_{f(x)}^{+}$ and $B_{f(x)}^{-}$, such that
$B_{f(x)}=B_{f(x)}^{+}\cup D_{f(x)}\cup B_{f(x)}^{-}$. By the continuity of
$f$, one of the half-balls $B_{f(x)}^{+}$ has a preimage open subset in
$\Delta^{n}$ whose boundary contains $x\in\Delta^{n}\subset\mathbb{R}^{n}$.
However, the other half-ball $B_{f(x)}^{-}$ does not have a preimage in
$\Delta^{n}$ whose boundary contains $x$. It follows that the number of
preimages of $x$ is larger than the number of preimages of points in
$B_{f(x)}^{-}$. This is a contradiction because $f$ is a folding map, and by
definition the number of preimages is the same everywhere in the interior of
the simplex. $\square$

##  Appendix B Deformations of maps

In order to know whether the dynamical properties of a given folding map
$f_{\ast}$ are preserved in its neighborhood in
$\mathsf{PMaps}(\Delta^{n},k)$, we need a systematic method to generate
deformations of $f_{\ast}$ and determine their properties. First, we endow
$\mathsf{PMaps}(\Delta^{n},k)$ with the $L^{2}$ distance associated with the
Euclidean volume form on the simplex, i.e.

$\|f_{\ast}-g\|_{L^{2}}=\sqrt{\int_{\Delta^{n}}\left[\sum_{i=1}^{n}(P_{i}(x)-Q_{i}(x))^{2}\right]\prod_{i=1}^{n}dx_{i}},$

where $f_{\ast},\,g\in\mathsf{PMaps}(\Delta^{n},k)$, and $\\{P_{i}(x)\\}_{i}$
and $\\{Q_{i}(x)\\}_{i}$ are the sets of defining polynomials for $f_{\ast}$
and $g$ respectively. Equipped with this metric we can generate deformations
$g\in B_{f_{\ast},\epsilon}$ of the folding map within a distance $\epsilon$,

$B_{f,\epsilon}=\\{g\in\mathsf{PMaps}(\Delta^{n},k)\colon\,\,\|f_{\ast}-g\|_{L^{2}}\leq\epsilon\\},$

by following the steps outlined in Fig. 8. In particular, we first generate
random interior points $r\in$ $\mathsf{PMaps}(\Delta^{n},k)$ and then take
convex combinations of those points with the folding map to obtain the desired
deformations $g=tr+(1-t)f_{\ast}\in B_{f_{\ast},\epsilon}$. The $t$ parameter
is chosen at random from the uniform distribution on the interval
$(0,t_{\epsilon}]$, with $t_{\epsilon}=\epsilon/\|f_{\ast}-r\|_{L^{2}}$.

---
Figure 8.  Schematic representation of our method to construct deformations of
a folding map. First, we construct random maps in the interior of the space of
polynomial maps (red dots). Then, we use the convexity of
$\mathsf{PMaps}(\Delta^{n},k)$ to construct random convex combinations of
interior points and the folding map such that the resulting deformations are
“close” to the folding map (blue dots).

As we showed in the proof of Theorem 1, $\mathsf{PMaps}(\Delta^{n},k)$ is a
bounded convex subset of the cone
$\times_{i=1}^{n}\mathcal{K}(\Delta^{n},k)\subset\oplus_{i=1}^{n}\Pi^{n}_{k}$.
Hence, in order to parametrize the interior of the space of polynomial maps we
need to first parametrize the convex cone of strictly positive polynomials
$\mathcal{K}_{+}(\Delta^{n},k)$, whose closure is the cone of non-negative
polynomials $\mathcal{K}(\Delta^{n},k)$. This is difficult to attain in
practice, because $\mathcal{K}_{+}(\Delta^{n},k)$ is not finitely generated
for values of $k$ larger than one. However, using Pólya’s criterion one can
construct an infinite sequence of finitely generated cones
$\\{\mathcal{K}_{N}(\Delta^{n},k)\\}_{N=1}^{\infty}$ that satisfy

(B.1)
$\mathcal{K}_{N}(\Delta^{n},k)\subseteq\mathcal{K}_{N+1}(\Delta^{n},k)\subseteq\mathcal{K}_{N+2}(\Delta^{n},k)\subseteq\cdots\to\mathcal{K}_{+}(\Delta^{n},k),$

and converge to $\mathcal{K}_{+}(\Delta^{n},k)$ in the limit of large $N$.
Therefore, in the numerical applications that concern us here, we use the cone
$\mathcal{K}_{N}(\Delta^{n},k)$ with large $N$ as an approximation of
$\mathcal{K}_{+}(\Delta^{n},k)$. This can be done more precisely as follows.
First, we consider the standard parametrization of the vector space of
polynomials $\Pi^{n}_{k}$, in which every polynomial $P(x)$ is expressed as a
linear combination of monomials

(B.2) $P(x)=\sum_{|\alpha|\leq
k}c_{\alpha_{1}\alpha_{2}\cdots\alpha_{n}}x_{1}^{\alpha_{1}}x_{2}^{\alpha_{2}}\cdots
x_{n}^{\alpha_{n}}.$

Then we construct $P_{H}(x)$, the homogenization of $P(x)$, which is obtained
by multiplying every term of degree $|\alpha|<k$ in Eq. (B.2) by the factor
$\left(\sum_{i=1}^{n+1}x_{i}\right)^{k-|\alpha|}$. Finally, we multiply
$P_{H}(x)$ by $(\sum_{i=1}^{n+1}x_{i})^{N}$ and expand the resulting
polynomial

$P_{H}(x)\times\left(\sum_{i=1}^{n+1}x_{i}\right)^{N}$

in the basis of homogenous monomials of degree $N+k$ in $n+1$ variables:

$\left\\{x_{1}^{\beta_{1}}x_{2}^{\beta_{2}}\cdots
x_{n}^{\beta_{n}}x_{n+1}^{\beta_{n+1}}\right\\}_{|\beta|=N+k}.$

The reader can notice that each term in this expansion corresponds to a linear
combination of the $c_{\alpha}$-coefficients that appear in Eq. (B.2). We then
apply Pólya’s theorem, which requires each coefficient of the expansion to be
positive [HLP52]. As the dimension of the space of homogeneous polynomials of
degree $N+k$ in $n+1$ variables is $(N+k+n)!/((N+k)!n!)$, imposing that each
coefficient in the expansion has to be positive gives rise to an equal number
of inequalities in $\Pi^{n}_{k}$. We define $\mathcal{K}_{N}(\Delta^{n},k)$ to
be the semi-algebraic subset that satisfies these inequalities. Furthermore,
in order to parametrize this convex cone one can determine its generators by
means of a vertex enumeration algorithm (e.g. [Avi00]) such that any point in
the cone can be expressed as a linear combination of the generators with
positive coefficients.

In order to learn how one can implement this algorithm, it is useful to
consider a particular example. Here, we consider deformations of the
polynomial two-fold in $\mathsf{PMaps}(\Delta^{2},2)$. In this case, the space
of quadrics has dimension six. To determine the plot shown in Fig. 6, we used
the value $N=8$ to construct the cone $\mathcal{K}_{8}(\Delta^{2},2)$. The
application of Pólya’s criterion gave rise to $66$ inequalities in
$\mathbb{R}^{6}\simeq\Pi_{2}^{2}$, which define a cone with $900$ generators
(we used the reverse search algorithm $\mathsf{lrs}$ to determine the
generators of the cone [Avi00]). We scaled each generator $e_{i}$ by a
positive real number $\lambda_{i}$, such that the maximum value of the
corresponding polynomial on the simplex was one. Given this basis of scaled
generators, we used the Dirichlet distribution with $901$ concentration
parameters

$\alpha_{1}=\alpha_{2}=\cdots=\alpha_{901}=1/1000,$

to generate random convex combinations of the scaled generators. These
combinations correspond to rays in $\mathcal{K}_{8}(\Delta^{2},2)$.
Furthermore, by choosing a random number sampled from the uniform distribution
on $(0,1)$ we generate random quadrics in the segment of the ray defined by
the origin of the cone and the random convex combination of the scaled
generators. We applied this algorithm twice to generate pairs of random
quadrics $\\{R^{\ast}_{1}(x),\,R^{\ast}_{2}(x)\\}$. Then, we computed the
maximum of the function $R^{\ast}_{1}(x)+R^{\ast}_{2}(x)$ on $\Delta^{2}$,
that we denote by $S_{max}$. It follows that sampling a random number
$t^{\ast}$ from the uniform distribution on $[0,1/S_{max}]$, defines an
interior map $r\in\mathsf{PMaps}(\Delta^{2},2)$ with defining polynomials

$R_{1}(x)=t^{\ast}R^{\ast}_{1}(x),\quad R_{2}(x)=t^{\ast}R^{\ast}_{2}(x).$

Finally, the construction of a deformation $g$ of the two-folding map
$f_{\ast}$ requires an additional convex combination of $f_{\ast}$ with the
map $r$ (see Fig. 8.).

## References

  * [AR64] R. L. Adler and T. J. Rivlin, Ergodic and Mixing Properties of Chebyshev Polynomials, Proceedings of the American Mathematical Society, Vol. 15, No. 5 (1964), pp 794-796.
  * [Avi00] D. Avis, lrs: A Revised Implementation of the Reverse Search Vertex Enumeration Algorithm, in: Polytopes - Combinatorics and Computation, G. Kalai & G. Ziegler eds., DMV Seminar Band 29, (2000), pp 177-198.
  * [Baa01] E. Baake, Mutation and recombination with tight linkage, Journal of Mathematical Biology, Vol. 42, Issue 5 (2001), 455-488.
  * [Be22] S.N. Bernstein, Mathematical problems in modern biology, Science in the Ukraine, 1 (1922) pp. 14–19.
  * [Fel88] J. Felsenstein, Phylogenies from molecular sequences: inference and reliability, Annu. Rev. genet. 22, 521-565 (1988).
  * [Gro07] M. Gromov, Mendelian dynamics and Sturtevant’s paradigm, Geometric and Probabilistic Structures in Contemporary Mathematics Series: Dynamics, AMS, Providence RI, Vol. 469 (2007), 227-242.
  * [HLP52] G. H. Hardy, J. E. Littlewood and G. Pólya, Inequalities, 2nd ed., Cambridge University Press, Cambridge, 1952.
  * [HS03] J. Hofbauer and K. Sigmund, Evolutionary game dynamics, Bull. Amer. Math. Soc., Vol. 40 (2003), 479-519.
  * [HW88] M. W. Hoffman and W. D. Withers, Generalized Chebyshev polynomials associated with affine Weyl groups, Transactions of the American Mathematical Society, Vol 308, No. 1 (1988), 91-104.
  * [Lyub92] Y. I. Lyubich, Mathematical structures in population genetics, Biomathematics 22. Springer-Verlag, Berlin, 1992.
  * [LH12] S. Lukić and J. Hey, Demographic inference using spectral methods on SNP data, with an analysis of the human out-of-Africa expansion, Genetics 192 619-639 (2012).
  * [MC05] Gilean A. T. McVean and Niall J. Cardin, Approximating the coalescent with recombination, Phil. Trans. R. Soc. B (2005) 360 1387-1393.
  * [May76] R. M. May, Simple mathematical models with very complicated dynamics, Nature 261, (1976) 459-467.
  * [NHB99] T. Nagylaki, J. Hofbauer and P. Brunovsky, Convergence of multilocus systems under weak epistasis or weak selection, Journal of Mathematical Biology, Vol. 38, Issue 2 (1999), 103-133.
  * [Sal97] L. Saloff-Coste, Lectures on finite Markov chains, Lectures on probability theory and statistics (1997): 301-413.
  * [Sma67] S. Smale, Differentiable dynamical systems, Bulletin of the American Mathematical Society, Vol. 73 (6) (1967), 747–817.
  * [SH13] V. Sousa and J. Hey, Understanding the origin of species with genome-scale data: modeling gene flow, Nat. Rev. Genet. (2013) doi:10.1038/nrg3446.
  * [Ves91] A. P. Veselov, Integrable maps, Russian Math. Surveys, Vol. 46, No. 5 (1991), 1-51.


Positive Rate Binary Interactive Error Correcting Codes Resilient to >1/2 Adversarial Erasures
Meghal Gupta E-mail:
Microsoft Research Rachel Yun ZhangE-mail:. She is supported by an Akamai Presidential Fellowship.

An interactive error correcting code ($\iECC$) is an interactive protocol with the guarantee that the receiver can correctly determine the sender's message, even in the presence of noise. This generalizes the concept of an error correcting code ($\ECC$), which is a non-interactive $\iECC$ that is known to have erasure resilience capped at $\frac12$. The work of [1] constructed the first $\iECC$ resilient to $> \frac12$ adversarial erasures. However, their $\iECC$ has communication complexity quadratic in the message size. In our work, we construct the first positive rate $\iECC$ resilient to $> \frac12$ adversarial erasures. For any $\epsilon > 0$, our $\iECC$ is resilient to $\frac6{11} - \epsilon$ adversarial erasures and has size $O_\epsilon(n)$.

§ INTRODUCTION

Consider the following task: Alice wishes to communicate a message to Bob such that even if a constant fraction of the communicated bits are adversarially tampered with, Bob is still guaranteed to be able to determine her message. This task motivated the prolific study of error correcting codes, starting with the seminal works of [5, 4]. An error correcting code encodes a message $x$ into a longer codeword $\ECC(x)$, such that the Hamming distance between any two distinct codewords is a constant fraction of the length of the codewords.

An important question in the study of error correcting codes is determining the maximal possible error resilience. It is known that in the adversarial bit-flip model, any $\ECC$ can be resilient to at most $\frac14$ corruptions, and in the adversarial erasure error model any $\ECC$ can be resilient to at most $\frac12$ corruptions.

This prompts the following natural question: Can we achieve better error resilience if we use interaction?

In their recent work [1], Gupta, Kalai and Zhang introduce the notion of an interactive error correcting code ($\iECC$), which is an interactive protocol with a fixed length and speaking order, such that Bob can correctly learn Alice's input $x$ as long as not too large a fraction of the total communication is erased.

They demonstrate that $\iECC$'s can in fact achieve a higher erasure resilience than standard error correcting codes. In particular, they design an $\iECC$ that is resilient to adversarial erasure of $\frac35 - \epsilon$ of the total communication.

Note that a classical error correcting code is an $\iECC$ in which Alice speaks in every round. Their result essentially shows that Bob talking occasionally instead of Alice actually improves the error resilience. It is not obvious that this should be the case – since Bob can only send feedback, while Alice can actually send new information, Bob's messages a priori seem a lot less valuable than Alice's. Nevertheless, they are able to leverage this to improve the erasure resilience past $\frac12$.

However, the size of their protocol is quadratic in the length of Alice's original message $x$. This leaves open the question of whether there exists an $\iECC$ achieving $> \frac12$ erasure resilience with size linear in the length of the original message. In this paper, we answer this question to the affirmative.

§.§ Our Results

Our main result is a positive rate $\iECC$ that achieves an erasure resilience of $\frac{6}{11} - \epsilon$ over the binary erasure channel.

For any $\epsilon > 0$, there exists an $\iECC$ over the binary erasure channel resilient to $\frac6{11}-\epsilon$ erasures, such that the communication complexity for inputs of size $n$ is $O_\epsilon(n)$ and the time complexity is $\poly_\epsilon(n)$.

We remark that our $\iECC$ achieves a lower erasure resilience than the quadratic sized $\iECC$ of [1], which is resilient to $\frac35 - \epsilon$ erasures. However, we believe that an $\iECC$ achieving both positive rate and $\frac35 - \epsilon$ erasure resilience can likely be constructed by combining ideas from this paper and [1]. Nevertheless, we leave open the existence of such an $\iECC$.

§.§ Overview of Ideas

In this overview, we briefly review the $\iECC$ of [1] then describe how to modify it to have a linear communication complexity.

The overarching goal of the original protocol, as well as ours, is to perform the following three steps.

* Bob learns that Alice's value of $x$ is one of two possible values. (This idea is known as list decoding, which achieves better noise resilience than unique decoding.)

* Bob conveys to Alice an index $i$ on which the two possible inputs differ.

* Alice sends the value of her input at index $i$.

Summary of the Protocol of [1]. The original protocol consists of many (say $\approx \frac{n}{\epsilon}$) chunks, where in each chunk Alice sends a message followed by Bob's reply. The protocol is designed so that each such chunk will make progress towards Bob's unambiguously learning Alice's input, as long as the adversary did not invest more than $\frac35 - \epsilon$ erasures in that chunk. At a high level, in the first chunk with $< \frac35 - \epsilon$ erasures, Bob narrows down Alice's input to at most two options. In every future chunk with $<\frac35 - \epsilon$ erasures, either Alice gets closer to learning the index $i$ on which the two options differ, or Bob fully determines $x$ by ruling out one of the two values of $x$, e.g. by learning the value of $x[i]$ or by uniquely decoding Alice's message. Alice keeps track of a counter $\cnt$ initially set to $0$ indicating her guess for $i$. The main purpose of Bob's messages is to increment Alice's counter to $i$.

At the beginning of the protocol, Alice sends $\ECC(x, \cnt)$ to Bob in every chunk. At the first point there are $<\frac35 - \epsilon$ erasures in a chunk, Bob will be able to list decode Alice's message to at most two options, say $(x_0, \cnt_0 = 0)$ and $(x_1, \cnt_1 = 0)$. This must happen because the relative message lengths of Alice and Bob will be such that the adversary cannot corrupt too much of Alice's message even if they corrupt none of Bob's message. Since we are in the setting of erasures, one of the two decodings must be Alice's true state, and in particular must contain Alice's true input.

At this point, Bob begins signaling to Alice to increment $\cnt$. His goal is to tell Alice to increment $\cnt$ until $\cnt = i$. He does this by only sending one of two codewords[In our protocol, Bob will send one of four codewords each message. This contributes to the lower erasure resilience of $\frac6{11} - \epsilon$.] every message that have relative distance $1$ apart.
This way, if Bob's message is not entirely erased, Alice learns what Bob tried to send. The key is that every time $< \frac35 - \epsilon$ of a chunk is corrupted,

we can guarantee both that Bob will decode Alice's message to two possible messages, and Alice uniquely decodes Bob's message,[It is also possible that instead Bob uniquely decodes Alice's message, but then he will have uniquely learned $x$.] so that Alice and Bob make progress towards Alice learning $i$.

Once Alice has discovered $i$, Bob signals for Alice to send the bit $x[i]$ for the rest of the protocol,[The reader familiar with the $\iECC$ of [1] may recall that in the case that the two Alices Bob sees have different values of $\cnt$, Bob may instruct Alice to send a different bit for the rest of the protocol, but we do not address this for now.] which allows him to distinguish whether Alice has $x_0$ or $x_1$.

Modifications to Achieve Positive Rate.

The communication complexity of the above protocol is $O(n^2)$. This comes from two parts: (1) $O(n)$ chunks are necessarily for Bob to communicate the index $i \in [n]$ to Alice via incrementation, and (2) Alice sends her length $n$ input in every chunk. We show how to lessen both requirements, thus making the final protocol linear in length.

First, for Bob to communicate $i \in [n]$ to Alice, instead of incrementing $\cnt$ by $1$ until it equals $i$, which requires $O(n)$ rounds of interaction, he builds $i$ bit by bit. That is, Bob writes $i$ out in binary, and then sends Alice each bit of this binary representation in sequence. This only requires $O(\log n)$ rounds of interaction.

Second, we show that instead of sending $x$ every message, it suffices for Alice to encode a shorter string that is different than the corresponding short string for any other $x'$ in most chunks. More precisely, consider an error correcting code $\ECC$ with the following property: set some $\alpha$ and for any $x \not= x'$,

\[
\ECC(x)[j \alpha, (j+1)\alpha - 1] \not= \ECC(x')[j\alpha, (j+1)\alpha - 1]
\]

for all but an $\epsilon$ fraction of values $j$. Then, Alice rotates through the sections, sending $\ECC(x)[j\alpha, (j+1)\alpha - 1]$ in the $\left( j~\text{mod}~\frac{|\ECC(x)|}{\alpha} \right)$'th chunk. Then, if Bob has narrowed down Alice's input to $x_0$ and $x_1$, he can simply ignore the $\epsilon$ fraction of chunks in which $\ECC(x_0)[j\alpha, (j+1)\alpha - 1] = \ECC(x_1)[j\alpha, (j+1)\alpha - 1]$. In the remainder of chunks, the segment $\ECC(x)[j\alpha, (j+1)\alpha - 1]$ is sufficient for Bob to distinguish between $x_0$ and $x_1$. If we were to let $\alpha\approx \log{n}$,[$\alpha = \Theta(\log n)$ is necessary since Alice also sends her current guess of $i$ each message, which has length $\log n$.] then our chunks are now only length $O(\log{n})$.

Combining the two modifications, we see that $\Theta(n)$ communication from Alice is necessary for Bob to narrow down Alice's input to two options, and then after that, Bob can convey $i$ to Alie in $O(\log n)$ chunks each of size $O(\log n)$. This results in an $\iECC$ with total communication $O(n + \log^2 n) = O(n)$.

We remark that our protocol has erasure resilience $\frac6{11} - \epsilon$. The limiting factor is in the construction of a protocol in which Bob builds $i$ bit by bit: our protocol requires Bob sending $4$ codewords with distance $\frac23$, which could possibly be improved to $2$ codewords with distance $\frac12$ to achieve $\frac35 - \epsilon$ erasure resilience, though we do not do that here. However, combining our second observation with the protocol from [1] would be enough to give an $\iECC$ with $\frac35 - \epsilon$ erasure resilience with communication $O(n \log n)$.

§ PRELIMINARIES AND DEFINITIONS

Before we dive into the technical part of our paper, we present important preliminaries on classical error correcting codes, and define an $\iECC$ formally and what it means for one to be resilient to $\alpha$-fraction of erasures.

Notation. In this work, we use the following notations.

* The function $\Delta(x, y)$ represents the Hamming distance between $x$ and $y$.

* The interval $[a,b]$ for $a,b\in \bbZ_{\ge 0}$ denotes the integers from $a$ to $b$ inclusive. The interval $[n]$ denotes the integers $1, \dots, n$.

* The symbol $\perp$ in a message represents the erasure symbol that a party might receive in the erasure model.

* When we say Bob $k$-decodes a message, we mean that he list decodes it to exactly $k$ possible messages Alice could have sent in the valid message space.

* The output of an $\ECC$ is $0$-indexed. All other strings are $1$-indexed.

§.§ Classical Error Correcting Codes

An error correcting code ($\ECC$) is a family of maps $\ECC = \{ \ECC_n : \{ 0, 1 \}^n \rightarrow \{ 0, 1 \}^{m(n)} \}_{n \in \bbN}$. An $\ECC$ has relative distance $\alpha > 0$ if for all $n \in \bbN$ and any $x \not= y \in \{ 0, 1 \}^n$,

\[
\Delta \left( \ECC_n(x), \ECC_n(y) \right) \ge \alpha m(n).
\]

Binary error correcting codes with relative distance $\approx \frac12$ are well known to exist with linear blowup in communication complexity.

For all $\epsilon > 0$, there exists an explicit linear error correcting code $\ECC_\epsilon = \{ \ECC_{\epsilon, n} : \{ 0, 1 \}^n \rightarrow \{ 0, 1 \}^{m} \}_{n \in \bbN}$ with relative distance $\frac12 - \epsilon$ and with $m = m(n) = O_\epsilon (n)$. Furthermore, all codewords other than $\ECC_{\epsilon, n}(0^n)$ are relative distance $\frac12 - \epsilon$ from $0^m$ and $1^m$ as well.

A relative distance of $\frac12$ is in fact optimal in the sense that as the number of codewords $N$ approaches $\infty$, the maximal possible relative distance between $N$ codewords approaches $\frac12$. We remark, however, that for small values of $N$, the distance can be much larger: for $N = 2$, the relative distance between codewords can be as large as $1$, e.g. the codewords $0^M$ and $1^M$, and for $N = 4$, the relative distance can be as large as $\frac23$, e.g. the codewords $(000)^M, (110)^M, (101)^M, (011)^M$. Our constructions leverage this fact that codes with higher relative distance exist for a small constant number of codewords.

We will also need the following important lemma about the number of shared bits between any three codewords in an error correcting code scheme that has distance $\frac12$.

For any error correcting code $\ECC_\epsilon = \{ \ECC_{\epsilon, n} : \{ 0, 1 \}^n \rightarrow \{ 0, 1 \}^{m} \}_{n \in \bbN}$ with relative distance $\frac12 - \epsilon$, and any large enough $n \in \bbN$, any three codewords in $\ECC_{\epsilon, n}$ overlap on at most $\left( \frac14 + \frac32 \epsilon \right) \cdot m$ locations.

Consider three codewords $c_1, c_2, c_3$. Suppose that all pairs are relative distance at least $\left( \frac12 - \epsilon \right)$. Let $c_1$ and $c_2$ share $f \le \left( \frac12 + \epsilon \right) \cdot m$ bits, and all three codewords share $e \le f$ bits. Then, note that

\begin{align*}
2m \cdot \left( \frac12 - \epsilon \right)
&\le \Delta(c_1, c_3) + \Delta(c_3, c_2) \\
&\le 2(f-e) + (m-f) \\
&= m + f - 2e \\
&\le \left( \frac32 + \epsilon \right) \cdot m - 2e,
\end{align*}

which means that

\[
e \le \left( \frac14 + \frac32 \epsilon \right) \cdot m,
\]

as claimed.

Lemma <ref> means that assuming that $< \frac34$ of a codeword is erased, the resulting message is list-decodable to a set of size $\le 2$, at least in theory. The following theorem says such a code exists with list-decoding being polynomial time, while also satisfying a couple other properties necessary in the protocol construction in Section <ref>.

For all $\epsilon>0$, any explicit (given with its encoding matrix) linear code $\ECC_\epsilon = \{ \ECC_{\epsilon, n} : \{ 0, 1 \}^n \rightarrow \{ 0, 1 \}^{m} \}_{n \in \bbN}$ with relative distance $(\frac12-\epsilon)$, can be efficiently decoded and list-decoded. That is, there exists a $\poly_\epsilon(n)$-time decoding algorithm $\DEC_\epsilon = \{ \DEC_{\epsilon,n} : \{ 0, 1 \}^m \rightarrow \mathcal{P}(\{ 0, 1 \}^n) \}_{n \in \bbN}$, such that for any $n \in \bbN$, $x \in \{ 0, 1 \}^n$, and corruption $\sigma$ consisting of fewer than $(\frac12 - \epsilon) \cdot m$ erasures,

\[
x = \DEC_{\epsilon,n}(\sigma \circ \ECC_{\epsilon, n}(x)).
\]

Moreover, for any corruption $\sigma$ consisting of fewer than $(\frac34 - \frac32\epsilon) \cdot m$ erasures,

\[
\left| \DEC_{\epsilon,n}(\sigma \circ \ECC_{\epsilon, n}(x)) \right| \le 2, \qquad x \in \DEC_{\epsilon,n}(\sigma \circ \ECC_{\epsilon, n}(x)).
\]

Our following theorem gives an $\ECC$ such that any two codewords differ on most segments of length $\alpha$.

For all $n \in \bbN, \epsilon > 0$, there exists $\alpha = \Theta_\epsilon(\log{n})$ such that, there exists an explicit linear code $\ECC_\epsilon = \{ \ECC_{\epsilon, n} : \{ 0, 1 \}^n \rightarrow \{ 0, 1 \}^{m} \}_{n \in \bbN}$ with $m = m(n) = O_\epsilon (n)$, satisfying the following property:
For all $n \in \bbN, \epsilon > 0$, it holds that $\alpha | m$ and for any $x \not= x' \in \{ 0, 1 \}^n$ and $j\in \{0\dots \frac{m}{\alpha}-1\}$,

\[
\ECC_{\epsilon, n}(x)[ j\alpha : (j+1)\alpha - 1 ] = \ECC_{\epsilon, n}(x')[ j\alpha : (j+1)\alpha - 1 ] \label{prop:jlogn}
\]

for at most $\frac{\epsilon m}{\alpha}$ values of $j$.

The paper of [2] and many others provide an explicit concatenated linear code where the outer code is a large alphabet ($|\Sigma|=O_\epsilon(1)$) linear code $\ECC'_{\epsilon, n} : \{ 0, 1 \}^n \rightarrow \Sigma^{m'}$ of distance $1-\epsilon$ and the inner code is a distance $\frac12$ Hadamard code $\Had : \Sigma \rightarrow \{ 0, 1 \}^r$, so that $m = m'r$.

We choose $\alpha'=\lfloor\frac{\epsilon}{2}\log{n}\rfloor$ and assume that $\alpha' | m'$ (this can be done by possibly padding the outer code with up to $\alpha'$ $0$'s; note that the code is still linear after padding with $0$'s).

Now, let $\alpha=\alpha'r$. Clearly, $\alpha | m$. To see why (<ref>) holds, notice that in order for some $j \in [0, \frac{m}{\alpha} - 1]$ to satisfy

\[
\ECC_{\epsilon, n}(x)[j\alpha:(j+1)\alpha-1] = \ECC_{\epsilon, n}(x')[j\alpha:(j+1)\alpha-1]
\]

it must be the case that

\[
\ECC'_{\epsilon, n}(x) \left[ j\alpha', (j+1)\alpha' - 1 \right]
= \ECC'_{\epsilon, n}(x') \left[ j\alpha', (j+1)\alpha' - 1 \right].
\]

Let $J$ be the number of such $j \in [0, \frac{m}{\alpha} - 1]$. Then since $\ECC'_{\epsilon, n}$ has distance $\ge 1 - \epsilon$ it holds that $J \le \frac{\epsilon m}{\alpha}$, as claimed.

Let $\ECC'_\epsilon = \{ \ECC'_{n, \epsilon} : \{ 0, 1 \}^n \rightarrow \{ 0, 1 \}^{m(n)} \}$ be a explicit linear code satisfying the properties of Theorem <ref> with $\alpha = \alpha_n = \Theta_\epsilon(\log n)$. For all linear $\ECC_\epsilon = \{ \ECC_{n, \epsilon} : \{ 0, 1 \}^{\alpha} \times \{ 0, 1 \}^{\beta} \rightarrow \{ 0, 1 \}^{p(n)} \}$ with relative distance $\frac12-\epsilon$ and for all $\beta$, the code defined by

\[
C(x) = \ECC_\epsilon(\ECC'_\epsilon(x)[0,\alpha-1], 0^\beta) || \dots || \ECC_\epsilon(\ECC'_\epsilon(x)[m-\alpha, m-1], 0^\beta)
\]

is a linear code with relative distance $\frac12-\frac32\epsilon$. In particular, assuming that less than $\frac34 - \frac94 \epsilon$ of $C(x)$ is erased, there is an efficient algorithm to obtain a set of size $2$ containing $x$.

Regardless of the choice of $\ECC$, since both $\ECC$ and $\ECC'$ are linear, it follows that $C$ is linear as well.

We now show that the relative distance between $C(x)$ and $C(x')$ is at least $\frac12-\frac32\epsilon$. By Lemma <ref>, at most $\epsilon$ chunks $\ECC'[j\alpha,(j+1)\alpha-1]$ are identical for $x$ and $x'$. For the remaining $1-\epsilon$ fraction of the chunks, the relative distance is at least $\frac12-\epsilon$, so the total relative distance is at least

§.§ Interactive Error Correcting Codes

We formally define our notion of an interactive error correcting code ($\iECC$). The two types of corruptions we will be interested in are erasures and bit flips. We first start by defining a non-adaptive interactive protocol.

A non-adaptive interactive protocol $\pi = \{ \pi_n \}_{n \in \bbN}$ is an interactive protocol between Alice and Bob, where in each round a single party sends a single bit to the other party. The order of speaking, as well as the number of rounds in the protocol, is fixed beforehand. The number of rounds is denoted $|\pi|$.

An interactive error correcting code ($\iECC$) is a non-adaptive interactive protocol $\pi = \{ \pi_n \}_{n \in \bbN}$, with the following syntax:

* At the beginning of the protocol, Alice receives as private input some $x \in \{ 0, 1 \}^n$.

* At the end of the protocol, Bob outputs some $\hat{x} \in \{ 0, 1 \}^n$.

We say that $\pi$ is $\alpha$-resilient to adversarial bit flips (resp. erasures) if there exists $n_0 \in \bbN$ such that for all $n > n_0$ and $x \in \{ 0, 1 \}^n$, and for all online adversarial attacks consisting of flipping (resp. erasing) at most $\alpha \cdot |\pi|$ of the total communication, Bob outputs $x$ at the end of the protocol with probability $1$.

There exists a family of error correcting codes $\ECC = \{ \ECC_n : \{ 0, 1 \}^n \rightarrow \{ 0, 1 \}^{M(n)} \}_{n \in \bbN}$ with relative distance $\frac12 - \epsilon$ and $M = O_\epsilon(n)$ that satisfies that for all $x, y \in \{ 0, 1 \}^n$,

\[
\left| \left\{ j \in [M/\log n] : \ECC(x) \left[ \frac{(j-1)M}{\log n} + 1 : \frac{jM}{\log n} \right] = \ECC(y) \left[ \frac{(j-1)M}{\log n} + 1 : \frac{jM}{\log n} \right] \right\} \right| \le \frac{\epsilon M}{\log n}.
\]

Note that this holds for a random code: for any two strings $x, y \in \{ 0, 1 \}^n$ the expected number of the $M/\log n$ strings that are equal is $2^{-\log n} \cdot M/\log n$ and the variance is $\approx 2^{-\log n} \cdot M/\log n$, so the probability that at least $\epsilon$ fraction of the $M/\log n$ strings are the same is at most $\approx 2e^{-\epsilon^2 n \log n \cdot M}$ by Chernoff. Union bounding over all $2^{2n}$ pairs of strings $x, y$, we get that the probability the above property holds is at least $\approx 1$.

Let $\ECC$ be a family of error correcting codes $\{ \ECC_n : \{ 0, 1 \}^{\log n} \times \{ 0, 1 \}^{\le \log n} \rightarrow \{ 0, 1 \}^{P(n)}$ with relative distance $\frac12 - \epsilon$ and $P = O_\epsilon(\log n)$.

Define $C : \{ 0, 1 \}^n \times (\{ 0, 1 \}^{\le \log n})^{M/\log n}$

\[
C(x; r_1, \dots, r_{M/\log n}) = \ECC \left( \ct\left[ 1 : \log n \right], r_1 \right) || \dots || \ECC \left( \ct\left[ M - \log n + 1 : M \right], r_{M/\log n} \right),
\]

where $\ECC$ is a family of error correcting codes $\{ \ECC_n : \{ 0, 1 \}^{\log n} \times \{ 0, 1 \}^{\le \log n} \rightarrow \{ 0, 1 \}^{P(n)}$ with relative distance $\frac12 - \epsilon$, and $\ct = \ECC'(x)$ where $\ECC' = \{ \ECC'_n : \{ 0, 1 \}^n \rightarrow \{ 0, 1 \}^N \}_{n \in \bbN}$ satisfies the properties listed in Lemmma <ref>.

Then for any $x \not= y \in \{ 0, 1 \}^n$ and any $(r_1, \dots, r_{M/\log n}), (s_1, \dots, s_{M/\log n}) \in (\{ 0, 1 \}^{\le n})^{M/\log n}$, it holds that

\[
\Delta(C(x; r_1, \dots, r_{M/\log n}), C(y; s_1, \dots, s_{M/\log n})) \ge \left( \frac12 - \frac32 \epsilon \right) \cdot \frac{MP}{\log n}.
\]

On all but at most $\frac{\epsilon M}{\log n}$ of the intervals $\left[ \frac{(j-1)M}{\log n} + 1, \frac{jM}{\log n} \right]$, $\ECC'(x)$ and $\ECC'(y)$ are not identical. For each such interval $\left[ \frac{(j-1)M}{\log n} + 1, \frac{jM}{\log n} \right]$,

\[
\Delta \left(
\ECC \left( \ECC'(x) \left[ \frac{(j-1)M}{\log n} + 1, \frac{jM}{\log n} \right], r_j \right),
\ECC \left( \ECC'(y) \left[ \frac{(j-1)M}{\log n} + 1, \frac{jM}{\log n} \right], s_j \right)
\right)
\ge \left( \frac12 - \epsilon \right) \cdot P,
\]

so that

\begin{align*}
\Delta(C(x; r_1, \dots, r_{M/\log n}), C(y; s_1, \dots, s_{M/\log n}))
\ge \frac{(1-\epsilon)M}{\log n} \cdot \left( \frac12 - \epsilon \right) \cdot P
\ge \left( \frac12 - \frac32 \epsilon \right) \cdot \frac{MP}{\log n}.
\end{align*}

[1]
Meghal Gupta, Yael Tauman Kalai, and Rachel Yun Zhang.
Interactive error correcting codes over binary erasure channels
resilient to $>\frac12$ adversarial corruption, 2021.

[2]
Venkatesan Guruswami and Madhu Sudan.
List decoding algorithms for certain concatenated codes.
In Proceedings of the Thirty-Second Annual ACM Symposium on
Theory of Computing, STOC '00, page 181–190, New York, NY, USA, 2000.
Association for Computing Machinery.

[3]
V. Guruswami.
List decoding from erasures: bounds and code constructions.
IEEE Transactions on Information Theory, 49(11):2826–2833,

[4]
R. W. Hamming.
Error detecting and error correcting codes.
The Bell System Technical Journal, 29(2):147–160, 1950.

[5]
C. E. Shannon.
A mathematical theory of communication.
The Bell System Technical Journal, 27(3):379–423, 1948.


* [GO18b] Lyudmila Grigoryeva and Juan-Pablo Ortega.  Universal discrete-time reservoir computers with stochastic inputs and linear readouts using non-homogeneous state-affine systems.  Journal of Machine Learning Research, 19(24):1–40, 2018.
  * [GO20] Lukas Gonon and Juan-Pablo Ortega.  Reservoir computing universality with stochastic inputs.  IEEE Trans. Neural Netw. Learn. Syst., 31(1):100–112, 2020.
  * [GS20] Claudio Gallicchio and Simone Scardapane.  Deep randomized neural networks.  In Recent Trends in Learning From Data: Tutorials from the INNS Big Data and Deep Learning Conference (INNSBDDL2019), pages 43–68. Springer, 2020.
  * [GSS21] Friedrich Götze, Holger Sambale, and Arthur Sinulis.  Concentration inequalities for polynomials in $\alpha$-sub-exponential random variables.  Electron. J. Probab., 26:Paper No. 48, 22, 2021.
  * [Ham20] Hans Hamburger.  Über eine Erweiterung des Stieltjesschen Momentenproblems.  Mathematische Annalen, 81(2-4):235–319, 1920.
  * [HD20] Moritz Helias and David Dahmen.  Statistical field theory for neural networks, volume 970 of Lecture Notes in Physics.  Springer, Cham, [2020] ©2020.
  * [HP91] Peter Hilton and Jean Pedersen.  Catalan numbers, their generalization, and their uses.  The mathematical intelligencer, 13:64–75, 1991.
  * [HS97] Sepp Hochreiter and Jürgen Schmidhuber.  Long short-term memory.  Neural computation, 9(8):1735–1780, 1997.
  * [Hwa04] Suk-Geun Hwang.  Cauchy’s interlace theorem for eigenvalues of hermitian matrices.  The American mathematical monthly, 111(2):157–159, 2004.
  * [JGP+19] Li Jing, Caglar Gulcehre, John Peurifoy, Yichen Shen, Max Tegmark, Marin Soljacic, and Yoshua Bengio.  Gated orthogonal recurrent units: On learning to forget.  Neural computation, 31(4):765–783, 2019.
  * [JH04] Herbert Jaeger and Harald Haas.  Harnessing nonlinearity: Predicting chaotic systems and saving energy in wireless communication.  Science, 304(5667):78–80, 2004.
  * [JL84] William B. Johnson and Joram Lindenstrauss.  Extensions of Lipschitz mappings into a Hilbert space.  In Conference in modern analysis and probability (New Haven, Conn., 1982), volume 26 of Contemp. Math., pages 189–206. Amer. Math. Soc., Providence, RI, 1984.
  * [Lat06] Rafał Latała.  Estimates of moments and tails of Gaussian chaoses.  The Annals of Probability, 34(6):2315 – 2331, 2006.
  * [LHEL21] Zhong Li, Jiequn Han, Weinan E, and Qianxiao Li.  On the curse of memory in recurrent neural networks: Approximation and optimization analysis.  In International Conference on Learning Representations, 2021.
  * [LKDB18] Floris Laporte, Andrew Katumba, Joni Dambre, and Peter Bienstman.  Numerical demonstration of neuromorphic computing with photonic crystal cavities.  Opt. Express, 26(7):7955–7964, Apr 2018.
  * [MJS07] Wolfgang Maass, Prashant Joshi, and Eduardo D Sontag.  Computational aspects of feedback in neural circuits.  PLoS Comput Biol, 3(1):e165, 2007.
  * [SCS88] H. Sompolinsky, A. Crisanti, and H.-J. Sommers.  Chaos in random neural networks.  Phys. Rev. Lett., 61(3):259–262, 1988.
  * [SH06] M.D. Skowronski and J.G. Harris.  Minimum mean squared error time series classification using an echo state network prediction model.  In 2006 IEEE International Symposium on Circuits and Systems, pages 4 pp.–3156, 2006.
  * [ST50] James Alexander Shohat and Jacob David Tamarkin.  The problem of moments, volume 1.  American Mathematical Society (RI), 1950.
  * [Sze36] Gabriel Szegö.  On some Hermitian forms associated with two given curves of the complex plane.  Trans. Amer. Math. Soc., 40(3):450–461, 1936.
  * [Tin20] Peter Tino.  Dynamical systems as temporal feature spaces.  J. Mach. Learn. Res., 21:Paper No. 44, 42, 2020.
  * [TV08] Terence Tao and Van Vu.  Random matrices: the circular law.  Communications in Contemporary Mathematics, 10(02):261–307, 2008\.
  * [VALT22] Pietro Verzelli, Cesare Alippi, Lorenzo Livi, and Peter Tiňo.  Input-to-state representation in linear reservoirs dynamics.  IEEE Trans. Neural Netw. Learn. Syst., 33(9):4598–4609, 2022.
  * [Ver18] Roman Vershynin.  High-dimensional probability, volume 47 of Cambridge Series in Statistical and Probabilistic Mathematics.  Cambridge University Press, Cambridge, 2018.  An introduction with applications in data science, With a foreword by Sara van de Geer.
  * [VMVV+14] Kristof Vandoorne, Pauline Mechet, Thomas Van Vaerenbergh, Martin Fiers, Geert Morthier, David Verstraeten, Benjamin Schrauwen, Joni Dambre, and Peter Bienstman.  Experimental demonstration of reservoir computing on a silicon photonics chip.  Nature communications, 5(1):3541, 2014.
  * [VSP+17] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Łukasz Kaiser, and Illia Polosukhin.  Attention is all you need.  Advances in neural information processing systems, 30, 2017.
  * [Wig55] Eugene P. Wigner.  Characteristic vectors of bordered matrices with infinite dimensions.  Annals of Mathematics, 62(3):548–564, 1955.
  * [Yap00] Chee Keng Yap.  Fundamental problems of algorithmic algebra.  Oxford University Press, New York, 2000.