text:
\begin{document}

\markboth{Berry and Sauer}{Journal of Difference Equations and Applications}

\title{{\itshape Convergence of Periodically-Forced Rank-Type Equations}}

\author{Tyrus Berry\thanks{Corresponding author. Email: tberry@gmu.edu} and Timothy Sauer \\\vspace{6pt}  {\em{George Mason University; Fairfax, VA  22030}; }}

\newcommand{\ov}[0]{\overline}

\maketitle

\begin{abstract} Consider a difference equation which takes the k-th largest output of $m$ functions  of the previous $m$ terms of the sequence.  If the functions are also allowed to change periodically as the difference equation evolves this is analogous to a differential equation with periodic forcing.  A large class of such non-autonomous difference equations are shown to converge to a periodic limit which is independent of the initial condition.  The period of the limit does not depend on how far back each term is allowed to look back in the sequence, and is in fact equal to the period of the forcing.\end{abstract}

{\em To appear in the Journal of Difference Equations and Applications}

\section{Introduction}
Recently, there has been substantial interest in max-type difference equations \cite{GL,Voulov1,Voulov2,Voulov3,BFS,KR,Yang,BGKL,Chen}.
For the  equation
\begin{equation} x_n = \max_{1\leq i \leq M} \{f_i(x_{n-i}) \} \end{equation}
with initial sequence $(x_1,...,x_M)$, it has been shown  \cite{sauerMax} that if the $f_i$ are contractive, then all solutions converge to a fixed point, i.e.
\begin{equation*} \lim_{n\to\infty} x_n = r_*. \end{equation*}
Moreover, the fixed point $r_*$ can be identified as the maximum of the individual fixed points
\begin{equation*} r_* = \max_{1\le i\le M}\{r_i\}, \end{equation*}
where $r_i$ is the unique fixed point $r_i=f_i(r_i)$.

In view of this result, it is reasonable to ask about convergence in the periodically-forced case. Assume  that  $f_1, \ldots, f_M$ vary periodically with the discrete time variable $n$. We will investigate limiting behavior of the nonautonomous difference equation
\begin{equation}\label{maxeq} x_n = \max_{1\leq i \leq M} \{f_i(x_{n-i},n) \} \end{equation}
where $f_i:\mathbb{R}\times \mathbb{N} \to \mathbb{R}$ are contractions with respect to the first variable and $P$-periodic with respect to the second.  That is, we assume there exists $\alpha <1$ such that
\begin{equation} \label{alpha}  |f_i(x,n)-f_i(y,n)| \leq \alpha|x-y|
\end{equation}
for all $i,n$, and
$ f_i(x,n+P) = f_i(x,n)$ for all $n$.
We think of $P$ as the {\em forcing period} and $M$ as the {\em memory length}.  We will show below that solutions of the contractive, $P$-periodic difference equation (\ref{maxeq}) converge to a unique periodic orbit with period $P$, for all initial sequences.

In \cite{sauerRank}, rank-type equations were proposed as a generalization of max-type equations. Consider the difference equation
\begin{equation}\label{rankeq} x_n = \underset{1\leq i \leq M}{\textup{$k$-rank}} \{f_i(x_{n-i}) \} \end{equation}
where $k$-rank denotes the $k$th largest value. It was shown in \cite{sauerRank} that with the same contractiveness hypotheses, solutions converge to $r_* = \text{ $k$-rank\  } \{ r_i\}$, the $k$th largest of the individual fixed points. In this paper, we further generalize this result to the periodically-forced case. Our main result is the following.

\begin{theorem} \label{theorem1}
Let $f_1,\ldots,f_M:\mathbb{R}\times \mathbb{N} \to \mathbb{R}$ be contractions with respect to the first variable and $P$-periodic with respect to the second. Let $1\leq k \leq M$. Then for any initial sequence, the solution of the difference equation
\begin{equation}\label{rankeq1} x_n =  \underset{1\leq i \leq M}{\textup{$k$-rank}} \{f_i(x_{n-i},n) \} \end{equation}
is asymptotically periodic with period $P$.
\end{theorem}

Unlike the autonomous case, in general we do not know how to find a formula for the periodic solution in terms of the individual dynamics of the $f_i$. In the fourth section of this article we relate some partial progress in this direction.

The conclusion of Theorem \ref{theorem1} fails if $\alpha\geq 1$ in (\ref{alpha}), even if $k=P=1$, the autonomous max-type case \cite{GL}. We mention two well-known examples where the condition on $\alpha$ does not hold.
\begin{example}
Consider the difference equation
\[x_n=\max\{-x_{n-1},-x_{n-2}\}.\]
Although this is an autonomous max-type equation ($k=P=1$), is easy to check that all solutions have period $3\neq P$. In this example, $\alpha=1$, and Theorem \ref{theorem1} does not apply.
\end{example}
\begin{example}
The first-order difference equation
\[x_n=\max\{1-2x_{n-1},2x_{n-1}-1\}\]
has bounded, nonperiodic (chaotic) solutions for almost every initial condition $x_0$. In this example, $\alpha=2$. As a dynamical system, this example is the upside-down tent map.
\end{example}
\begin{example}\label{examp0}
For a straightforward application of Theorem \ref{theorem1}, fix positive integers $P$ and $k\leq M$, and denote by $A$ and $B$ two $P\times M$ matrices of real numbers, where the entries $|A_{ij}|<1$. Define the difference equation
\begin{equation}
x_n = \textup{$k$-rank} \{A_{\overline{n}1}x_{n-1}+B_{\overline{n}1}, A_{\overline{n}2}x_{n-2}+B_{\overline{n}2}, \ldots, A_{\overline{n}M}x_{n-M}+B_{\overline{n}M}\}
\end{equation}
where $\overline{n}\equiv  1+(n-1 \mod P)$. Thus the coefficients cycle through the rows of the matrices $A$ and $B$, forcing the equation with period $P$. Theorem \ref{theorem1} implies that all solutions converge asymptotically to a period $P$ solution.
\end{example}

\begin{example}\label{ex1} Let $A$ be a $P\times M$ matrix of positive real numbers and consider the recurrence
\begin{equation} \label{eq31}
 x_n = \max\{A_{\ov n1}x_{n-1}^{\alpha_1},\ldots,A_{\ov nM}x_{n-M}^{\alpha_m} \}
 \end{equation}
where $-1 < \alpha_i < 1$ and $\overline{n}= 1+(n-1 \mod P)$.
It follows from Theorem \ref{theorem1} (applied to $y_n=\ln x_n$) that the recurrence converges to a unique period $P$ orbit.   Note that this result is independent of the memory length $M$. The periodicity of the limit only depends on the periodicity of the forcing.

In sufficiently simple cases, we can say more about the convergent solution.  In Section 4, we further pursue the special case $M=P=2$ of (\ref{eq31}). A closed form for the globally attracting solution of (\ref{eq31})  is given by
\begin{eqnarray*} \lim_{n\to\infty} x_{2n} &=& \max \left\{ A_{12}A_{21}^{\frac{\alpha_1}{1-\alpha_2}},A_{11}^{\frac{\alpha_1}{1-\alpha_1^2}}A_{12}^{\frac{1}{1-\alpha_1^2}}, A_{22}^{\frac{1}{1-\alpha_2}} \right\} \\
\lim_{n\to\infty} x_{2n+1} &=& \max \left\{ A_{11}A_{22}^{\frac{\alpha_1}{1-\alpha_2}},A_{12}^{\frac{\alpha_1}{1-\alpha_1^2}}A_{11}^{\frac{1}{1-\alpha_1^2}}, A_{21}^{\frac{1}{1-\alpha_2}} \right\}.
 \end{eqnarray*}
 (See Example \ref{ex44} for details.)
The complexity of this solution contrasts with the simplicity of the case $P=1$, which is simply $x_n \to \max\left\{A_i^{\frac{1}{1-\alpha_i}}\right\}$ (see \cite{sauerMax,Sun,Stev}).
\end{example}

\begin{example} Let $A_1, A_2, A_3$ be real numbers less than $0.15$, $B_1,B_2, B_3$ be arbitrary real numbers, and $P$ be a positive integer. Then it follows from Theorem \ref{theorem1} that every solution of
\[ x_n = \textup{median} \left\{ e^{A_1\sin(B_1+2\pi n/P) - x_{n-1}^2}, e^{A_2\sin(B_2+2\pi n/P) - x_{n-2}^2},
e^{A_3\sin(B_3+2\pi n/P) - x_{n-3}^2}\right\} \]
is asymptotically periodic with period $P$.  Note that median is synonymous with 2-rank. It is easily checked that the condition $A_i<0.15< \frac{1}{2} - \frac{\ln 2}{2}$ implies that the functions $f_i(x)=\exp(A_i\sin(B_i+2\pi n/P)-x^2)$ are contractive, so the convergence is implied by Theorem \ref{theorem1}.
\end{example}

Our main convergence result Theorem \ref{theorem1} will be proved as a special case of a more general result, Corollary \ref{mainresult}, which applies to a class of difference equations called sup-contractive. The next example is covered under the sup-contractive hypothesis, which is more general than a fixed $k$-rank.

\begin{example}\label{ex3} Let $f_1,...,f_M : \mathbb{R} \to \mathbb{R}$ be contractions.  Let $P \leq M$ be a positive integer and denote  $\ov{n}=1+(n\mod P)$.  We will show in Example \ref{ex37} that every solution of
\[ x_n = \frac{1}{2} \left( \max_{1\leq i\leq M}\{f_i(x_{n-i})\} - \underset{1\leq i\leq M}{\overline{n}\textup{-rank}}
\{f_i(x_{n-i})\} \right) \]
is asymptotically periodic with period $P$.
\end{example}

In Section 2 below, we develop some facts about functions that are contractive in the sup-norm.  An important result is Lemma \ref{ranknonexp} which will imply that (\ref{rankeq1}) is sup-contractive.  In Section 3, we show that the general class of sup-contractive recurrence equations converges as desired.  Finally, in Section 4, we return to the max-type equation (\ref{maxeq}), and find the value of the asymptotically convergent periodic orbit for some specific values of M and P.

\begin{section}{Sup-Contractive Functions}
The convergence proofs in the next section apply to a wide class of functions $G: \mathbb{R}^a \to \mathbb{R}^b$, that includes compositions of contractive functions as in (\ref{alpha}) with the max and $k$-rank functions. We will be interested in the operator norm of functions $G$, where the norm used is  the sup norm. Namely, assume there exists $L \geq 0$ such that
\[ ||G(x) - G( y)||_{\infty} \leq L ||x - y||_{\infty} \]
for all $ x, y \in \mathbb{R}^a$. We can think of $L$ as the Lipschitz constant of $G$ in the sup norm, and will refer to it as the {\it sup-Lipschitz constant} in the following sections. When $L = 1$ we will call $G$ {\it sup-non-expansive} and when $L<1$ we will call $G$ {\it sup-contractive}.  Note that a contraction on $\mathbb{R}$ is sup-contractive since the infinity norm on $\mathbb{R}$ is simply the absolute value.

\begin{example} Let $G:\mathbb R^M\to \mathbb R$ be the function $G(x) = c||x||_p$ where $c >0$. Since
\[ |c ||x||_p - c ||y||_p| = c| ||x||_p - ||y||_p | \leq c||x-y||_p \leq c M^{1/p} ||x-y||_{\infty}, \]
$G$ is sup-non-expansive when  $c = M^{-1/p}$ and sup-contractive when  $c < M^{-1/p}$.
\end{example}

\begin{example} For a fixed $z \in \mathbb{R}^M$, define $G:\mathbb R^M\to \mathbb R$ by $G(x) = z^T x$. Since
\[ | z^T x - z^T y| = | z^T( x- y)| = \sum_{i=1}^M |z_i||x_i-y_i| \leq || x-  y||_{\infty} \sum_{i=1}^M |z_i| = || x -  y||_{\infty} || z||_1, \]
$G$ is sup-non-expansive when $|| z||_1= 1$, and sup-contractive when $|| z||_1 < 1$.
\end{example}

\begin{example} For $1\leq k\leq M$, let $R_k : \mathbb{R}^M \to \mathbb{R}$ be the function that returns the $k$-th largest entry of the $M$-dimensional input vector.  Customarily, $R_k$ is called the ``$k$-rank" function. It includes the $\max$ $(k=1)$ and $\min$ $(k=M)$  as special cases. Not surprisingly, the max function is sup-non-expansive. Somewhat more surprising is that this property holds for all $k$, as shown in the next Lemma.
\begin{lemma} \label{ranknonexp} The $k$-rank function $R_k$ is sup-non-expansive.
\end{lemma}
\begin{proof} Let $ x, y \in \mathbb{R}^M$, and assume, without loss of generality, that $R_k( y) \leq R_k( x)$.  Set $y_{r(i)}$ be the $i$-th largest component of $ y$ and let $x_{s(i)}$ be the $i$-th largest component of $ x$.  Then:
\[ y_{r(1)} \leq y_{r(2)} \leq \cdots \leq y_{r(k)} = R_k( y) \leq R_k( x) = x_{s(k)} \leq x_{s(k+1)} \leq \cdots \leq x_{s(M)} \]
Now examine the list of natural numbers $I = (r(1),...,r(k),s(k),...,s(M))$.  We note that $I$ has length $M+1$ but each entry is chosen from $\{1,...,M\}$.  Thus by the pigeonhole principle at least two of the listed numbers must be the same.  Note that all the $r(1),...,r(k)$ are distinct and all the $s(k),...,s(M)$ are distinct, thus there must be some $1\leq i \leq k$ and some $k \leq j \leq M$ such that $r(i) = s(j) = t$.  Therefore
\[ y_t = y_{r(i)} \leq R_k( y) \leq R_k( x) \leq x_{s(j)} = x_t, \]
which implies that
\[ |R_k( x) - R_k( y)| = R_k( x) - R_k( y) \leq x_t - y_t = |x_t - y_t| \leq  ||x- y||_{\infty} \]
Since the absolute value is the infinity norm on $\mathbb{R}$,
 $R_k$ is sup-non-expansive.
 \end{proof}
\end{example}

The sup-Lipschitz constants are multiplied under composition, based on the next Lemma.
\begin{lemma}\label{supcompositions} Assume $f_1,...,f_k: \mathbb{R}^M \to \mathbb{R}$ have sup-Lipschitz constants less than $L_1$ and $f: \mathbb{R}^k \to \mathbb{R}$ has sup-Lipschitz constant $L_2$, then $g( x) = f(f_1( x),f_2( x),...,f_k( x))$ has sup-Lipschitz  constant no larger than $L_1L_2$.
\end{lemma}
\begin{proof} Let $ x,  y \in \mathbb{R}^M$ then:
\begin{eqnarray*} |g( x) - g( y)| &=& | f(f_1( x),...,f_k( x))- f(f_1( y),...,f_k( y))| \\
 &\leq& L_2 \max_{1\leq i \leq k} |f_i( x) - f_i( y)| \\ &\leq& L_1L_2 ||x-y||_{\infty}
\end{eqnarray*}
\end{proof}
Note that if $f$ is sup-contractive, and all $f_i$ are sup-non-expansive then $g$ is sup-contractive.  Similarly if $f$ is sup-non-expansive, and all $f_i$ are sup-contractive then $g$ is sup-contractive.  Finally, the next Lemma shows that if we bundle together functions into a vector, the sup-Lipschitz constant cannot grow.
\begin{lemma} \label{lem26} Assume $f_1,...,f_k:\mathbb{R}^M \to \mathbb{R}$ have sup-Lipschitz constants at most $L$.  Then $f:\mathbb{R}^M\to \mathbb{R}^k$ defined by
\[ f( x) = (f_1( x),f_2( x),...,f_k( x))\] has sup-Lipschitz constant at most $L$.
\end{lemma}
\begin{proof} Let $ x, y \in \mathbb{R}^M$ then:
\[ || f( x) -  f( y)||_{\infty} = \max_{1 \leq i \leq k} \{ |f_i( x) - f_i( y)|\} \leq L || x -  y||_{\infty} \]
\end{proof}
Note that if $k=M$ and $L<1$ then in fact $ f$ is a contraction on $\mathbb{R}^M$, so bundling  functions is a way to build contractions.  Finally, note that these statements could be generalized to allow $f_i$ to be vector valued but this is not necessary for what follows.
\end{section}

\begin{section}{General Case}
We now return to our original recurrence equation
\[ x_n = \underset{1\leq i \leq M}{\textup{$k$-rank}} \{f_i(x_{n-i},n) \}\]
under the assumption that each $f_i:\mathbb{R}\to \mathbb{R}$ is a contraction.
Since $\textup{$k$-rank}$ is sup-non-expansive by Lemma \ref{ranknonexp}, and each $f_i$ is sup-contractive, the composition of the two is sup-contractive by Lemma \ref{supcompositions}.  This is a special case of the equation
\[ x_n = G_{n}(x_{n-1},...,x_{n-M}) \]
where $G_1,...G_{P}$ are any sup-contractive functions and $G_{n+P} = G_n$ for all $n$.  Setting $s = PM$, we will first show that for all initial conditions the limit is periodic of period $s$ by a contraction mapping argument.  We will then use the $P$-periodicity of $G_n$ to show that the limit is actually periodic of period $P$.

With slight abuse of notation we can consider $G_n$ as a function of all $s$ variables, although it depends only on the first M:
\begin{equation}\label{Gdef} x_n = G_{n}(x_{n-1},...,x_{n-s})= G_{n}(x_{n-1},...,x_{n-M}) \end{equation}
where $G_{n+s} = G_{n+MP} = G_n$ for all $n$.  We define
\[ \ov x_n = (x_{ns+1},...,x_{ns+s}) \]
and let $\ov x_0$ be the initial condition.  We want to show that we can write
\[ \ov x_{n+1} = \ov F(\ov x_n) \]
where $\ov F: \mathbb{R}^s \to \mathbb{R}^s$ is a contraction in the infinity norm on $\mathbb{R}^s$.  Let $\ov y = (y_1,...,y_s)$ and define the following functions for $k = 1,...,s$.
\begin{eqnarray*} F_1(\ov y) &=& G_1(y_{s},y_{s-1},...,y_{1}) \\
 F_k(\ov y) &=& G_k(F_{k-1}(\ov y),...,F_{1}(\ov y),y_s,y_{s-1},...,y_k) \end{eqnarray*}
Note that this is an inductive definition since $F_k$ depends on $F_1,...,F_{k-1}$.  Finally, we define
\begin{equation}\label{Fdef} \ov F(\ov y) = (F_1(\ov y),...,F_s(\ov y))\end{equation}
The recursive nature of the definition of $\ov F$ requires the following Lemma to show that $\ov{F}$ encapsulates the evolution of the sequence $\{x_n\}$.
\begin{lemma} Let $G_1,...,G_s:\mathbb{R}^s \to \mathbb{R}$ and let $\{x_i\}$ be defined by \eqref{Gdef}.  Define $\overline{F}$ by \eqref{Fdef}.  Then for all $n$ we have $\ov x_{n+1} = \ov F(\ov x_n)$. \end{lemma}
\begin{proof} This is equivalent to showing that $x_{ns+s+k} = F_k(\ov x_n)$ for $k = 1,...,s$ which is equivalent to:
\[ F_k(\ov x_n) = G_k(x_{ns+s+k-1},...,x_{ns+k}) \]
Note that when $k=1$ we have:
\[ F_1(\ov x_n) = G_1(x_{ns+s},...,x_{ns+1}) = x_{ns+s+1} \]
By definition of $F_1$ and $G_1$.  Then for $2 \leq k \leq s$ we can proceed by induction.  Assuming $x_{ns+s+\overline{k}} = F_{\ov k}(\ov x_n)$ for all $\ov k < k$ we have:
\begin{eqnarray*} F_k(\ov x_n) &=& G_k( F_{k-1}(\ov x_n),...,F_{1}(\ov x_n) ,x_{ns+s},...,x_{ns+k} ) \\
 &=& G_k(x_{ns+s+k-1},...,x_{ns+s+1} ,x_{ns+s},...,x_{ns+k} ) =  x_{ns+s+k}
 \end{eqnarray*}
Which completes the proof.
\end{proof}
Now that we have defined $\ov F$ we need to show that it is a contraction.
\begin{theorem}\label{contractiontheorem} Let $G_1,...,G_s:\mathbb{R}^s \to \mathbb{R}$ be sup-contractive and define $\ov{F}$ as in \eqref{Fdef}.  Then $\ov F: \mathbb{R}^s \to \mathbb{R}^s$ is a contraction with respect to the infinity norm.
\end{theorem}
\begin{proof} By Lemma \ref{lem26} it suffices to show that each $F_i$ is sup-contractive.  Note that the projection function $\pi_i(\ov y) = y_i$ is sup-non-expansive since:
\[ |\pi_i(\ov x) - \pi_i(\ov y)| = |x_i - y_i| \leq ||\ov x - \ov y||_{\infty} \]
Thus $F_1 = G_1(\pi_s(\ov y),...,\pi_1(\ov y))$ is sup-contractive by Lemma \ref{supcompositions} because $G_1$ is sup-contractive and $\pi_i$ is sup-non-expansive.  For $1<i\leq s$ we proceed by induction.  Assume $F_j$ is sup-contractive for all $j < i$, then:
\[ F_i(\ov y) = G_i(F_{i-1}(\ov y),...,F_{1}(\ov y),\pi_s(\ov y), \pi_{s-1}(\ov y),...,\pi_i(\ov y)) \]
So $F_i$ is a composition of the sup-contractive function $G_i$ with sup-contractive functions $F_{i-1},...,F_1$ and sup-non-expansive projection functions $\pi_s,...,\pi_i$.  Thus by Lemma \ref{supcompositions}, $F_i$ is sup-contractive, which completes the induction and thus shows that $\ov F$ is a contraction.
\end{proof}
\begin{corollary}\label{PMperiodicity}  Let $G_1,...,G_P:\mathbb{R}^M \to \mathbb{R}$ be sup-contractive and let $G_{n+P} = G_P$ for all $n$.  Given an initial condition $(x_1,...,x_M)$ let
\[ x_n = G_n(x_{n-1},...,x_{n-M}) \]
then $x_n$ converges to a unique $PM$-periodic orbit independent of initial conditions.
\end{corollary}
\begin{proof} By Theorem \ref{contractiontheorem} we can construct a contraction mapping $\ov F: \mathbb{R}^s \to \mathbb{R}^s$ such that $\ov x_{n+1} = \ov F(\ov x_n)$.  Thus, by the Contraction Mapping Theorem, $\ov F$ has a unique fixed point $x^* \in \mathbb{R}^s$ and for any $\ov x_0 \in \mathbb{R}^s$ we have $\lim_{n\to\infty}\ov x_n = x^*$.
\end{proof}
Note that $s$ may not be the prime period (smallest possible period).  So we now show that in fact $x_n$ must be asymptotically periodic of period $P$.
\begin{theorem}\label{Pperiodicity} Let $G_1,...,G_P:\mathbb{R}^M \to \mathbb{R}$ be sup-contractive and define $\ov{F}$ as in \eqref{Fdef}.  Let $x^*$ be the unique fixed point of $\ov F$.  Then $ x^*$ is periodic of period $P$.
\end{theorem}
\begin{proof} Define a shift operator by $S(\ov x_n) = (x_{ns+1+P},...,x_{ns+s+P})$.  We will use the fact that $G_i$ is actually periodic of period $P$ to show that the function $S$ commutes with $\ov F$.  Note that
\begin{eqnarray*} S(\ov F(\ov x_n)) &=& S(\ov x_{n+1}) = (x_{ns+s+1+P},...,x_{ns+s+s+P})  \\
  \ov F(S(\ov x_n)) &=& \ov F(x_{ns+1+P},...,x_{ns+s+P})
\end{eqnarray*}
First we examine the first component of $\ov F(S(\ov x_n))$:
\begin{eqnarray*} ( \ov F(S(\ov x_n)))_1 &=& G_1(x_{ns+s+P},...,x_{ns+1+P}) \\
 &=& G_{1+P}(x_{ns+s+P},...,x_{ns+1+P}) = x_{ns+s+1+P}
 \end{eqnarray*}
where the second equality comes from the fact that $G_i$ is periodic of period P.  This shows that the first components of $S(\ov F(\ov x_n))$ and $\ov F(S(\ov x_n))$ are the same.  We proceed inductively to show that all the components are the same.  Assume that $(S(\ov F(\ov x_n)))_j = (\ov F(S(\ov x_n)))_j$ for all $j<i$.  Then:
\begin{eqnarray*} ( \ov F(S(\ov x_n)))_i &=& G_i(F_{i-1}(S(\ov x_n)),...,F_1(S(\ov x_n)),x_{ns+s+P},...,x_{ns+s+P + i-s}) \\
&=& G_{i+P}(x_{ns+s+i-1+P},...,x_{ns+s+1+P},x_{ns+s+P},...,x_{ns+s+P + i-s}) \\
&=& x_{ns+s+i+P}
\end{eqnarray*}
where we have again used the P-periodicity if $G_i$ to conclude $G_i = G_{i+P}$.  This shows that
\[ S(\ov F(\ov x_n)) =  \ov F(S(\ov x_n)) \]
So inductively we have $S(\ov x_n) = \ov F^n(S(\ov x_0))$.  Let $x^*$ be the unique fixed point of $\ov F$ and define two sequences, the first with $x_0 = x^*$ and the second with $y_0 = S(x^*)$.  Note that
\[  \lim_{n\to\infty} y_n = \lim_{n\to\infty} \ov F^n(\ov y_0) = x^* \]
since all initial conditions converge to $x^*$, and at the same time:
\[  \lim_{n\to\infty} y_n = \lim_{n\to\infty} \ov F^n(\ov y_0) = \lim_{n\to\infty} \ov F^n(S(x_0)) = \lim_{n\to\infty} \ov S(\ov F_n(x_0))= \lim_{n\to\infty} S(x_0) = S(x^*) \]
So we conclude that $S(x^*) = x^*$ and thus $x^*$ is periodic of period $P$.
\end{proof}
\begin{corollary}\label{mainresult} Let $G_1,...,G_P:\mathbb{R}^M \to \mathbb{R}$ be sup-contractive and let $G_{n+P} = G_P$ for all $n$.  Given an initial condition $(x_1,...,x_M)$ let
\[ x_n = G_n(x_{n-1},...,x_{n-M}) \]
then $x_n$ converges to a unique $P$-periodic orbit independent of initial conditions.
\end{corollary}
\begin{proof} By Theorem \ref{Pperiodicity} there exists a unique $x^* \in \mathbb{R}^s$ which is $P$-periodic such that $\lim_{n\to\infty} \ov x_n = x^*$.  Thus $x_n$ is asymptotically periodic of period P.
\end{proof}
We conclude that $x_n$ approaches a unique periodic orbit, for any initial condition, whose period is equal to the forcing period $P$.  The periodicity of the rank-type equation (\ref{rankeq1}) is now an easy Corollary.
\begin{corollary}\label{rankResult} For $i=1,...,M$ let $f_i(x,n):\mathbb{R}\times\mathbb{N}\to\mathbb{R}$ be contractive in $x$ and $P$-periodic in $n$.  Given an initial condition $(x_1,...,x_M)$ and $k \in \{1,...,M\}$ let
\[ x_n = \textup{k-rank}\{f_i(x_{n-i},n)\} \]
then $x_n$ converges to a unique $P$-periodic orbit independent of initial conditions.
\end{corollary}
\begin{proof} Recall that by Lemma \ref{ranknonexp}, the $\textup{k-rank}$ function is sup-non-expansive, and each $f_i(x,n)$ is sup-contractive in $x$ so by Lemma \ref{supcompositions}, the composition is sup-contractive.  By Corollary \ref{mainresult}, $x_n$ is asymptotically periodic of period P.
\end{proof}
We can now return to the equation from Example \ref{ex3}.
\begin{example} \label{ex37}
Let $f_1,...,f_M : \mathbb{R} \to \mathbb{R}$ be contractions.  Let $P \leq M$ be a positive integer and set $\ov{n}=1+(n\mod P)$.  Let
\[ x_n = G_n(x_{n-1},...,x_{n-M}) = \frac{1}{2} \left( \max\{f_i(x_{n-i})\} - \overline{n}\textup{-rank}\{f_i(x_{n-i})\} \right) \]
Recall that by Lemma \ref{ranknonexp} the rank functions are all sup-non-expansive.  Thus $\ov n\textup{-rank}\{f_i(x_{n-i})\}$ is a composition of a sup-non-expansive function with the sup-contractive functions $f_i$, and thus the composition is sup-contractive by Lemma \ref{supcompositions}.  Furthermore, setting $\ov z = (1/2,-1/2)$, we see that $||\ov z||_1 = 1$ so $f(x) = z^T x$ is sup-non-expansive. Therefore
\begin{eqnarray*} G_n(y_1,...,y_M) &=& z^T \left(\max\{f_i(y_{i})\},\overline{n}\textup{-rank}\{f_i(y_{i})\right) \\
&=& \frac{1}{2} \left( \max\{f_i(y_{i})\} - \overline{n}\textup{-rank}\{f_i(y_{i})\} \right)
\end{eqnarray*}
is sup-contractive for all $n$.  By Corollary \ref{mainresult}, $x_n$ is asymptotically periodic with period $P$.
\end{example}
\end{section}

\begin{section}{Finding the Periodic Limit}
We now return to the max-type equation \eqref{maxeq} and rank-type equation \eqref{rankeq1} and attempt to find closed-form solutions.  The closed-form solution to the autonomous contractive rank-type equation was first given in \cite{sauerRank}, and we are able to reprove this result as a special case.  However, we will see that finding a closed formula for the limit under periodic forcing is in general more difficult.
We will need two lemmas about contractions on $\mathbb{R}$.
\begin{lemma} \label{lem1} Let $f: \mathbb{R} \to \mathbb{R}$ be a contraction with fixed point $r$ then $x > r$ implies $f(x) < x$, and $x < r$ implies $f(x) > x$.
\end{lemma}
\begin{proof}
Since $f(r)=r$ and $f$ is a contraction, $|f(x)-r|<|x-r|$. If $x>r$, then
\[f(x)-r\leq |f(x)-r|<|x-r|=x-r,\]
so $f(x)<x$. If $x<r$, then
\[r-f(x)\leq |f(x)-r|<|x-r|=r-x,\]
so $f(x)>x$.
\end{proof}
\begin{lemma} \label{lem21} Let $f,g: \mathbb{R} \to \mathbb{R}$ be contractions, and assume that $r_1, r_2, r_3, r_4$ satisfy $f(r_2)=r_1<r_3$ and $g(r_1)=r_2<r_4$. Then either $f(r_4)<r_3$ or $g(r_3)<r_4$.
\end{lemma}
\begin{proof}
Assume $g(r_3)\geq r_4$. Then
\[ r_4-r_2=r_4-g(r_1)\leq g(r_3)-g(r_1) < r_3-r_1.\]
The contractivity of $f$ yields
\[f(r_4)-r_1\leq |f(r_4)-r_1|=|f(r_4)-f(r_2)|< |r_4-r_2|=r_4-r_2< r_3-r_1\]
which implies that $f(r_4)<r_3$.
\end{proof}

First, we can use Lemma \ref{lem1} to find the explicit solution of the autonomous equation \eqref{rankeq} where each $f_i$ is a contraction with fixed point $r_i$.  The following result, first proved in \cite{sauerRank}, shows that every initial condition converges to the constant solution $\textup{$k$-rank}\{r_i\}$.
\begin{theorem}\label{periodone}
 Let $f_i:\mathbb{R}\to\mathbb{R}$ be a contraction with fixed point $r_i$ for $1\leq i \leq M$ and set $ x_{n} = \textup{$k$-rank} \{ f_i(x_{n-i})\}$.  Then $\lim_{n\to\infty} x_n = \textup{k-rank} \{r_i\}$.
\end{theorem}
\begin{proof} Note that this recurrence is autonomous and thus $P=1$ so by Corollary \ref{rankResult} every initial condition converges to a unique constant (period one) solution $x^*$.  Thus it remains only to show that $x^* = \textup{$k$-rank}\{r_i\}$ is a fixed point of the recurrence.  Let $\sigma: \{1,...,M\} \to \{1,...,M\}$ be a permutation such that:
\[ r_{\sigma(1)} \leq r_{\sigma(2)} \leq \cdots \leq r_{\sigma(M)} \]
So the constant solution should be $x^* = r_{\sigma(k)}$.  Assume $x_1 = x_2 = \cdots = x_M = x^*$.  Then:
\begin{eqnarray*} x_{M+1} &=& \textup{$k$-rank}\{f_1(x^*),...,f_{\sigma(k)}(x^*),...,f_M(x^*)\}  \\
&=& \textup{$k$-rank}\{f_1(x^*),...,x^*,...,f_M(x^*)\}
\end{eqnarray*}
We want to show that $x_{M+1} = x^*$.  Note that when $i < k$ we have $r_{\sigma(i)} \leq r_{\sigma(k)} = x^*$ so by Lemma 3.1 we have $f_{\sigma(i)}(x^*) \leq x^*$.  Similarly when $i > k$ we have $r_{\sigma(i)} \geq r_{\sigma(k)} = x^*$ so by Lemma 3.1 we have $f_{\sigma(i)}(x^*) \geq x^*$.  Thus we have:
\[ f_{\sigma(1)}(x^*),...,f_{\sigma(k-1)}(x^*) \leq x^* \leq f_{\sigma(k+1)}(x^*),...,f_{\sigma(M)}(x^*) \]
Thus $x_{M+1} = x^*$ so  $x^*$ is a fixed point of the recurrence and therefore it is the unique limit for any initial condition.
\end{proof}
Theorem \ref{periodone} gives the complete solution to the rank-type equation \eqref{rankeq1} in the period one case.  We now turn to the case of period two forcing and restrict our attention to the max-type equation \eqref{maxeq}, for which it is possible to find a closed-form solution.    This solution gives insight into the complexity of solutions to \eqref{maxeq} when the forcing period is large.
For $M=P=2$,  the equation (\ref{rankeq1}) becomes the period two recurrence
\begin{eqnarray}\label{tte1}
 x_{2i} &=& \max \{ f_1(x_{2i-1}) , f_2(x_{2i-2}) \} \nonumber \\
 x_{2i+1} &=& \max \{ g_1(x_{2i}) , g_2(x_{2i-1}) \}
\end{eqnarray}
where $f_1, f_2, g_1$ and $g_2$ are contractions.

We can denote the fixed points
\begin{eqnarray*}
 f_1(g_1(r_1)) &=& f_1(r_2) = r_1 \\
 g_1(f_1(r_2)) &=& g_1(r_1) =  r_2 \\
  f_2(r_3) &=& r_3 \\
  g_2(r_4) &=& r_4.
 \end{eqnarray*}

\begin{theorem} \label{explicit} The period-two orbit
\begin{eqnarray} \label{maxsol}
x_{2i} &=& \max\{f_1(\max\{r_2,r_4\}),r_3\}\nonumber \\
x_{2i+1} &=& \max\{g_1(\max\{r_1,r_3\}),r_4\}
\end{eqnarray}
for $i\ge 0$ is the attracting period-two orbit of difference equation (\ref{tte1}).
\end{theorem}

\begin{proof} Due to Corollary \ref{rankResult}, it suffices to show that the formula (\ref{maxsol}) gives a period-two orbit. First note the following table can be obtained easily:
\begin{center}
  \begin{tabular}{ | c | c || c | c |}
    \hline
    $\max\{r_1,r_3\}$ & $\max\{r_2,r_4\}$ & $x_{2i}$ & $x_{2i+1}$ \\ \hline\hline
    $r_1$ & $r_2$ & 					$r_1$ & $r_2$ \\ \hline
    $r_3$ & $r_4$ & 					$\max\{r_3,f_1(r_4)\}$ &  $\max\{r_4,g_1(r_3)\}$ \\ \hline
    $r_1$ & $r_4$ & 					$\max\{r_3,f_1(r_4)\}$ & $r_4$ \\ \hline
    $r_3$ & $r_2$ & 					$r_3$ & $\max\{r_4,g_1(r_3)\}$ \\
    \hline
  \end{tabular}
\end{center}
For example, in the first row we have:
\[ x_{2i} = \max\{f_1(\max\{r_2,r_4\}),r_3\}  = \max\{f_1(r_2),r_3\} = \max\{ r_1, r_3 \} = r_1 \]
The rest follow from similar simple logic.  Now, we can make a new table considering all the possibilities for $x_{2i}$, $x_{2i+1}$, and $x_{2i+2} = \max\{f_1(x_{2i+1}),f_2(x_{2i})\}$.
\begin{center}
  \begin{tabular}{ | c | c || c | c || c |}
    \hline
    $x_{2i}$ & $x_{2i+1}$ &				$f_1(x_{2i+1})$ &		$f_2(x_{2i})$ &	$x_{2i+2}$ \\ \hline\hline
    $r_1$ & $r_2$ &				$r_1$ &			$f_2(r_1)$ &	$r_1$ \\ \hline
    $r_3$ & $r_4$ &		 		$f_1(r_4)$ &		$r_3$ &		 $r_3$\\ \hline
    $f_1(r_4)$ & $r_4$ &			$f_1(r_4)$ &		$f_2(f_1(r_4))$& $f_1(r_4)$ \\ \hline
    $r_3$ & $ g_1(r_3)$ &			$f_1(g_1(r_3))$ &	$r_3$&		$r_3$ \\ \hline
  \end{tabular}
\end{center}

 Since columns three and four are clear, we must justify the final column using Lemma \ref{lem1}, (note that the final column is simply the max of columns three and four).  In the first row, since $r_1 > r_3$ and $r_3$ is the fixed point of $f_2$, Lemma \ref{lem1} implies that $f_2(r_1) < r_1$.  In the second row, note that $r_3  = \max\{ f_1(r_4), r_3\}$ so $r_3 > f_1(r_4)$.  In the third row, note that $f_1(r_4)  > r_3$ and since $r_3$ is the fixed point of $f_2$ by Lemma \ref{lem1} we conclude that $f_2(f_1(r_4)) < f_1(r_4)$.  Finally, in the fourth column note that this combination only occurs when $r_3 > r_1$ and $r_1$ is the fixed point of $f_1 \circ g_1$ so by Lemma \ref{lem1} we have $f_1(g_1(r_3)) < r_3$.  Thus the table shows that $x_{2i+2}=x_{2i}$.

Note that one combination, $x_{2i} = f_1(r_4)$ and $x_{2i+1} = g_1(r_3)$ is missing from the table. If this combination occurred, then $f_1(r_2)= r_1<r_3\leq f_1(r_4)$ and $ g_1(r_1)=r_2 < r_4 \leq g_1(r_3)$, which contradicts Lemma \ref{lem21}.  Thus this combination is impossible, and the four rows of the table represent all possibilities.

We now construct an analogous table to show that $x_{2i+3}=x_{2i+1}$, which will complete the proof.
\begin{center}
  \begin{tabular}{ | c | c || c | c || c |}
    \hline
    $x_{2i+2}$ & $x_{2i+1}$ &				$g_1(x_{2i+2})$ &		$g_2(x_{2i+1})$ &	$x_{2i+3}$ \\ \hline\hline
    $r_1$ & $r_2$ &				$r_2$ &			$g_2(r_2)$ &	$r_2$ \\ \hline
    $r_3$ & $r_4$ &		 		$g_1(r_3)$ &		$r_4$ &		 $r_4$\\ \hline
    $f_1(r_4)$ & $r_4$ &			$g_1(f_1(r_4))$ &		$r_4$& $r_4$ \\ \hline
    $r_3$ & $ g_1(r_3)$ &			$g_1(r_3)$ &	$g_2(g_1(r_3))$&		$g_1(r_3)$ \\ \hline
  \end{tabular}
\end{center}
In the first row, $r_2 > r_4$ so by Lemma \ref{lem1}, $g_2(r_2) < r_2$.  In the second row, $r_4  = \max\{g_1(r_3),r_4\}$ so $r_4 > g_1(r_3)$.  In the third row, $r_4 > r_2$ so by Lemma \ref{lem1} $g_1(f_1(r_4) < r_4$.  In the final row, $g_1(r_3) > r_4$ so by Lemma \ref{lem1} we have $g_2(g_1(r_3)) < g_1(r_3)$.  This completes the justification of this table, and finishes the proof.
\end{proof}
\begin{example}
\label{ex44}
We return to the equation from Example \ref{ex1} with $M = P = 2$:
\begin{eqnarray*}
x_{2i}&=& \max\{A_{11}x_{2i-1}^{\alpha_1},A_{21}x_{2i-2}^{\alpha_2} \}\\
x_{2i+1} &=& \max\{A_{12}x_{2i}^{\alpha_1},A_{22}x_{2i-1}^{\alpha_2} \}
\end{eqnarray*}
where $A_{jk} > 0$ and $-1<\alpha_1,\alpha_2<1$. We rewrite this equation by taking the natural log of each term (since $\ln$ is monotonic) to get
\begin{eqnarray*}
y_{2i}&=& \max\{\ln A_{11}+\alpha_1y_{2i-1},\ln A_{21}+\alpha_2 y_{2i-2}\}\\
y_{2i+1} &=& \max\{\ln A_{12}+\alpha_1 y_{2i},\ln A_{22}+\alpha_2 y_{2i-1} \}.
\end{eqnarray*}
In terms of Theorem 4.2 we have:
\begin{eqnarray*}
r_1 &=& \frac{\ln A_{12} + \alpha_1 \ln A_{11}}{1-\alpha_1^2}\\
r_2 &=& \frac{\ln A_{11} + \alpha_1 \ln A_{12}}{1-\alpha_1^2}\\
r_3 &=& \frac{\ln A_{22}}{1-\alpha_2}\\
r_4 &=& \frac{\ln A_{21}}{1-\alpha_2}
\end{eqnarray*}
Then the period two limit is  defined in Theorem 4.2, and we can exponentiate to get the period two limit of $x_n$.  Thus we have
\begin{eqnarray*}
x_{2i} &=& \max \left\{ A_{12}A_{21}^{\frac{\alpha_1}{1-\alpha_2}},A_{12}A_{11}^{\frac{\alpha_1}{1-\alpha_1^2}}A_{12}^{\frac{\alpha_1^2}{1-\alpha_1^2}}, A_{22}^{\frac{1}{1-\alpha_2}} \right\} \\
x_{2i+1}&=& \max \left\{ A_{11}A_{22}^{\frac{\alpha_1}{1-\alpha_2}},A_{11}A_{12}^{\frac{\alpha_1}{1-\alpha_1^2}}A_{11}^{\frac{\alpha_1^2}{1-\alpha_1^2}}, A_{21}^{\frac{1}{1-\alpha_2}} \right\}
\end{eqnarray*}
as the periodic limit.  This solution is consistent with the autonomous case  in Example 2.6 of \cite{sauerMax}. Setting $A_1=A_{11}=A_{12}$ and $A_2=A_{21}=A_{22}$, the solution is found to be simply $x_n \to \max\left\{ A_1^{\frac{1}{1-\alpha_1}}, A_2^{\frac{1}{1-\alpha_2}} \right\}$.
\end{example}
\end{section}

\begin{section}{Discussion}
We have shown in Corollary \ref{rankResult} that solutions of periodically-forced rank-type difference equations are asymptotically periodic of the forcing period. The same is true of a class of more general equations called sup-contractive, according to Corollary \ref{mainresult}. In some simple cases, we were able to identify explicit solutions, as in Theorem \ref{periodone} and Theorem \ref{explicit}. The solutions appear to be significantly more complicated for larger period $P$ and memory $M$ than treated here, and explicit formulas for the solutions remain to be found.

\bigskip

{\bf Acknowledgments.} This research was partially supported by NSF under DMS-0508175 and DMS-0811096.

\end{section}


text:
\begin{document}

\title{Inferring Causal Networks of Dynamical Systems \\ through Transient Dynamics and Perturbation}

\author{George Stepaniants}
\email[]{gstepan@mit.edu}

\affiliation{Department of Mathematics, Massachusetts Institute of Technology, Cambridge, MA\\
Department of Mathematics, University of Washington, Seattle, WA}

\author{Bingni W. Brunton}
\email[]{bbrunton@uw.edu}

\affiliation{Department of Biology, University of Washington, Seattle, WA}

\author{J. Nathan Kutz}
\email[]{kutz@uw.edu}

\affiliation{Department of Applied Mathematics, University of Washington, Seattle, WA}

\begin{abstract}

Inferring causal relations from time series measurements is an ill-posed mathematical problem, where typically an infinite number of potential solutions can reproduce the given data. We explore in depth a strategy to disambiguate between possible underlying causal networks by perturbing the network, where the actuations are either targeted or applied at random. The resulting transient dynamics provide the critical information necessary to infer causality. Two methods are shown to provide accurate causal reconstructions:  Granger causality (GC) with perturbations, and our proposed {\em perturbation cascade inference} (PCI). Perturbed GC is capable of inferring smaller networks under low coupling strength regimes. Our proposed PCI method demonstrated consistently strong performance in inferring causal relations for small (2--5 node) and large (10--20 node) networks,  with both linear and nonlinear dynamics. Thus the ability to apply a large and diverse set of perturbations/actuations to the network is critical for successfully and accurately determining causal relations and disambiguating between various viable networks.

\end{abstract}

\maketitle

\section{Introduction}

The ability to determine causal relationships in complex, dynamical networks from time-series measurements alone is an important open challenge in the engineering, biological and physical sciences.  This task is challenging because inferring causal networks from time series observations alone is an ill-posed mathematical problem, and a potentially infinitude of solutions may accurately reproduce the given data.  Despite decades of research effort, the large and diverse set of mathematical methods that have been developed still have limitations in accurately reproducing causal network structures, especially for nonlinear dynamical systems on networks~\cite{lusch2016, sugihara2012,monster2017,daniella,emily,krakovska,angulo2017}.

In certain applications, one has the ability to actuate the dynamical network of interest and generate additional information about the unknown network structure.  In this work, we exploit this capacity to perturb a dynamical network, in targeted or random ways, to extract an accurate directed graph of the true underlying dynamical system.  We demonstrate the accuracy and efficiency of our \emph{perturbation cascade inference} (PCI) algorithm on a variety of test problems, showing in which cases such network perturbation strategies can resolve this ill-posed inference problem.

Wiener first proposed a statistical notion of causality~\cite{wiener1956} by noting that $Y$ causes $X$ if knowing the past of $Y$ improves the prediction of $X$,  compared to using the past of $X$ alone.   Many innovations, modifications and reformulations of the causal inference problem have since been proposed~\cite{imbensrubin2015, pearl2009_1, morganwinship2014, holland1986}, with various statistical regularizations used to achieve potential solutions of this ill-posed problem.  The Nobel prize winning work of Granger~\cite{granger1969, granger1980} built upon Wiener's definition and formalized time series causality through linear regression of stochastic processes.  However, determining causal relationships in large networks of nonlinear dynamical systems still remains problematic with not only {\em Granger Causality} (GC)~\cite{lusch2016, sugihara2012}, but also with {\em convergent cross mapping} (CCM)~\cite{monster2017}, and other sparsity-promoting techniques (e.g. Ref.~\cite{daniella} Fig. 1, 2 and Ref.~\cite{emily} Fig. 4). Krakovsk\'a et al~\cite{krakovska} provide useful inference methods on a very limited set of systems and fail to produce even approximately accurate directed graphs for more complex network dynamics.

\begin{figure}[t]
\includegraphics[width=9cm]{Images/inferenceLimitations.jpeg}
\caption{One simple example to demonstrate the fundamentally ill-posed nature of network inference. Specifically, two Kuramoto networks can easily produce identical time-series data, and it is impossible to  disambiguate between them from this data alone.
We show the best-case predictions of GC. When nodes 2 and 3 are perturbed (i.e. their initial conditions are allowed to vary), GC is able to perfectly resolve both graphs.}
\label{fig:netinfProblems}
\end{figure}

\begin{figure*}[t]
\includegraphics[width=15cm]{Images/NetworkInferenceSummary.png}
\caption{An illustration of two approaches to network inference. Given the true network shown on the left, one approach is to passively collect time series observations of each network node (top center) and use time series analysis to infer the network structure shown on the right. If perturbations of the network nodes are allowed, another approach is to systematically force the network nodes and study their dynamics after each perturbation to infer the network structure.}
\label{fig:inferenceSummary}
\end{figure*}

Fig.~\ref{fig:netinfProblems} illustrates how inferring networks from observational data is an ill-posed problem. In particular, identical time series data can be collected from several different network topologies, making disambiguation of underlying network structures impossible on the given data.  In such cases, it is necessary to push the networked dynamical system into regimes where its dynamics are no longer degenerate; in other words, we need observations in regimes where it is not possible for two different networks to produce the same observed trajectories.

Leveraging the ability to actuate or perturb components of a networked system can lead to significant insights about its network topology and dynamics~\cite{timme2014}.  Perturbation inference methods subject a networked dynamical system to external driving forces and measure its collective response to reverse-engineer the underlying structure.  For instance, drug perturbations and gene knockouts are applied extensively in the study of gene regulatory networks~\cite{ud2016optimal, meinshausen2016methods}, and activations of genetically modified neurons in optogenetics enable the identification of specific neuronal connectivities~\cite{lepperod2018inferring}. A large body of research has also investigated inference of networked systems subject to small perturbations about a stable point, particularly in the context of gene regulatory network reconstruction~\cite{bansal2006, bansal2007, gardner2003}.  By observing response dynamics of a system from a series of systematic perturbations, various inference methods have been able to model the time dynamics and underlying network topology of coupled oscillator systems, cellular signaling networks and competitive economic markets~\cite{timme2007, molinelli2013, kuypers2012, delabays2020network}.

The present work discusses approaches to reverse engineer the network connectivity  by perturbing, or actuating, the networked dynamical system as shown in Fig.~\ref{fig:inferenceSummary}.  We use simulated nonlinear dynamical systems, for which the ground-truth network structure is known, to study the network reconstruction quality of the various inference methods in its passive, unperturbed state, and under the influence of actuation or perturbation.  We propose two methods for uniquely disambiguating data for inference:  (a) judiciously perturbed GC with sampling only of the transient dynamics, and (ii) an active inference approach called \textit{Perturbation Cascade Inference} (PCI).
PCI infers the structure of a dynamical network through systematic perturbations of its components.  This approach is similar in spirit to {\em perturb and observe} methods for system identification for control~\cite{billings1980identification}.   In order to evaluate the performance of network inference algorithms in both the passive and active regimes, we simulate respectively the unperturbed and perturbed time dynamics of coupled spring-mass and Kuramoto oscillators on random Erdos-Renyi networks with varying edge densities, coupling strengths and sizes.  We show that by either targeted perturbations, or random global network actuation, we can accurately reconstruct the causal network with a sufficiently large number of perturbations.

\section{Background}
Before developing the proposed PCI algorithm, we briefly review the mathematical architecture of one of the first time-series causal inference techniques.
This section also serves to establish some basic notions of causal relationships.
Granger Causality (GC) infers the causal relationship between two multivariate time series $\mathbf{X}_t, \mathbf{Y}_t$ by fitting a vector autoregressive (VAR) model to test if $\mathbf{X}_t$ conditional on its own past does not depend on the past of $\mathbf{Y}_t$. Let $\mathbf{X}_t, \mathbf{Y}_t$ be $n \times 1$ vectors. A $p$th order VAR model for $\mathbf{X}_t$ has the form
\begin{equation}\label{eq:1}
\mathbf{X}_t = \sum_{k=1}^p \mathbf{A}_k\mathbf{X}_{t - k} + \bm{\epsilon}_t,
\end{equation}
where the $n \times n$ matrices $\mathbf{A}_k$ are the regression coefficients that specify how $\mathbf{X}_t$ depends on its past and the $n \times 1$ vector $\bm{\epsilon}_t$ contains the residuals or prediction error of the model.  GC tests whether past values of $\mathbf{Y}_t$ improve the prediction of $\mathbf{X}_t$ as compared to the past values of $\mathbf{X}_t$ alone. The joint regression of $\mathbf{X}_t$ on its own past and the past of $\mathbf{Y}_t$ is given by
\begin{equation}\label{eq:2}
\mathbf{X}_t = \sum_{k=1}^p \mathbf{A}'_k\mathbf{X}_{t - k} + \sum_{k=1}^p \mathbf{B}'_k\mathbf{Y}_{t - k} + \bm{\epsilon}'_t,
\end{equation}
where $\mathbf{A}'_k, \mathbf{B}'_k$ are $n \times n$ matrices and the residual $\bm{\epsilon}'_t$ is an $n \times 1$ vector. Using the residuals from both VAR models in (\ref{eq:1}) and (\ref{eq:2}), we can construct the log-likelihood ratio
$\mathcal{F}_{\mathbf{Y} \to \mathbf{X}} = \ln ({|\Sigma'|}/{|\Sigma|})$
where $\Sigma' = \Cov(\bm{\epsilon}'_t)$ and $\Sigma = \Cov(\bm{\epsilon}_t)$ are the covariance matrices of the residuals of both VAR models. This log-likelihood ratio is called the \textit{G-causality} from $\mathbf{Y}$ to $\mathbf{X}$ and it characterizes the causal influence of $\mathbf{Y}$ on $\mathbf{X}$. $\mathcal{F}_{\mathbf{Y} \to \mathbf{X}}$ is the test statistic for the null hypothesis of zero causality
\begin{equation}\label{eq:4}
\mathbf{B}_1' = \mathbf{B}_2' = ... = \mathbf{B}_p' = 0.
\end{equation}
We say that $\mathbf{Y}$ \textit{Granger causes} $\mathbf{X}$ if the G-causality $F_{\mathbf{Y} \to \mathbf{X}}$ is statistically significant.

In unconditional G-causality described above, we can erroneously infer that $\mathbf{Y}$ causes $\mathbf{X}$ even if there is no causal relationship, since both variables are dependent on a third latent, confounding variable $\mathbf{Z}$. The method of \textit{conditional G-causality} eliminates such spurious causalities by ``conditioning out" all potential common dependencies. In the context of network inference, we can apply conditional G-causality to determine the existence of a causal relationship between every pair of nodes in a network. This approach is known as \textit{pairwise-conditional G-causality} and is the basis of the Multivariate Granger Causality (MVGC) toolbox developed by Barnett \& Seth \textit{et al.} which we use in this paper \cite{barnettseth2014}.

Results using GC inferences from the MVGC toolbox on time series observations of unforced mass-spring networks are shown in Fig.~\ref{fig:GCHarmonicUnforced}. Since GC performs an autoregressive linear fit to the observed time series data, it is expected to successfully infer the causal relationships in a linear dynamical system. Accordingly, GC makes highly accurate inferences on the mass-spring system for a wide range of network sizes when the coupling strengths (spring constants) are close to 1. The performance of this method rapidly degrades if the coupling strength is significantly smaller or larger than this value.  However, we are primarily interested in networked nonlinear dynamics, as most complex systems of interest are rarely linear.  Figure~\ref{fig:inferenceAccuracy} considers a network of nonlinear oscillators (Kuramoto oscillator network) and shows that MVGC fails to correctly identify the causal structure under most conditions~\cite{lusch2016}. This example motivates our need to develop a more robust method to extract causal relations.

\begin{figure}[t]
\includegraphics[width=9cm]{Images/GCHarmonicUnforced.jpeg}
\caption{\textbf{(a)} An example of a 5 node mass-spring network with the endpoints attached to two fixed walls. \textbf{(b)} The performance of GC on such mass-spring networks varies with network sizes and coupling strengths. For coupling strengths near 1, GC is capable of inferring the network structure of the mass-spring system with 80-85\% accuracy for network sizes between 2 to 10 nodes. GC is only able to achieve accuracy above 85\% when the network size is 2. Furthermore, if we increase or decrease the coupling strength of the system away from 1, the inference accuracy of GC drops below 80\% for mass-spring networks of all sizes.}
\label{fig:GCHarmonicUnforced}
\end{figure}

\begin{figure*}[t]
    \centering
    \vspace*{-.5in}
    \includegraphics[width=15cm]{Images/inferenceAccuracy.jpeg}
    \vspace*{-.6in}
    \caption{Comparison of Granger inferred networks to the true underlying Kuramoto oscillator network with varied simulation endtimes and coupling strengths. Seemingly large edge prediction accuracies do not resemble the true underlying network structure. Moreover, the eigenvalues of the network connectivity matrix are not well characterized. Therefore, it would not be possible to use an inferred network that is 75-80\% accurate to simulate the time dynamics of a system.}
    \label{fig:inferenceAccuracy}
\end{figure*}

There are a number of alternative approaches to GC for computing causality.  For instance, Convergent Cross Mapping ( CCM,~\cite{sugihara2012} assumes that the time-series data is sampled from a dynamical attractor and uses time-delay embeddings to construct attractor manifolds from the time-series observations of variables ${\bf X}$ and ${\bf Y}$. The correlation between the true state of ${\bf Y}$ and the predictions of ${\bf Y}$ given the ${\bf X}$-manifold determine the strength of the causal relationship from ${\bf Y}$ to ${\bf X}$.  Using information-theoretical approaches\cite{rahimzamani, krakovska}, it is possible to reconstruct causal graphs by computing statistics such as transfer entropy or directed information between the time series measurements of every pair of variables in a system. Similar to CCM, such approaches make no assumptions about the data generation process and often utilize time-lagged representations in the analysis.  Finally, a variety of model-based approaches have been proposed to fit the time series data through time-lagged regressions and ordinary differential equations. As described in \cite{daniella}, such methods often promote sparsity in the network parameters in order to restrict the space of possible solutions (e.g. obtain networks with the fewest number of edges that closely predict the system dynamics). Yet other methodologies such as Dynamic Causal Modeling \cite{friston} and Bayesian networks \cite{yu2004} have been developed to solve this ill-posed problem in application to neuronal dynamics and gene regulatory networks. Each method imposes a different regularization to extract a solution to the causal inference problem. However, these methods still face the same problem for causal inference in networked, nonlinear dynamics problems as illustrated in Fig.~\ref{fig:inferenceAccuracy}. In this paper, we show how the systematic perturbation of nodes in the network, either targeted or random, allow the resulting transient time dynamics to disambiguate the causal relations.

\section{Transients for Inference}
We will show how perturbations and transients can be used to reverse-engineer the structure of a time-dependent network.  Importantly, a distinct and diverse set of perturbations such as individual actuations of nodes or random kicks to the system are typically required in order to fully disambiguate the causal pathways.  Suppose we have a mass-spring network of five nodes with \textit{directed interactions} like the one depicted in Fig. \ref{fig:GCHarmonicUnforced}(a).

To learn about the connectivity of this system, we could perturb the blue node in the network away from its equilibrium position and observe how all the other oscillators respond to this perturbation over time.  After the initial mass is released a longitudinal wave travels through the network, displacing the masses and allowing us to observe the order in which they become displaced. Under the \textit{strict assumption} that all the coupling strengths between coupled oscillators are equal, we in fact observe that the red and yellow nodes become perturbed directly after the blue node, which subsequently activates the purple node. Note that there is no directed path from the blue to the green node, hence the green node remains unperturbed. This highlights, at least in part, the features exploited by algorithmic structures to infer causal relations.

\subsection{GC and Transients}

Our first use of perturbations of a networked dynamical systems is with the standard GC framework for inference.  But instead of simply giving GC time-series measurements, we perturb the system by giving it different starting \textit{initial conditions}, and use the transient information to disambiguate the structured network.  As has already been shown, GC by itself fails to produce meaningful results in nonlinear systems unless time series lengths and coupling strengths of the system are chosen judiciously (see Fig.~\ref{fig:inferenceAccuracy}).

In the case of the mass-spring system, the damped harmonic oscillators eventually tend to their equilibrium positions and in the case of the Kuramoto model, the oscillators synchronize and rotate with the same angular velocity. After energy imparted on the system is lost, the dynamics of the coupled oscillators become uninformative for the purposes of network reconstruction as shown in Fig. \ref{fig:kuramotoGCTransients}. The perturbation inference approach discussed here takes advantage of the transient time dynamics of a system directly after it is perturbed, while disregarding the complex nonlinear interactions that can arise later in time.

Fig.~\ref{fig:kuramotoGCTransients} shows an example of a networked, nonlinear dynamical system of Kuramoto oscillators.  We have already shown in Fig.~\ref{fig:inferenceAccuracy} that GC fails when simply presented with a given set of time-series data. But by considering the transient dynamics, and in fact, only using the first 5 seconds of transient dynamics, the GC accuracy is 85\% in reconstructing the true network. This shows that the transients contain most of the information required for inferring causality. By sampling for longer times, the dynamics fall into an attractor which obfuscates the true network architecture.  A more careful analysis of the Kuramato oscillator system is considered in Fig.~\ref{fig:kuramotoGCTransientsStrength}. The analysis of this plots shows the dependency of the sampling time on the coupling strength, which is directly related to the transient time exhibited in the dynamics.  Thus the time sampled of the dynamics, and its transient behavior, is critical in determining an accurate representation of causal relations. A  more detailed assessment of GC when the oscillator system is physically perturbed (as opposed to restarted with new initial conditions) is done in Sec.~IV.

\begin{figure}
\includegraphics[width=8cm]{Images/GCKuramotoTransientPrediction.jpeg}
\vspace*{-.5in}
\caption{GC network prediction of a 5 node Kuramoto oscillator network with a coupling strength of 5 when length of time series data analyzed is stopped at 5, 10, 15, 20 or 25 seconds. Each of the colored lines is a trajectory of one of the five oscillators. Note that the network inferences are best when the length of the simulation passed into GC is close to the transient time length of 5 seconds.}
\label{fig:kuramotoGCTransients}
\end{figure}

\begin{figure}
\includegraphics[width=8cm]{Images/GCKuramotoEndtime.jpeg}
\caption{\textbf{(a)} Plot of transient time length of Kuramoto network vs coupling strength and network size. The transient time length of the system depends inversely on the coupling strength, and this dependence stays approximately the same for all networks larger than 5 nodes. Right graph plots the transient time length vs the inverse coupling strength of Kuramoto networks ranging from 1 to 20 nodes in size. A cluster of linear regions can be seen in the graph with a slope of approximately 20 and this cluster corresponds to all networks 5 nodes or larger. Therefore, for moderately large networks the transient time length is approximately 20 over the coupling strength of the system.
\textbf{(b)} GC accuracy over Kuramoto networks with sizes $n = 5, 10, 15$ and varying simulation endtimes and coupling strengths. Black line in heatmap plots corresponds to $T(n, K) = \frac{4.5n}{K}$ (optimal simulation time length is equal to 4.5 times the network size over the coupling strength). Note that for all network sizes, GC achieves a maximum accuracy around $T(n, K)$.}
\label{fig:kuramotoGCTransientsStrength}
\end{figure}
\newpage

\subsection{Perturbation Cascade Inference}

The PCI method developed here physically forces nodes in the network to infer its underlying structure. It relies on the observation that different nodes are activated in time as information and perturbations spread across a dynamical network.  PCI learns the length of the shortest path (e.g. distance) from every node in the network to the initially perturbed node.  Indeed, PCI learns the distances of every node in the network from each perturbed node, and applies this information to reconstruct the underlying graph connectivity. To implement this approach on the mass-spring and Kuramoto oscillator networks studied in this paper, we apply impulse forcings at particular network nodes and study the transient time dynamics of the oscillators in the short time window after the perturbation was applied. We correlate the oscillator trajectories to sort all the oscillators in order of their activations in time and use these node orderings to determine the probabilities of all edges in the network. It is important to stress that this methodology is only applicable when the coupling strengths between all pairs of neighboring nodes in the oscillator network are approximately equal. Figure \ref{fig:inferenceWithPerts} shows a diagram of this method on a four node network.

\begin{figure}
    \includegraphics[width=9cm]{Images/ToyPerturbationProblem}
    \caption{A simple example illustrates the strategy to infer the structure of a network of four coupled nodes using systematic perturbations. We make the important assumption that when a node is activated, it perturbs only its nearest neighbors at the next point in time. In the first perturbation, N1 is activated which causes N2 and N3 to become perturbed which in turn cause N4 to become perturbed. Originally, all edges in the network are given a probability of 0.5. This cascade of perturbations provides new evidence that N1 causes N2 and N3, either N2 or N3 causes N4, and N1 does not cause N4. Therefore, the edges from N1 to N2 and N3 are updated to 1, the edges from N2 and N3 to N4 are updated to 0.66, and the edge from N1 to N4 is updated to 0.05. This process is repeated for the second perturbation of N2 and the third perturbation of N3 with the adjacency matrix of edge probabilities updated each time. The final probability matrix is thresholded at 0.5 which results in the fully accurate network prediction shown in the top right.}
    \label{fig:inferenceWithPerts}
\end{figure}

To reverse-engineer the structure of a network, we must understand how to probabilistically infer its edges from information about the distance of every node from a set of perturbed nodes. Let's assume that after perturbing node $p$ in the network $\mathbf{A}$ we infer the length of the shortest directed path from $p$ to any node $v \in [n]$, denoted by $\dist(p, v)$. By convention, if there is no directed path from $p$ to $v$ then we write $\dist(p, v) = \infty$. Now we define the set of all nodes which have a directed shortest path of length $k$ from node $p$
\begin{equation}
    D_k(p) = \{v: \dist(p, v) = k\}.
\end{equation}
By this definition, $D_1(p)$ is the set of all the nearest neighbors of $p$ (i.e. nodes linked by a directed edge from $p$) and $D_\infty(p)$ contains all the nodes which are disconnected from $p$. Notice that any node $y \in D_{k+1}(p)$ must be connected to some node $x \in D_k(p)$ by an outgoing edge from $x$ (i.e. caused by $x$).  Furthermore, any node $x \in D_\infty(p)$ cannot be caused by/linked to any node in $D_k(p)$ where $k < \infty$. Therefore, if by perturbing node $p$ we can learn $\dist(p, v)$ for every node $v$ in the network, then we discover that:
\begin{enumerate}
    \item All nodes caused by $p$ lie in $D_1(p)$.\\
    \item Every node in $D_{k+1}(p)$ is caused by a node in $D_k(p)$.\\
    \item Every node in $D_k(p)$ for $k < \infty$ does not cause any of the nodes in $D_\infty(p).$
\end{enumerate}

As we perturb more nodes, we would like to assign a probability to each edge in the network that it actually exists in the underlying graph. We start with an uninformative prior where every edge in the network $x \to y$ has a $\mathbb{P}(x \to y) = 0.5$ probability of existing. Then for every subsequent perturbation, we update the network edge probabilities as follows:

Suppose we have performed a perturbation of node $p$ in the network and have learned the sets $D_1(p), D_2(p), ..., D_\infty(p)$ by sorting all of the oscillators in order of their activation times. Now we take any nodes $x, y$ from our network and update the probability that the edge $x \to y$ exists as follows:\\

\begin{algorithm}[t]
\SetAlgoLined
\KwResult{Returns edge probability matrix A}
 \tcp{A(x, y) probability x causes y}
 A = ones(n, n)/2\;
 \BlankLine
 perturbNodes = [1, 2, ...]\;
 \For{$p \in \mathrm{perturbNodes}$}{
  perturb node $p$\;
  observe sets $D_1(p), ..., D_\infty(p)$\;
  \For{$x \in [n]$}{
   \tcp{Case 1}
   \If{$x \in D_\infty(p)$}{
    continue\;
   }
   \BlankLine
   \tcp{Case 2}
   $k = \dist(p, x)$\;
   \For{$y \notin \bigcup_{m=1}^{k+1} D_m(p) \cup D_\infty(p)$}{
    $A(x, y) \leftarrow A(x, y) / 10$\;
   }
   \BlankLine
   \tcp{Case 3}
   \For{$y \in D_{k+1}(p)$}{
    $c = \prod_{v \in D_k(p)} (1 - A(v, y))$\;
    $A(x, y) \leftarrow A(x, y) / (1 - c)$\;
   }
  }
 }
 \caption{Perturbation Cascade Inference \label{alg:one}}
\end{algorithm}

\textbf{Case 1:}
If $x \in D_\infty(p)$ then $x$ is not caused by $p$ and we have not learned any information about the edge $x \to y$ so $\mathbb{P}(x \to y)$ is not updated.\\

\textbf{Case 2:} If $x \in D_k(p)$ for some $0 \leq k < \infty$ and $y \in D_m(p)$ where $k + 1 < m \leq \infty$ then we know from the discussion above that $x$ does not cause $y$ so we penalize the prior probability $\mathbb{P}(x \to y)$ by dividing it by 10.\\

\textbf{Case 3:} If $x \in D_k(p)$ and $y \in D_{k + 1}(p)$, where $0 \leq k < \infty$, we apply Bayes's rule to update the probability $\mathbb{P}(x \to y)$. Since $y \in D_{k + 1}(p)$, we know that $y$ is caused by at least one of the nodes in $D_k(p)$ (not necessarily by $x$). Each node $v \in D_k(p)$ has \textit{a prior} probability $\mathbb{P}(v \to y)$ that it causes $y$. Using Bayes's rule, we know that the conditional probability that the edge $x \to y$ exists given the information that at least one of the edges $v \to y$ exists for $v \in D_k(p)$ is

\begin{eqnarray}\label{eq:13}
&&\mathbb{P}\Big(\!x \to y \ \!\Big|\!\ \exists v \!\in\! D_k(p) \!\text{ s.t. }\! v \to y \!\Big) \nonumber \\
&&\hspace*{.3in} = \! \frac{\mathbb{P}\Big(\exists v \!\in\! D_k(p) \!\text{ s.t. }\! v \to y \ \!\Big|\!\ x \to y\Big)\mathbb{P}(x \to y)}{\mathbb{P}\Big(\!\exists v \!\in\! D_k(p) \text{ s.t. } v \to y\!\Big)} \nonumber \\
&& \hspace*{.3in} =\! \frac{\mathbb{P}(x \to y)}{1 - \prod_{v \in D_k(p)}\! \Big(\!1 \!-\! \mathbb{P}(v \to y)\!\Big)}.
\end{eqnarray}

These computations produce a probability of a causal connection between one node and another. The entire algorithm is summarized in Algorithm~\ref{alg:one}.

In the discussion above, we did not mention how to infer the distance sets $D_k(p)$ after perturbing node $p$ in the network. The construction of these sets is an entirely data-driven problem that heavily depends on the time-dynamics of the system being studied. In the case of the simulated mass-spring and Kuramoto oscillator networks, we use correlations and variances between oscillator trajectories respectively to estimate which nodes lie in each distance set $D_k(p)$ by studying at what times they become activated.

The effectiveness of PCI is illustrated in the following section.

\section{Computational Results:  Coupled Oscillators}
\label{sec:results}

\begin{figure*}
\includegraphics[width=15cm]{Images/PertVarySizeForcingStrengths.jpeg}
\caption{Accuracy of PCI, GC, and CCM on perturbed mass-spring and Kuramoto oscillator ER networks (0.5 connection probability) with varying sizes, perturbation forcing magnitudes, and coupling strengths. For the mass-spring system, PCI achieves $\approx$80\% accuracy for all networks of at most 10 nodes and this accuracy increases slightly if the force magnitude by which the nodes are perturbed is increased from 10N to 50N. For larger network sizes, PCI requires lower coupling strengths in order to reach the same predictive accuracy. As an example, for the 20 node Erdos-Renyi networks, PCI achieves this accuracy only when all springs in the mass-spring network have spring constants equal to 0.1. For the Kuramoto oscillator system PCI requires large coupling and forcing strengths $(>20)$ in order to correctly infer networks of 2 nodes or larger. Therefore, if the coupling and forcing of the system are large, then PCI successfully infers Kuramoto oscillator networks with an accuracy of 80-90\%. GC is capable of inferring network structure of a 10 node mass-spring system with 85\% accuracy when the coupling strengths are 0.1 or smaller. Similar to PCI, its performance drops below 80\% for larger coupling strengths and network sizes and is not significantly affected by the magnitude with which the nodes are perturbed. GC consistently gets below 75\% accuracy on Kuramoto oscillator networks with 3 nodes or larger. Note that even with these perturbed systems, CCM is incapable of accurately predicting any structure and for networks of larger than 3 nodes, achieves an accuracy of close to 50\%.}
\label{fig:pertSFS}
\end{figure*}

In this section, we evalute the performance of GC, CCM, and PCI on mass-spring and the nonlinear Kuratomo oscillator systems where the oscillators are \textit{physically forced}.

Fig.~\ref{fig:pertSFS} gives an extensive comparison of the models as a function of the network size, coupling strength between nodes and number of forcings (perturbations). The networks are randomly generated Erdos-Renyi (ER) graphs with a 0.5 probability of connection. Fig.~\ref{fig:pertSFS} shows that CCM fails to accurately infer causal connections even when perturbations are allowed, much like what has been found in \cite{monster2017,monster2018}. Perturbed GC has a range of parameter space where accurate causal inference can be established, although its validity is especially limited for the nonlinear oscillators. PCI has strong potential for accurate causal inference, especially as the number of perturbations to the system is increased.

\begin{figure*}
\includegraphics[width=15cm]{Images/PCIAccVSPerts.jpeg}
\vspace*{-.5in}
\caption{The accuracy of PCI on inferring 20 node mass-spring networks, when nodes are perturbed randomly or in order of their outdegree or centrality. We perform this experiment on three different network types: Erdos-Renyi, Erdos-Renyi with cluster, and scale-free. We observe one case (namely on the clustered ER graphs), where perturbing less than half of the network nodes (8/20) results in an optimistic 85\% inference accuracy.}
\label{fig:PCIAccVSPerts}
\end{figure*}

Given the success of the PCI method, we explore its use on larger 20 node mass-spring networks as shown in Fig.~\ref{fig:PCIAccVSPerts}. In these experiments, we infer the structure of an Erdos-Renyi network with 0.5 connection probability, a 0.5 probability Erdos-Renyi network with a fully-connected cluster/clique of 8 nodes, and a scale-free network built by the Barab\'{a}si-Albert model~\cite{albert2002statistical}. For each graph, we try three different orders of node perturbations: random order, by decreasing order of outdegree, and by decreasing order of outcloseness. In a real experiment, we would never have a full understanding of the degree or closeness statistics of all nodes in the network. However, this experiment is meant to simulate the ``best" possible scenario for inference if we have prior knowledge about the importance or centrality of nodes in the graph.
We do not consider either perturbed GC or CCM as their performance does not scale to 20 node networks and does not significantly depend on the order of node perturbations.

PCI is shown to be a viable technique for extracting accurate causal relations provided enough perturbations of the system are allowed. The causal relations are nicely extracted for ER networks, and with ER networks containing clusters and scale-free structure. In particular, on the ER cluster network we observe that if we perturb in decreasing order of degree or closeness instead of randomly we can reach 85\% prediction accuracy with only 8 out of the 20 nodes perturbed. This shows that in some instances, PCI can infer certain network structures with a small set of properly placed perturbations, similar to the way that messages efficiently spread throughout a network when they are transmitted from highly-connected hubs.

\section{Conclusion}

In this paper, we showed that perturbing components of the network, or by observing many unique transient dynamics, a large number of potential networks can be disambiguated, achieving reconstruction of unique and accurate network structures. Specifically, we demonstrated that our PCI method, along with GC with perturbations and transients, is capable of inferring linear (mass-spring) and nonlinear (Kuramoto oscillator) networked dynamical systems.
Our proposed PCI method demonstrated consistently strong performance in inferring causal relations for small (2--5 node) and large (10--20 node) networks for both linear and nonlinear systems. Perturbed GC is capable of inferring smaller networks under low coupling strength regimes, while methods such as CCM did not infer the structure of any oscillator networks. Our analysis suggests that the ability to apply a large and diverse set of perturbations/actuations to the network, in either a targeted or random way, is critical for successfully and accurately determining causal relations and disambiguating between various viable networks.

Beyond data analysis and modeling, the problem of network inference asks a foundational question: How can a complex system be described by a small, interpretable set of causal relationships among its components?

In many scientific and engineering applications, studying transient dynamics and impulse responses have long provided insight on the desired causal structure.
We note, however, that to obtain reasonable network inferences and interpret them, much care must be taken to choose and observe only the important nodes of multi-component networks rather than analyzing the entire system in all of its complexity.
Importantly, there remain many open fundamental mathematical questions for future study.
These include how to optimally place activations and observations of network nodes for inference, what types of network dynamics contain unique information about their connectivity, and above all, whether it is possible to formulate a consistent and unifying theory of causality for time dependent systems.

\section*{Acknowledgements}
GS would like to acknowledge support from the Mary Gates Research Foundation and the National Science Foundation Graduate Research Fellowship under Grant No. 174530.
JNK acknowledges support from the Air Force Office of Scientific Research (AFOSR) grant FA9550-17-1-0329.
BWB acknowledges support from the Washington Research Foundation and Air Force Research Lab (AFRL) grant FA8651-16-1-0003.

All code and additional experiments can be found at \href{https://github.com/sgstepaniants/netinf}{https://github.com/sgstepaniants/netinf}.

\section*{Appendix}

The dynamical models used to evaluate causal relations are of two type:  linear and nonlinear.  The linear model is a standard mass-spring system where masses interact through springs and Hooke's law.  The nonlinear system is a set of nonlinear oscillators known as the Kuramoto system. Details are given in each section below.

\subsection{Mass-Spring System}
Mass-spring networks have numerous applications in a variety of disciplines including modeling of deformable objects in computer graphics \cite{nealen2005}, molecular dynamics in complex polymer materials \cite{zhang2014} and organ simulations for surgical procedures \cite{meier2005}. The system of ODEs that governs the movement of $n$-coupled oscillators in a directed mass-spring network is
\begin{equation}\label{eq:5}
m\frac{d^2\mathbf{x}}{dt^2} = k\mathbf{M}\mathbf{x} - c\frac{d\mathbf{x}}{dt} + \mathbf{f}(t),
\end{equation}
where we assume the mass $m$ and damping constant $c$ of every oscillator is the same and the spring constant $k$ is also the same for every spring in the network. $\mathbf{f}(t)$ is the forcing function for every oscillator. The binary matrix $\mathbf{M}$ defines the connectivity of the network. We impose fixed boundary conditions such that the first and last oscillators are attached to fixed walls by springs. For all other oscillators in the network, they cannot be directly connected to the boundary walls. Therefore, $\mathbf{M}$ has dimension $n + 1 \times n + 1$ and has the form
\begin{equation}
\mathbf{M} =
\begin{bmatrix}
    0 & 1 & 0 & \hdots & 0\\ \cline{2-4}
    1 & \multicolumn{3}{|c|}{\multirow{3}{*}{\scalebox{2}{$\mathbf{A}$}}} & \vdots\\
    0 & \multicolumn{3}{|c|}{} & 0\\
    \vdots & \multicolumn{3}{|c|}{} & 1\\ \cline{2-4}
    0 & \hdots & 0 & 1 & 0
\end{bmatrix},
\end{equation}
where the $n \times n$ binary adjacency matrix $\mathbf{A}$ represents the connectivity of all of the nodes in the network with the exception of the fixed boundary walls whose connectivity is predetermined. If $A_{ji} = 1$ then oscillator $i$ is connected to oscillator $j$ by a `directed spring' (i.e. oscillator $i$ can force oscillator $j$) and if $A_{ji} = 0$ then oscillator $i$ has no direct influence on oscillator $j$. Since causality is a directed relationship (i.e. $X$ causes $Y$ but $Y$ might not cause $X$), we do note impose restrictions that $\mathbf{A}$ must be symmetric.

\subsection{Kuramoto System}
The Kuramoto model proposed by Yoshiki Kuramoto is one of the most well-studied systems of nonlinear coupled oscillators \cite{kuramoto1975, kuramoto2003, strogatz2000kuramoto}. It is a canonical system for studying quasiperiodic dynamics, synchronization, and chaos and has found practical applications in a variety of areas in physics \cite{strogatz1998, strogatz1996, kourtchatov, jiang}, biology \cite{liu1997, buck1988, walker1969}, and medical sciences \cite{peskin1975, michaels1987}. Equation (\ref{eq:7}) describes the dynamics of $n$ coupled Kuramoto oscillators with a forcing term.
\begin{equation}\label{eq:7}
\frac{d\theta_i}{dt} = \omega_i + \frac{K}{n}\sum_{j=1}^n A_{ij}\sin(\theta_j - \theta_i) + f_i(t), \quad 1 \leq i \leq n.
\end{equation}

The dynamics of the $i$th oscillator is described by its angle $\theta_i$ which has a natural frequency $\omega_i$. Each $A_{ij}$ is an entry of the $n \times n$ binary adjacency matrix $\mathbf{A}$ which represents the connectivity of all $n$ oscillators in the network. All oscillators in the network are coupled to their adjacent neighbors by the same coupling strength $K$. For any Kuramoto network there exists a bifurcation value $K_c$ where for all $0 \leq K < K_c$ the oscillator trajectories are unsynchronized and for $K > K_c$ clusters of oscillators synchronize and eventually phase lock. This bifurcation value depends on various properties including the topology of the network and the distribution from which the natural frequencies $\omega_i$ are sampled from. Also, the Kuramoto model exhibits chaotic dynamics for networks of four or more nodes \cite{maistrenko2005}. Therefore, in all our experiments we test a sufficiently wide range of coupling strengths and network sizes to analyze how network inferences vary in synchronized, unsynchronized, chaotic, and non-chaotic regimes.

\bibliography{references.bib}

\end{document}

config: arxiv


TITLE: Are there integer solutions such that the sum of the roots of x and y equal the root of a prime p?
QUESTION [2 upvotes]: How could one prove whether or not integer pairs (x, y) exist such that $\sqrt{x} + \sqrt{y} = \sqrt{p}$, where p is prime? `

REPLY [0 votes]: A slightly less elementary approach. Consider the polynomial whose roots are the pairwise sums of the roots of the minimal polynomials of $\sqrt x$ and $\sqrt y$: $$(z - \sqrt x - \sqrt y)(z - \sqrt x + \sqrt y)(z + \sqrt x - \sqrt y)(z + \sqrt x + \sqrt y) \\
= ((z - \sqrt x)^2 - y) ((z + \sqrt x)^2 - y) \\
= (z^2 + x - y - 2z\sqrt x) (z^2 + x - y + 2z\sqrt x) \\
= (z^2 + x - y)^2 - 4z^2x \\
= z^4 - 2(x+y) z^2 + (x-y)^2$$
If $\sqrt p$ is a root then its minimal polynomial $z^2 - p$ is a factor. Since we have two monic polynomials, that implies that $p$ is a factor of $(x-y)^2$. The tail of Elliot's argument wraps it up.
We could instead have gone on to do a long division, finding that the remainder is $$p^2 - 2(x+y)p + (x-y) ^2 =0$$ as in Fimpellizieri's answer.


TITLE: Can I use the Polynomial Interpolation Error Theorem to conclude that $p_n(x)$ tends to $f(x)$ as $n$ goes to $\infty$?
QUESTION [2 upvotes]: I was doing some revision questions on the polynomial interpolation error theorem, $$f(x)-p_n(x) = \frac{f^{n+1}\pi_{n+1}(x)}{(n+1)!}$$  After doing out the proof of this theorem, I'm asked if you can use this theorem to conclude that $p_n(x)$ tends to $f(x)$ as $n$ goes to $\infty$?
My initial thoughts are that it can be used as $(n+1)!$ will go to infinity meaning $\frac{f^{n+1}\pi_{n+1}(x)}{(n+1)!}$ will go to zero but I'm not sure I'm right. Anyone able to tell me if I'm wrong?

REPLY [0 votes]: You cannot conclude that, the sequence of Lagrange polynomials of $f$ does not always converge pointwise towards $f$, the reason being that $\|f^{n+1}\|_{\infty}$ may explode faster that $(n+1)!$. This is known as the Runge's phenomenon.
By the way, the error term at $x$ of $f\colon]a,b[\rightarrow\mathbb{R}$ is given by:
$$\frac{f^{n+1}(\xi_x)\pi_{n+1}(x)}{(n+1)!},$$
where $\xi_x\in]a,b[$, you forgot to write $\xi_x$.


TITLE: The sixth power integral moment of automorphic L-function attached to Maass Forms
QUESTION [5 upvotes]: It is known that the sixth power integral moment of  automorphic L-function attached to Cusp Forms has been proved by M. Jutila, that is $\int_{0}^{T}|L(1/2+it,f)|^{6}dt \ll T^{2+\varepsilon}$.
 And it should also be true for Maass Forms philosophically, are there some references about this?

REPLY [3 votes]: I have seen several papers in this direction, but don't know exactly what is known. Perhaps the following reference is helpful:
M. P. Young obtains in his article The second moment of $GL(3)\times GL(2)$ L-functions at special points in $2009$ such estimates for a family of $L$-functions $L(φ×u_j,s)$ where $u_j$ runs over the family of Hecke-Maass cusp forms on $SL(2,\mathbb{Z})$, $\phi$ is a fixed $SL(3,\mathbb{Z})$ Maass form (Rankin-Selberg convolution of $\phi$ with $u_j$)
; at the special points $\frac{1}{2}+it_j$. The results are stated in Theorems $3.1$ and $3.2$. In particular,
$$
\sum_{t_j\le T}|L(u_j, \frac{1}{2}+it_j|^6 \ll T^{2+\varepsilon}.
$$


TITLE: Finding the $ n $-th derivative of $ {\cos^{n}}(x) $.
QUESTION [2 upvotes]: Problem. Find the $ n $-th derivative of $ {\cos^{n}}(x) $.

What I’m doing is substituting $ t = \cos(x) $ and then finding the $ n $-th derivative of the new function, but I’ve a feeling that this is wrong. Could anyone please point out the correct method? Thank you!

REPLY [3 votes]: Doing it directly is likely to lead to a messy computation.
From trigonometry, we know that
\begin{align*}
2^{n-1}\cos^n x = \cos nx +\binom{n}{1} \cos(n-2)x + \binom{n}{2}\cos(n-4)x + \cdots
\end{align*}
Now, we calculate $\dfrac{d^n}{dx^n}(\cos kx)$ for an arbitrary $k$ as follows:
\begin{align*}
\frac{d}{dx}(\cos kx) &= -k\sin kx = k \cos\left(\frac{\pi}{2}+kx\right) \\
\frac{d^2}{dx^2}(\cos kx) &= \frac{d}{dx}\left( k \cos\left(\frac{\pi}{2}+kx\right)\right) \\
&= -k^2 \sin\left(\frac{\pi}{2}+kx\right) \\
&= k^2 \cos\left(2\cdot\frac{\pi}{2}+kx\right)
\end{align*}
Hence, more generally,
\begin{align*}
\frac{d^n}{dx^n}(\cos kx) &=k^n \cos\left(n\cdot\frac{\pi}{2}+kx\right)
\end{align*}
Thus,
\begin{align*}
2^{n-1}\frac{d^n}{dx^n}(\cos^n x) &= n^n \cos\left(n \cdot \frac{\pi}{2} + nx\right)+\binom{n}{1}(n-2)^n \cos\left(n \cdot \frac{\pi}{2} +(n-2)x\right)+\binom{n}{2}(n-4)^n \cos\left(n \cdot \frac{\pi}{2} +(n-4)x\right)\cdots
\end{align*}


TITLE: Finding mass of a region given density
QUESTION [1 upvotes]: I am trying to solve the following exercise:

Consider the region $\Omega$, in the first quadrant, bounded by the curves:
\begin{align}
x^2 + y^2 = 4; \;\;\; x^2 + y^2 = 9; \;\;\; x^2 - y^2 = 1; \;\;\; x^2 - y^2 = 4
\end{align}
By a change of variables calculate de mass of $\Omega$ if the density at each point $(x,y)\in \Omega$ is equal to the product of the distances of the point to the coordinate axis.

My attempt is the following. First, I define de density function $\rho\colon \Omega \rightarrow \mathbb{R}$ by
\begin{align*}
\rho(x,y) = x\cdot y.
\end{align*}
Then, I applied the substitution $x = r\cdot \cos(\theta)$ and $y = r \cdot \sin(\theta)$. For the first two equations I get $r = 2$ and $r=3$. My problem lies in the equations $x^2 - y^2 = 1$ and $x^2 - y^2 = 4$, I don't know how to continue from here. Of course, what I want is to find the limits of the integral
\begin{align*}
m(\Omega) = \iiint_{\Omega} x\cdot y.
\end{align*}

REPLY [1 votes]: If $\require{cancel}u=x^2+y^2$ and $v=x^2-y^2$, then $4\leqslant u\leqslant9$ and $1\leqslant v\leqslant4$. Besides,$$xy=\sqrt{x^2y^2}=\sqrt{\frac{u+v}2\times\frac{u-v}2}=\frac{\sqrt{u^2-v^2}}2.$$Finally, $(x,y)=\left(\sqrt{\frac{u+v}2},\sqrt{\frac{u-v}2}\right)$. But if$$g(u,v)=\left(\sqrt{\frac{u+v}2},\sqrt{\frac{u-v}2}\right),$$then$$\bigl|\det Jg(u,v)\bigr|=\frac1{4\sqrt{u^2-v^2}}.$$So, since $\Omega=g([4,9]\times[1,4])$, you have\begin{align}\iint_\Omega xy\,\mathrm dx\,\mathrm dy&=\iint_{g([4,9]\times[1,4])}xy\,\mathrm dx\,\mathrm dy\\&=\iint_{[4,9]\times[1,4]}\frac{\sqrt{u^2-v^2}}2\left|\det Jg(u,v)\right|\,\mathrm du\,\mathrm dv\\&=\int_4^9\int_1^4\frac{\cancel{\sqrt{u^2-v^2}}}2\frac1{4\cancel{\sqrt{u^2-v^2}}}\,\mathrm dv\,\mathrm du\\&=\frac{15}8.\end{align}


text:
\section{Introduction}
\label{sec:introduction}

In this paper, we will consider the safety of data-trained memoryless feedback
controllers\footnote{RL-trained feed-forward neural networks, for example.} in
the following context: we will assume an autonomous vehicle that is described
by the Kinematic Bicycle Model (KBM) dynamics (a good approximation for
four-wheeled vehicles \cite{KongKinematicDynamicVehicle2015}), and a safety
property that the autonomous vehicle should avoid a stationary, fixed-radius
disk in the plane. In particular, we propose ShieldNN, an algorithm to design
Rectified Linear Unit (ReLU) safety networks for this scenario: a ShieldNN
network composed in series with \emph{any} memoryless feedback
controller\textsuperscript{1} makes the composition of the two controllers
provably safe by the aforementioned criterion. This structure itself
distinguishes ShieldNN from most other work on safe data-trained controllers:
instead of designing a single safe controller, ShieldNN uses the KBM dynamics
to design a \emph{controller-agnostic} NN that corrects -- in real-time --
unsafe control actions generated by \emph{any} controller. In other words,
ShieldNN designs a ``safety-filter'' NN that takes as input the instantaneous
control action generated by a controller (along with the state of the system)
and outputs a safe control for the KBM dynamics; this safety-filter NN thus
replaces \emph{unsafe} controls generated by the original controller with safe
controls, whereas \emph{safe} controls generated by the original controller are
passed through unaltered -- i.e., unsafe controls are ``filtered'' out. A block
diagram of this structure is depicted in \cref{fig:shieldnn_diagram}.

The benefits of this approach are manifest, especially for data-trained
controllers. On the one hand, existing controllers that have been designed
\emph{without} safety in mind can be made safe by merely incorporating the
safety filter in the control loop. In this scenario, the safety filter can also
be seen as a countervailing factor to controllers trained to mimic experts: the
expert learning can be seen as a design for ``performance'', and the safety
filter is added to correct unanticipated unsafe control behavior as needed.
\begin{wrapfigure}{r}{0.49\textwidth}
	\centering \includegraphics[width=0.49\textwidth]{shieldnn.pdf}
	\caption{Block diagram of ShieldNN in the control loop for a four-wheeled
	autonomous vehicle.}
	\label{fig:shieldnn_diagram}
	\vspace{-12pt}
\end{wrapfigure}
On the other hand, the controller-agnostic nature of our proposed filter means
that ShieldNN itself may be incorporated into training. In this way, the safety
filter can be seen to function as a kind of ``safety expert'' during training,
and this can potentially improve sample efficiency by eliminating training runs
that end in unsafe states.

The main theoretical contribution of this paper is thus the development of the
ShieldNN algorithm. The central pillar of ShieldNN is the notion of a barrier
function, because a barrier function allows the safety problem for a
\emph{feedback} controller to be recast as a set membership problem for the
\emph{outputs} of said controller. In particular, this recasting reduces the
safety-filter design problem into one of designing a \emph{prediction-style} NN
whose outputs are constrained to lie in a specific set. As a prerequisite for
ShieldNN then, we propose a novel class of candidate barrier functions for the
KBM dynamics that is characterized by three real-valued parameters (one of
which is the safety radius). The necessity for a \emph{class} of candidate
barrier functions stems from the difficulty in analytically designing a barrier
function for the KBM dynamics.

Thus, ShieldNN is functionally divided into two parts: a verifier, which
soundly verifies a particular choice of barrier function (from the
aforementioned class), and a synthesizer, which designs the actual ShieldNN
filter.

Furthermore, we also validated these theoretical results with experimental
validation on four-wheeled vehicle models. In particular, we apply ShieldNN
safety-filters both before and after RL training for an autonomous vehicle
simulated in CARLA \cite{Dosovitskiy17}. Our results show that incorporating
ShieldNN dramatically improved the safety of the vehicle: it reduced the number
of obstacle collisions by 99.4\%-100\% in our safety experiments. We also
studied the effect of incorporating ShieldNN during training: for a constant
number of episodes, 28\% less reward was observed when ShieldNN wasn’t used
during training. This suggests that ShieldNN has the further property of
improving sample efficiency during RL training.

\noindent \textbf{Related Work.}
Motivated by the lack of safety guarantees in Deep RL, recent works in the literature of safe RL have focused on
designing new RL algorithms that can take safety into account. The work in this area can be classified into three
categories. The works in the first category focus on how to modify the training algorithm to take into account safety constraints. Representative examples of this work include reward-shaping~\cite{saunders2018trial}, Bayesian and robust regression~\cite{liu2019robust,berkenkamp2016bayesian,pauli2020training}, and policy optimization with constraints~\cite{gaskett2003reinforcement,moldovan2012safe,turchetta2016safe,wen2020safe}. Unfortunately, such approaches do not
provide provable guarantees on the safety of the trained controller. The second category of literature focuses on using
ideas from control theory to augment the RL agent and provide safety guarantees. Examples of this literature include the use of Lyapunov methods~\cite{berkenkamp2017safe,chow2019lyapunov,chow2018lyapunov}, safe model predictive control~\cite{koller2018learning}, reachability analysis~\cite{akametalu2014reachability,govindarajan2017data,fisac2018general}, barrier certificates~\cite{wang2018safe,cheng2019end,wabersich2018scalable,srinivasan2020synthesis,taylor2020control,li2019temporal,cheng2020safe,wang2018safe,robey2020learning}, and online learning of uncertainties~\cite{shi2019neural}. Unfortunately, such methods suffer from either being computationally expensive, specific to certain controller structures or training algorithms or require certain assumptions on the system model. Finally, the third category focus on applying formal verification techniques (e.g., model checking) to verify formal safety properties of pretrained RL agents. Representative examples of this category are the use of SMT-like  solvers~\cite{sun2019formal,dutta2018output,liu2019algorithms} and hybrid-system verification~\cite{fazlyab2019efficient,xiang2019reachable,ivanov2019verisig}. These techniques only assess the safety of a given RL agent instead of designing a safe agent.