# IGN : Implicit Generative Networks

Haozheng Luo, Tianyi Wu, Colin Feiyu Han, Zhijun Yan  Computer Science
Northwestern University
{robinluo2022,tianyiwu2023,colinfhan<EMAIL_ADDRESS>
###### Abstract

In this work, we build recent advances in distributional reinforcement
learning to give a state-of-art distributional variant of the model based on
the IQN. We achieve this by using the GAN model’s generator and discriminator
function with the quantile regression to approximate the full quantile value
for the state-action return distribution. We demonstrate improved performance
on our baseline dataset - 57 Atari 2600 games in the ALE. Also, we use our
algorithm to show the state-of-art training performance of risk-sensitive
policies in Atari games with the policy optimization and evaluation.

###### Index Terms:

Reinforcement Learning, Implicit Quantile Networks, Generative adversarial
network

##  I Introduction

Distributional reinforcement learning [1] focus on the reward with intrinsic
randomness within the reinforcement learning framework. Distributional
reinforcement learning aims to model the distribution over returns and to use
these distributions to evaluate and optimize a policy. Any distributional
reinforcement learning algorithm is characterized by two aspects: the
parameter of the return distribution and the distance metric or loss function
being optimized. Distributional reinforcement learning with Quantile
regression [2] gives a considerable direction for Distributional reinforcement
learning. The work based on the Deep Q-network [3] became a popular direction
in the reinforcement learning research work, especially in the game theory. It
helps the distribution over returns is modelled explicitly instead of only
estimating the mean. In recent years, a lot of work shows a good improvement
in the game theory with the reinforcement learning policy, such as Implicit
Quantile Networks (IQN) [4] and Fully Parameterized Quantile Function (FQF)
[5]. They give a great improvement on the policy, especially risk-sensitive
policy. However, those works still show limited work. How to use the
Generative Network to help improve the performance of the Distributional
reinforcement learning has become a controversial problem. Because our
training resource is 57 Atari 2600 games in the ALE, which is high-quality
images, it gives us an idea using the transfer learning to improve the
performance of the algorithm in our work. We introduced a new model called
Implicit Generative Networks (IGN) based on the IQN. The Generative Network
can work by sending the embeddings of the source and target task to the
discriminator and using the generator to generate the learning thread of the
model. The resulting loss is then (inversely) backpropagated through the
encoder. This method is also regarded as one kind of transfer learning [6]
which can significantly reduce training time. In our work, we freeze both the
critic and generator of the original GAN’s low-level layers. Also, we designed
a method that gives out an auto-encoder constraint in order to compatible the
internal representations of the critic and the generator. This assumption
explains the accelerate in training time. Also, the generator can provide us
with a clear boundary of coverage to help the model acceleration the training.
In our experiment, we will use 57 Atari 2600 games in the ALE as the dataset
to compare with our baselines, such as IQN and FQF, to show start-of-art
performance on policy. The main contributions are as follows:

  * •

Using the generator and discriminator to rewrite the quantile function in the
IQN model

  * •

Build our new algorithm IGN with the base of the IQN model

  * •

Show our algorithm state-of-art performance with the baselines in the Atari
dataset.

##  II Related Work

###  II-A Dataset

There are a lot of datasets we can use to show our models’ performance and
compare them. In our work, we want to compare our model with baselines with
several classic datasets. One of them is Atari games datasets (see Fig. 1)
[7], which is one of the classic datasets in reinforcement learning in game
theory. We use the 57 Atari 2600 games in the ALE, and it can help show our
contribution to our algorithm in the risk-sensitive policies, in which the
Atari 2600 games have 200k frames collected using the OpenAI baseline and
implemented based on Advantage Actor Critic (A2C). Also, they use sticky
actions to make the problem become more challenging. There is a $25\%$
probability that the agent’s previous action is executed. It will greatly help
us show our algorithm performance on risk-sensitive policy.

Figure 1: Atari Game

###  II-B Quantile Networks and Transfer Learning

Before our work, there was a lot of research already work in reinforcement
learning. Implicit Quantile Networks (IQN) [4] is one of the models that shows
a good result. It brings a huge improvement on the risk-sensitive policy
training. However, the model based on the Quantile distribution make a huge
fragile on the performance of the model, and it should cost huge time for
training the model. We expected to design a state-of-art algorithm that has
more robustness and has more effectiveness in training the policy. One of the
possible methods is using transfer learning. Transfer learning [8] is an
optimization method that can improve the performance of the model or rapid
progress in the second task. In our research, transfer learning can help us
have a higher start and higher slope from the second training. Also, it will
have great help in the policy evaluation. There are some state-of-art transfer
learning methods, such as VGG [9], ResNet [10] Model. From several model
pooling, there is one algorithm using Gradient method is GAN model, which can
help our model have a more competitive algorithm with origin IQN model

###  II-C Existing Offline Policy Evaluation Methods

Most existing offline policy evaluation methods under the above setting with
mean criteria can be grouped into three categories. The first category is
direct methods, which directly estimate the state-action value function [11,
12], also known as the Q-function. The second category of offline policy
evaluation methods is motivated by marginal importance sampling [13]. The last
category of methods combines direct and marginal importance sampling methods
to construct some doubly robust estimator, which is also motivated by the
efficient influence function like equation 1. [14]. In which, $Q^{\pi}$ and
$\omega^{\pi}$ are two nuisance functions.

$\frac{1}{T}\sum_{t=0}^{T-1}\omega(S_{t},A_{t})(R_{t}+\\\
\gamma\sum_{a^{\prime}\in\mathcal{A}}\pi(a^{\prime}|S_{t+1})Q^{\pi}(S_{t+1},a^{\prime})-Q^{\pi}(S_{t},A_{t}))\\\
+(1-\gamma)\mathbb{E}_{S_{0}\sim\mathbb{G}}[\sum_{a_{0}\in\mathcal{A}}\pi(a_{0}|S_{0})Q^{\pi}(S_{0},a_{0})]-\mathcal{V}(\pi)$
(1)

$\mathcal{V}(\pi)=(1-\gamma)\int_{s\in\mathcal{S}}\sum_{a\in\mathcal{A}}\pi(a|s)Q^{pi}(s,a)\mathbb{G}(ds)$
(2)

###  II-D Gradient Method

In our work, we hope to use the generative method to improve our model. The
most popular method on the generative method is Generative Adversarial
Networks (GAN) [15]. The GAN model shows a great improvement for our model in
the gradient periods and make our model have high efficiency to get the best
reward score.

GANs can train two competing models (typically NNs)[15]. The generator takes
noise $z\sim P_{z}$ as input and generates samples according to some
transformation, $G_{\theta}(z)$. The discriminator takes samples from both the
generator output and the training set as input and aims to distinguish between
the input sources. Goodfellow[15] measured the discrepancy between the
generated and the real distribution using the Kullback-Leibler divergence.
However, this approach was improved by Arjovsky[16] by using the Wasserstein-1
distance [17]. WassersteinGANs exploit the Kantorovich-Rubinstein duality[18],

$W_{1}(\mathbb{P}_{r},\mathbb{P}_{g})=\sup_{f\in
1-Lip}\\{\mathbb{E}_{x\sim\mathbb{P}_{r}}f(x)-\mathbb{E}_{x\sim\mathbb{P}_{g}}f(x)\\}$
(3)

where $1-\text{Lip}$ is the class of Lipschitz functions with Lipschitz
constant 1, in order to approximate the distance between the real
distribution, $\mathbb{P}_{r}$, and the generated one, $\mathbb{P}_{g}$. The
GAN objective is then to train a generator model $G_{\theta}(z)$ with noise
distribution $P_{z}$ at its input, and a critic $f\sim 1-\text{Lip}$,
achieving

$\min_{G_{\theta}}\max_{f\sim
1-{Lip}}\\{\mathbb{E}_{z\sim\mathbb{P}_{r}}f(x)-\mathbb{E}_{z\sim\mathbb{P}_{z}}f(G_{\theta}(z))\\}$
(4)

The WGAN-GP[16] (Fig. 2(b)) is a stable training algorithm for equation 4 [17]
that employs stochastic gradient descent and a penalty on the norm of the
gradient of the critic with respect to its input.

Figure 2: GAN configurations. [17] (a) Bellman-GAN; (b) WGAN [16]

##  III Method

###  III-A Task Definition

The problem is driven by recent development in distributional reinforcement
learning. [4] Consider a single trajectory ${(S_{t},A_{t},R_{t})}_{t}\geq 0$
where $(S_{t},A_{t},R_{t})$ denotes the state-action-reward triplet collected
at time t. We use $S\subseteq\mathbb{R}_{p}$ and $\mathcal{A}$ denote the
state and action space, respectively. We assume A is discrete and finite, and
rewards $R_{t}$ are uniformly bounded by $R_{max}$. Consider the offline
setting. Then we observed data consisting of N trajectories, corresponding to
N independent and identically distributed copies of
${(S_{t},A_{t},R_{t})}_{t}\geq 0$. For any i = $1,\cdots,N,$ data collected
from the $i^{th}$ trajectory can be summarized by
${(S_{i,t},A_{i,t},R_{i,t},S_{i,t+1})}0\leq t<T_{i}$ , where $T_{i}$ denotes
the termination time. For simplicity, we assume $T_{i}$ = T for i =
$1,\cdots,N.$

A policy defines the agent’s way of choosing the action at each decision time.
We focus on the stationary policy. Specifically, at each time point t, a
stationary policy $\pi$ maps the current state value into a probability mass
function over the action space, i.e., $\pi(a|s)$ denotes the probability of
choosing action a given the state value s.

In our work, we will work on the policy evaluation and policy optimization. In
policy evaluation, we want to see the performance of the model what would have
happened in the absence of a policy. We will use the Pong environment agent as
a fixed agent, which we will use as one of the offline policy, and we will use
our model to train the online policy and see the performance of the model
training. In our work, we will compare the reward scores and efficiency of the
timestamps because comparing with the complex model with a simple model in the
resource using is not a fit game, we will compare the timestamps( measured
with the training cycle time, not the real time). In our future work, if there
is one model that requires less real time is one of the potential research
directions in our work. For the policy optimization, we will train a model
with different environments to evaluate our model performance in timestamps
and reward gain with other baselines models.

###  III-B Distributional Reinforcement Learning

We introduce a state that uses quantile regression to approximate risk-
sensitive strategies, in other words, an optimization method for reinforcement
learning via a full quantile function of the behavioural reward distribution.

We always model the agent and the environment as a Markov decision model[19]
in classical reinforcement learning designs. Here, we define the rule
$\pi(\cdot|x)$ as a mapping from states to distributions constantly changing
with actions. For any agent that follows the rule $\pi$, we define the
discounted sum of future rewards as a function $Z$. Clearly, $Z^{\pi}(x,i)$
varies with state and action as following

$Z^{\Pi}(x,i)=\sum\limits_{t=0}^{\infty}\beta^{t}R(x_{t},a_{t})$ (5)

, where $\beta$ is a discount factor and $\beta\in(0,1)$. The action-value
function can be defined as $Q^{\pi}(x,i)=\mathbb{E}[Z^{\pi}(x,i)]$, in which
$\mathbb{E}$ is the Bellman equation operator[20].

$Q^{\Pi}(x,i)=\mathbb{E}[R(x,i)]+\beta\mathbb{E}_{P,\Pi}[Q^{\Pi}(x^{\prime},i^{\prime})]$
(6)

Azer’s research[21] on reinforcement learning proposed a distributed
reinforcement learning theory that values the benefit distribution
$Z^{\pi}(x,i)$ rather than the action-value function $Q^{\pi}(x,i)$. According
to the analysis[22], it can be found that distributed reinforcement learning
methods have substantial advantages in terms of sample complexity, final
performance, and hyperparameter processing robustness.

###  III-C Quantile Regression Approaching

According to Bellemare’s study[23], the distributed Bellman operator is a
contraction in the p-Wasserstein metric. There are two main theories for
minimizing the Wasserstein metric. One is distributed reinforcement learning
using quantile regression and shows that the resulting predictive distributed
Bellman operator is a contraction of the $\infty$-Wasserstein metric by
properly selecting the quantile target[24]. Another study shows that the
original class of classification algorithms is a contraction in the Cramér
distance, the L2 metric on the cumulative distribution function.[25]

QR-DQN minimizes the Wasserstein distance distribution Bellman objective by
estimating the quantile function at precisely chosen points. This estimation
uses quantile regression, which has been shown to converge to the true
quantile function value when minimized using stochastic approximation. In QR-
DQN, the stochastic returns are approximated by a homogeneous mixture of N
Diracs,

$Z_{\theta}(x,i):=\frac{1}{N}\sum\limits^{N}_{i=1}a_{\theta_{i}(x,i)}$ (7)

in which, $\theta_{i}$ is assigned a fixed quantile target $\tau$. We can use
Huber quantile regression loss to train this quantile target $\tau$ with
threshold $\delta$.

$\rho^{\delta}_{\tau}(a_{ij})=|\tau-\mathbb{I}\\{a_{ij}<0\\}|\frac{L_{\delta}(a_{ij})}{\delta}\text{,
with}$ (8) $L_{\delta}(a_{ij})=\left\\{\begin{aligned}
\frac{1}{2}a_{ij}^{2}&&|a_{ij}|\leq\delta\\\
\delta(|a_{ij}|-\frac{1}{2}\delta)&&|a_{ij}|>\delta\end{aligned}\right.$ (9)

In practice, we use the IQN as the baseline and implement the gradient method
with the GAN model and calculate the GAN score as the gradient score and
calculate the new gradient with the Huber loss (see equation 9) to get the
reward score and use this new reward score into the agent.

$f^{\pi}_{Y}\left(y,s,a\right)=\mathbb{E}\left[\sum\limits_{a^{\prime}\in\mathcal{A}}\pi(a^{\prime}|S_{t+1})f^{\pi}_{Y}(\frac{y-R_{t}}{\gamma},S_{t+1},a^{\prime})|S_{t}=s,A_{t}=a\right]$
(10)

###  III-D Distributional Bellman Equation

Define a random variable $Y=\sum\limits_{t=0}^{+\infty}\gamma^{t}R_{t}$, whose
support is
$\left[-\frac{R_{max}}{1-\gamma},\frac{R_{max}}{1-\gamma}\right]$ by
assumption. Given a target policy $\pi$, our goal is to estimate the
distribution of $Y$ given any state $s\in\mathcal{S}$ and $a\in\mathcal{A}$.
We use $\mathcal{U}$ to denote the product space of cumulative rewards, state
and actions, i.e.,
$\mathcal{U}=\left[-\frac{R_{max}}{1-\gamma},\frac{R_{max}}{1-\gamma}\right]\times\mathcal{S}\times\mathcal{A}$.
For any $y\in\mathbb{R}$, define either a conditional density or a conditional
probability mass function of $Y=y$ given a state-action pair $(s,a)$ as

$f^{\pi}_{Y}\left(y|S_{0}=s,A_{0}=a\right)\triangleq
f^{\pi}_{Y}\left(y,s,a\right),$

where $f^{\pi}_{Y}$ indicates actions are selected according to $\pi$. In
order to estimate $f^{\pi}_{Y}$ given any $(s,a)$, we have the distributional
Bellman equation that under Assumptions (A1) and (A2), for any $t\geq 0$,
$y\in\mathbb{R}$, $s\in\mathcal{S}$ and $a\in\mathcal{A}$, the following
equation holds.

Motivated by the distributional Bellman operator is $\pi$-contraction under a
modified Wasserstein distance, and the recently developed generative
adversarial networks (GANs), we propose a Wasserstein GAN-based distributional
policy evaluation method to estimate $f^{\pi}_{Y}$ will discuss the benefits
of using GANs. For any continuous function $h$ defined over
$\mathbb{R}\times\mathcal{S}\times\mathcal{A}$, the following equation holds.

$\displaystyle\mathbb{E}\left[\frac{1}{T}\sum\limits_{t=0}^{T-1}\mathbb{E}^{\pi}\left\\{h\left(\sum\limits_{t^{\prime}=t}^{+\infty}\gamma^{t^{\prime}-t}R_{t^{\prime}},S_{t},A_{t}\right)|S_{t},A_{t}\right\\}\right]=\gamma\mathbb{E}\left[\frac{1}{T}\sum\limits_{t=0}^{T-1}\mathbb{E}^{\pi}\left\\{h\left(\gamma\sum\limits_{t^{\prime}=t+1}^{+\infty}\gamma^{t^{\prime}-t-1}R_{t^{\prime}}+R_{t},S_{t},A_{t}\right)|S_{t+1},R_{t},S_{t},A_{t}\right\\}\right]$
(11)

Motivated by equation 11, we propose to use generative adversarial networks to
estimate $f^{\pi}_{Y}$ of $Y$ given state-action pair
$(s,a)\in\mathcal{S}\times\mathcal{A}$. Specifically, we aim to learn a
generator $\mathbb{G}^{\pi}$ that takes $(s,a)$ a uniform random variable $Z$,
i.e.,$Z\sim$ uniform (0, 1), as the input, and the output is a pseduo sample
$\tilde{Y}$ as one realization of target conditional distribution
$f^{\pi}_{Y}$. This is also inspired by the conditional GANs introduced by
Mirza and Osindero[26]. In typical GANs, the generator is trained by
minimizing some divergence between the conditional distribution of $Y|s,a$ and
$\tilde{Y}|s,a$. Alternatively, a generator $\mathbb{G}^{\pi}(Z,S,A)$ is
trained such that the joint distribution of $(Y(S,A),S,A)$ is the same as
$(\mathbb{G}^{\pi}(Z,S,A),S,A)$, where $Y(S,A)$ is a random variable that has
the same distribution of $Y$ given $s\in\mathcal{S}$ and $a\in\mathcal{A}$.
However, different from conventional distribution learning tasks, neither $Y$
nor $Y(S,A)$ is not observed in our problem. Nevertheless, the distributional
Bellman equation and the equation11 indicate us to obtain $\mathbb{G}^{\pi}$
by solving the following min-max optimization problem, which has the similar
spirit of Wasserstein-GANs[27].

$\mathop{\mathrm{minimize}}\limits_{\mathbb{G}}\mathop{\mathrm{sup}}\limits_{h:\parallel
h\parallel_{\text{Lip}}\leq 1}\\\
\left\\{\gamma\mathbb{E}\left[\frac{1}{T}\sum\limits_{t=0}^{T-1}\sum\limits_{a\in\mathcal{A}}\pi(a|S_{t+1})h(\gamma\mathbb{G}(Z_{t},S_{t+1},a)+R_{t},S_{t},A_{t})\right]\mathbb{E}\left[\frac{1}{T}\\\
\sum\limits_{t=0}^{T-1}h(\mathbb{G(Z_{\textit{t}},S_{\textit{t}},A_{\textit{t}}),S_{\textit{t}},A_{\textit{t}})}\right]\right\\}$

where $\\{Z_{t}\\}_{0\leq t\leq(T-1)}$ are set to independently follow
uniform(0,1), and $h$ is a function mapping from $\mathbb{R}^{p+2}$ to
$\mathbb{R}$.

In practice, as in conventional GANs, we restrict $h$ to some set of neural
networks $\mathcal{R}_{D}$ in addition, which can provide tractable
computation. In addition, by powerful deep neural networks, we can handle
potential high-dimensional state space $\mathcal{S}$. Similarly, we also use
some set of neural networks $\mathcal{R}_{G}$ to approximate the generator
$\mathbb{G}^{\pi}$.

###  III-E Estimation and Algorithm

Our goal is to find a generator $\mathbb{G}^{\pi}(\bullet,S,A)$ to approximate
the conditional distribution of $Y$ given $S$ and $A$. This is, in general, a
very challenging problem, especially when $S$ is a high-dimensional state,
e.g., the number of states in $\mathcal{S}$ is huge. Thus, suggested by Haas
and Richter[28], we impose some structure assumptions on the underlying true
distribution. In particular, we assume the distribution of $(Y(S,A),S,A)$ lie
in a $d_{\mathbb{G}}$-dimensional subspace with $d_{\mathbb{G}}<p+2$. Based on
this consideration, we define the following function class for the generator.

Let $\mathcal{G}(d_{\mathbb{G}},\beta,K)$ be the class of all measurable
functions $\mathbb{G}:\mathbb{R}^{p+2}\rightarrow\mathbb{R}$ such that any
component only depends on $d_{\mathbb{G}}$ arguments and lies in
$C^{\beta}(\mathcal{U},K)$.

Given the offline dataset
$\mathcal{D}_{N}=\\{(S_{i,t},A_{i,t},R_{i,t},S_{i,t+1})\\}_{0\leq
t\leq(T-1),1\leq i\leq N}$, we use the empirical average to approximate the
objective function and solve the following optimization to estimate the
generator $\mathbb{G}^{\pi}$.

$\mathop{\mathrm{minimize}}\limits_{\mathbb{G\in\mathcal{R}(L_{\mathbb{G}},\textbf{p}_{\mathbb{G}},\text{s}_{\mathbb{G}})}}\mathop{\mathrm{sup}}\limits_{h\in\mathcal{R}(L_{h},p_{h},s_{h}),\parallel
h\parallel_{\text{Lip}}\leq 1}$ (12)
$\left\\{\gamma\frac{1}{NT}\sum_{i=1}^{N}\sum_{t=0}^{T-1}\sum_{a\in\mathcal{A}}\pi\left(a|S_{i,t+1}\right)h(\gamma\mathbb{G}\left(Z_{i,t},S_{i,t+1},a\right)+R_{i,t},S_{i,t},A_{i,t})\frac{1}{NT}\sum_{i=1}^{N}\sum_{t=0}^{T-1}h\left(\mathbb{G}\left(\tilde{Z}_{i,t},S_{i,t},A_{i,t}\right),S_{i,t},A_{i,t}\right)\right\\}$
(13)

where $\left\\{Z_{i,t}\right\\}_{1\leq i\leq N;0\leq t\leq(T-1)}$ and
$\left\\{\tilde{Z}_{i,t}\right\\}_{1\leq i\leq N;0\leq t\leq(T-1)}$ are two
independent samples from uniform (0, 1). The optimal solution of the
optimization above always exists because the inner maximization problem is
Lipschitz with respect to $\mathbb{G}$ and $\mathbb{G}$ is Lipschitz with
respect to all parameters in
$\mathcal{R}(L_{\mathbb{G}},p_{\mathbb{G}},s_{\mathbb{G}})$ with bounded
support. Similarly, the inner maximization problem also has an optimal
solution. Stochastic gradient descent algorithms as the state-of-the-art can
be implemented to solve the above optimization problem. We denote the final
minimizer as by $\hat{\mathbb{G}}^{\pi}$, which can generate a sequence of
pseudo samples $\\{\tilde{Y}^{\pi}_{m}(s,a)\\}_{1\leq m\leq M}$ for some
integer $M$ for approximating the distribution of $Y$ given state-action pair
$(s,a)$ under the target policy $\pi$. This provides a flexible way for
estimating interesting quantities related to the target policy $\pi$.

In all our examples, if $\hat{\mathbb{G}}^{\pi}$ is consistent and $M$ (or
$K$) diverges to infinity, $\hat{\mathcal{V}}(\pi)$, $\hat{q}^{\pi}_{\tau}(Y)$
and $\widehat{Var}^{\pi}(Y)$ are all consistent.

Input: A generator $\hat{\mathbb{G}}^{\pi}$ and $(s,a)$.

for m = 1 to M: do

(a) Sample $\textit{Z}\sim$ uniform (0,1)

(b) Compute the pseduo sample
$\tilde{Y}_{i}(s,a)=\hat{\mathbb{G}}^{\pi}(Z,s,a)$.

Output $\\{Y_{i}(s,a)\\}_{1\leq m\leq M}$.

Algorithm 1 Generate pseudo samples to approximate $d^{\pi_{old}}$

##  IV Implementation

Consider the neural network structure used by the IQN [4]. Let $\psi$ be the
function of the functions computed by the convolutional layers and f be the
function $\mathbb{R}^{d}\to\mathbb{R}^{|A|}$ with the subsequent fully-
connected layers mapping $\psi(x)$ to the estimated action-values. Also, there
is one additional function $\phi:[0,1]\textrightarrow R^{d}$ computing an
embedding for the sample point $\tau$.For our network, we use the same
functions $\psi$, $\phi$ and f as in IQN. We make an improvement on the
Quantile loss with our new methods.

Our training algorithm in IGN after the agent generates the states and
actions, we use the original DQN model as the generator and use the three
hidden layers neural network and ReLU activation function as the discriminator
to train the model in critic with the algorithm 2 [17]. In the calculation of
the Quantile loss, we put the result of the quantile result into the
discriminator and use the equation 14 to get the GAN loss. Our IGN model
optimizes with the gradient result of the GAN loss. Finally, we will use the
Huber loss quantile result like an equation 9 to update the discriminator
parameter in the IGN models .

Furthermore, We expected that N, the number of samples in $\tau\sim$ U([0,
1]), would increase the sample complexity of IQN, with larger values leading
to faster learning. Also, when we set the $\tau$ with N = 1, we would reach
DQN performance. This would confirm the theory that many distributional RL
algorithms’ higher performance is due to their role as auxiliary loss
functions, which would vanish if N = 1 was used.

To achieve model evaluation work, we make evaluation agent into our IQN model
as the fixed policy. We use the fixed policy to get the result of the fixed
network and use the new policy training in the online policy. We use the same
quantile loss calculation function in the policy evaluation.

Also we use policy optimization and policy evaluation to show the performance
of our IGN with the baselines. We implement the data collector function with
Tensor-board logs. We also implement the Wasserstein-1 distance [17] and
discriminator loss calculation function in the visualization method. This will
be discussed more on the experiment part.

Train critic(repeat $n_{c}ritic$ times)

(a) Sample
$\\{(s_{t}^{(i)},a_{t}^{(i)},r_{t}^{(i)},s_{t+1}^{(i)},a_{t+1}^{(i)})\\}^{m}_{i=1}$
from replay pool

(b) Sample both $\\{z^{(i)}\\}$ and $\\{z^{\prime(i)}\\},i=1,\cdots,m,$ from
$P_{z}$ and $\epsilon^{(i)},i=1,\cdots,m,$ from U$[0,1]$

(c)
$x_{\theta}^{(i)}=G_{\theta}(z^{i}|s_{t}^{(i)},a_{t}^{(i)}),x_{\theta}^{\prime(i)}=r_{t}^{(i)}+\gamma
G_{\theta}(z^{\prime(i)}|s_{t+1}^{(i)},a_{t+1}^{(i)})$

(d)
$\tilde{x}^{(i)}\leftarrow\epsilon^{(i)}x_{\theta}^{(i)}+(1-\epsilon^{(i)})x_{\theta}^{\prime(i)}$

(e)
$g_{\omega}\leftarrow\frac{1}{m}\bigtriangledown_{\omega}\sum_{i=1}^{m}[f_{\omega}(x_{\theta}^{(i)})-f_{\omega}(x_{\theta}^{\prime(i)})+\lambda(\parallel\bigtriangledown_{\tilde{x}}f_{\omega}(\tilde{x}_{\theta}^{(i)})\parallel-1)^{2}]$

(f) $\omega\leftarrow Adam(\omega,g_{\omega},\alpha)$

Algorithm 2 Train critic

$g_{\theta}=-\frac{1}{m}\bigtriangledown_{\theta}\sum_{i=1}^{m}[f_{\omega}(x_{\theta}^{(i)})-f_{\omega}(x_{\theta}^{\prime(i)})]$
(14)

##  V Experiment

In order to show the state-of-art of our IGN model, we will show the
performance of policy evaluation and policy optimization in our experiment. In
our work, we compare IGN performance with the final reward and the timestamp
the final reward get. And all of the work is based on the dataset - 57 Atari
2600 games in the ALE. We focus on two main questions, as follows.

  1. 1.

Whether the fixed model in IGN will use fewer timestamps in policy evaluation

  2. 2.

Whether our IGN model will use fewer timestamps to achieve the final reward
than IQN, FQF.

Firstly, let us discuss policy optimization. We compare our method with IQN,
FQF and other baselines from the efficiency of the model and best reward
getting with the policy. And We measure Mean, Q-value and W-distance from
eight different models (IQN, IGN, Seaquest, Qbert, Pong, Breakout, MsPacman,
Freeway) and comparing with the result with our model.

Then in policy evaluation, we use several policies to train the specific
environment with the fixed policy. In our work, we want to show whether the
fixed policy can make new environment training become more efficient. We have
Q-distribution for a specific evaluated policy based on different training
steps (timestamps).

Also, we introduce the method of Monte Carlo Estimation to determine the
policy. Monte Carlo Estimation is being used when state values are not enough
since we do not have the model of the environment. The make use of this method
will estimate action values instead of state values; for here, action values
include LEFT, RIGHT, LEFTFIRE, RIGHTFIRE, NOOP and FIRE. The goal for this
method is to have an estimation of the optimal action values. The result
graphs show different actions with different scores. For comparison, Q-online
estimation is also graphed for six actions.

##  VI Result

Firstly, we do the experiment to show the performance of our model work in the
policy evaluation. The graphs on the left side will have Q-Fixed and Q-Online
on the graph, where Q-Fixed is the Q-function of the fixed policy. The graphs
on the right side will have Wasserstein-1 distance and Wasserstein-2 distance
between Q-Fixed and Q-Online. The example graph is based on IQN and IGN for
the Pong environment (Fig. 3, Fig. 4). For other fixed policies, graphs will
be in the appendix. In our experiment, It is easy to realize our IGN can
reduce the W-distance in the policy evaluation quicker than the IQN.

Figure 3: Pong environment for Mean Q-value and W-distance graph based on IQN
Figure 4: Pong environment for Mean Q-value and W-distance graph based on IGN

In order to see how the model evaluation work in our IGN model, we use the
fixed policy ( Pong environment and IGN pretrained policy) to see how the
training work in the several timestamps. In each graph, the left figure
represents the state of the game; the top-mid number represents the training
steps where k means thousands; the mid indicate current action; the right
graph will show Q-distribution of Q-Fixed and Q-online where Q-value on the
X-axis and density on the Y-axis. The example graph is being shown here for
150k training steps(Fig. 5) and 1700K training steps (Fig. 6), other results
with different training steps will be in the appendix.

Figure 5:  Q-distribution based on 150K training steps Figure 6:
Q-distribution based on 1700K training steps

Also, we use the Monte Carlo Estimation and Q-online Estimation to visualize
the estimation distribution of the reward in our policy evaluation work. The
Fig. 7 shows Monte Carlo Estimation on six actions and the Fig. 8 shows
Q-online Estimation on six actions. The x-axis of both figures shows the
reward of the actions getting in estimation. The y-axis shows the number of
estimates getting a specific reward score for actions with a unit of thousands
(k).

Figure 7: Example graph for Monte Carlo Estimation on 6 actions Figure 8:
Example graph for Q-online Estimation on 6 actions

Finally, we will use the tensor-board to show the result of the policy
optimization result and compare it with other policy methods to show our model
is states-of-art. The graphs will have the training process of IGN and IQN
with blue and red color. The left graph shows the returned test result, and
the right graph shows the returned train result. (Fig. 9). In the graph, we
found our model IGN can get the same final reward score as the IQN model but
use fewer timestamps to arrive at the highest reward score than the IQN model.
It indicates our model can have efficiency in the model training.

Figure 9: Example graph for Policy Optimization

##  VII Conclusion and Future work

In our work, we make a new model, IGN, which shows a state-of-art performance
in distributional reinforcement learning. To increase the algorithm’s data
efficiency, IGN may be trained using as little as one sample from each state-
action value distribution or as many as computational restrictions allow. In
addition, IGN enables us to broaden the scope of control rules to include a
wide range of risk-sensitive policies linked to distortion risk measurements.
Also, we show the state-of-art performance of our IGN with other baselines on
the Atari-57 benchmark. However, our model still has some limitations.
Firstly, although our model needs less timestamp in the policy optimization
and evaluation, the model still needs a lot of real-time training. Because we
use the GAN model, it requires a higher cost than the IQN, FQF model. We hope
we can find one method which has a similar or better performance than IGN but
less cost required. Also, our model still has some sensitivity with the
robustness of the noise in the training process work. We are considering using
the transformer[29] and Kernel prepossessor [30] to see whether the model
performance will become better.

## References

  * [1] M. J. Sobel, “The variance of discounted markov decision processes,” _Journal of Applied Probability_ , vol. 19, no. 4, p. 794–802, 1982.
  * [2] W. Dabney, M. Rowland, M. G. Bellemare, and R. Munos, “Distributional reinforcement learning with quantile regression,” 2017.
  * [3] V. Mnih, K. Kavukcuoglu, D. Silver, A. Graves, I. Antonoglou, D. Wierstra, and M. Riedmiller, “Playing atari with deep reinforcement learning,” _arXiv preprint arXiv:1312.5602_ , 2013.
  * [4] W. Dabney, G. Ostrovski, D. Silver, and R. Munos, “Implicit quantile networks for distributional reinforcement learning,” 2018.
  * [5] D. Yang, L. Zhao, Z. Lin, T. Qin, J. Bian, and T. Liu, “Fully parameterized quantile function for distributional reinforcement learning,” 2020.
  * [6] Y. Frégier and J.-B. Gouray, “Mind2mind: transfer learning for gans,” in _International Conference on Geometric Science of Information_. Springer, 2021, pp. 851–859.
  * [7] V. Kurin, S. Nowozin, K. Hofmann, L. Beyer, and B. Leibe, “The atari grand challenge dataset,” 2017.
  * [8] F. Zhuang, Z. Qi, K. Duan, D. Xi, Y. Zhu, H. Zhu, H. Xiong, and Q. He, “A comprehensive survey on transfer learning,” 2020.
  * [9] K. Simonyan and A. Zisserman, “Very deep convolutional networks for large-scale image recognition,” 2015.
  * [10] K. He, X. Zhang, S. Ren, and J. Sun, “Deep residual learning for image recognition,” 2015.
  * [11] D. J. Luckett, E. B. Laber, A. R. Kahkoska, D. M. Maahs, E. Mayer-Davis, and M. R. Kosorok, “Estimating dynamic treatment regimes in mobile health using v-learning,” _Journal of the American Statistical Association_ , 2019.
  * [12] C. Shi, S. Zhang, W. Lu, and R. Song, “Statistical inference of the value function for reinforcement learning in infinite horizon settings,” _arXiv preprint arXiv:2001.04515_ , 2020.
  * [13] Q. Liu, L. Li, Z. Tang, and D. Zhou, “Breaking the curse of horizon: Infinite-horizon off-policy estimation,” _Advances in Neural Information Processing Systems_ , vol. 31, 2018.
  * [14] M. R. Kosorok and E. B. Laber, “Precision medicine,” _Annual review of statistics and its application_ , vol. 6, pp. 263–286, 2019.
  * [15] I. J. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozair, A. Courville, and Y. Bengio, “Generative adversarial networks,” 2014.
  * [16] I. Gulrajani, F. Ahmed, M. Arjovsky, V. Dumoulin, and A. C. Courville, “Improved training of wasserstein gans,” _Advances in neural information processing systems_ , vol. 30, 2017.
  * [17] D. Freirich, T. Shimkin, R. Meir, and A. Tamar, “Distributional multivariate policy evaluation and exploration with the bellman gan,” in _International Conference on Machine Learning_. PMLR, 2019, pp. 1983–1992.
  * [18] C. Villani, _Optimal transport: old and new_. Springer, 2009, vol. 338.
  * [19] J. Von Neumann and O. Morgenstern, “Theory of games and economic behavior,” 2007\.
  * [20] R. E. Bellman, “Dynamic programming,” 1957.
  * [21] M. G. Azar, R. Munos, and B. Kappen, “On the sample complexity of reinforcement learning with a generative model,” _arXiv preprint arXiv:1206.6461_ , 2012.
  * [22] G. Barth-Maron, M. W. Hoffman, D. Budden, W. Dabney, D. Horgan, D. Tb, A. Muldal, N. Heess, and T. Lillicrap, “Distributed distributional deterministic policy gradients,” _arXiv preprint arXiv:1804.08617_ , 2018\.
  * [23] M. G. Bellemare, W. Dabney, and R. Munos, “A distributional perspective on reinforcement learning,” in _International Conference on Machine Learning_. PMLR, 2017, pp. 449–458.
  * [24] W. Dabney, M. Rowland, M. Bellemare, and R. Munos, “Distributional reinforcement learning with quantile regression,” in _Proceedings of the AAAI Conference on Artificial Intelligence_ , vol. 32, no. 1, 2018.
  * [25] M. Rowland, M. Bellemare, W. Dabney, R. Munos, and Y. W. Teh, “An analysis of categorical distributional reinforcement learning,” in _International Conference on Artificial Intelligence and Statistics_. PMLR, 2018, pp. 29–37.
  * [26] M. Mirza and S. Osindero, “Conditional generative adversarial nets,” _arXiv preprint arXiv:1411.1784_ , 2014.
  * [27] M. Arjovsky, S. Chintala, and L. Bottou, “Wasserstein generative adversarial networks,” in _International conference on machine learning_. PMLR, 2017, pp. 214–223.
  * [28] M. Haas and S. Richter, “Statistical analysis of wasserstein gans with applications to time series forecasting,” _arXiv preprint arXiv:2011.03074_ , 2020.
  * [29] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, Ł. Kaiser, and I. Polosukhin, “Attention is all you need,” _Advances in neural information processing systems_ , vol. 30, 2017.
  * [30] H. Wang, H. Luo, and Y. Wang, “Mbgdt:robust mini-batch gradient descent,” 2022\. [Online]. Available: https://arxiv.org/abs/2206.07139


# Detection of coupling in Duffing oscillator systems

Martin Brešar<EMAIL_ADDRESS>Department of Systems and Control, Jožef
Stefan Institute, Jamova cesta 39, SI-1000 Ljubljana, Slovenia  Jožef Stefan
International Postgraduate School, Jamova cesta 39, SI-1000 Ljubljana,
Slovenia  Pavle Boškoski<EMAIL_ADDRESS>Department of Systems and
Control, Jožef Stefan Institute, Jamova cesta 39, SI-1000 Ljubljana, Slovenia
Martin Horvat<EMAIL_ADDRESS>Faculty of mathematics and physics,
University of Ljubljana, Jadranska cesta 19, SI-1000 Ljubljana, Slovenia

###### Abstract

In complex dynamical systems, the detection of coupling and its direction from
observed time series is a challenging task. We study coupling in coupled
Duffing oscillator systems in regular and chaotic dynamical regimes. By
observing the conditional mutual information (CMI) based on the Shannon
entropy, we successfully infer the direction of coupling for different system
regimes. Moreover, we show that in the weak coupling limit the values of CMI
can be used to infer the coupling parameters by computing the derivative of
the conditional mutual information with respect to the coupling strength,
called the information susceptibility. The complete numerical implementation
is available at https://repo.ijs.si/mbresar/duffing-cmi.

††preprint: AIP/123-QED

> Complex systems encountered in nature can be divided into simpler subsystems
> that are coupled to each other, a prime example being a biological system.
> The problem of the detection of coupling from measured time series of
> subsystems has been of great interest in different branches of science. In
> this article, we study methods for the detection of coupling based on
> information transfer on a mathematical model. We show that information
> transfer estimation can be used to infer the exact value of the model
> coupling parameter.

##  I Introduction

In nature we encounter many complex systems that are divided into better known
subsystems which are coupled to each other. The subsystems can be
nonautonomous, i.e. time-dependent. In experimental practice, we typically
only have measurements of observables of the identified individual subsystems
in the form of signals. Studying these signals naturally raises a question of
what we can infer from them about the system structure, i.e., subsystem
properties and their coupling. In this paper, we focus predominantly on the
last property.

More precisely, we study the effects of time-independent coupling between two
nonautonomous systems through information transfer or flow between the
trajectories of the systems determined by the conditional mutual information
Thomas M. Cover (2006). As an example of such systems, we use the well known
Duffing oscillators Holmes (1979), i.e., periodically driven nonlinear
harmonic oscillators. Although information flow is usually observed in the
context of detecting the direction of coupling, the focus here is on the
absolute values of information flow in a mathematical model. The goal is to
solve the inverse problem of inferring the coupling parameter in different
coupled nonautonomous oscillator systems from the calculated information flow.

Estimating the conditional mutual information has been successfully applied in
numerous studies, including biological systems. For example, information flow
has been used in the analysis of EEG recordings of epileptic patients Paluš
_et al._ (2001), in rats under the influence of different anaesthetics Musizza
_et al._ (2007), and in the human cardiorespiratory system Paluš and Vejmelka
(2007). Studies have also been done on various theoretical systems, such as
simple phase oscillators Paluš and Stefanovska (2003), discrete Henon maps
Paluš  _et al._ (2001), Rossler and Lorenz systems Paluš and Vejmelka (2007),
and different chaotic, stochastic, and nonautonomous systems Clemson and
Stefanovska (2014). It has been shown that information flow can correctly
detect both the existence and the direction of coupling in two coupled, but
not yet synchronised systems. This means that the coupling must be
sufficiently small. While it has been shown that the information flow between
systems increases as the coupling increases, this relationship has not been
thoroughly investigated.

Another measure for quantifying information transfer between systems is the
transfer entropy proposed by Schreiber (2000). It has been shown that the
transfer entropy in a specific (although generally useful) setup is
mathematically equivalent to the conditional mutual information, see e.g.
(Paluš and Vejmelka, 2007, Appendix A). In a more general setup, the transfer
entropy can suffer from a lack of statistical significance of distributions
(due to high parameter dimensionality) used in calculations for short input
signals.

First, we study the numerical properties of the information flow calculation,
where the convergence was found to be faster than expected. Next, we consider
the behavior of the information flow when different pairs of Duffing systems
are coupled. Addressing the issue of the dependence of the information flow on
the coupling strength, we study the inverse problem, i.e., inferring the
coupling parameter from the information flow. The result is an estimate of the
information susceptibility as a derivative of the information flow. For
details, see Section V.3. It was found that the inverse problem is solvable in
the limit of strong coupling and for most cases of driven chaotic systems in
the limit of weak coupling. Finally, the range of the coupling parameter where
this is possible is given for individual pairs of systems.

##  II Inferring the direction of coupling

Consider two nonautonomous dynamical systems with trajectories denoted by
$\mathbf{x}(t)$ and $\mathbf{y}(t)$ and let the time evolution of such coupled
systems be given by equations

$\displaystyle\frac{d\mathbf{x}(t)}{dt}=\mathbf{f}_{x}(\mathbf{x},t)+\mathbf{g}_{x}(\mathbf{x},\mathbf{y},t)$
(1)
$\displaystyle\frac{d\mathbf{y}(t)}{dt}=\mathbf{f}_{y}(\mathbf{y},t)+\mathbf{g}_{y}(\mathbf{y},\mathbf{x},t),$
(2)

where $\mathbf{x}(t)$ and $\mathbf{y}(t)$ can be multidimensional and are
therefore written in bold letters. In the case where both $\mathbf{g}_{x}$ and
$\mathbf{g}_{y}$ are nonzero, the coupling is bidirectional. In our studies of
the inverse problem the coupling functions ${\bf g}_{x}$ and ${\bf g}_{y}$ are
time indepedent. If only one is nonzero, the coupling is unidirectional. In
the subsequent analysis, only unidirectional coupling is analysed. When
$\mathbf{g}_{y}=0$ and $\mathbf{g}_{x}\neq 0$, we call $\mathbf{x}(t)$ the
driven system and $\mathbf{y}(t)$ the driving system.

The trajectories of both subsystems $\mathbf{x}(t)$ and $\mathbf{y}(t)$ are
observed at equally spaced times $t_{i}=t_{0}+i\Delta t$ for $i\in\mathbb{Z}$.
Thus, we obtain time series of these trajetories
$(\mathbf{x}(t_{i}):i\in\mathbb{Z})$ and $(\mathbf{y}(t_{i}):i\in\mathbb{Z})$.
The methods presented in the following are used to determine the direction
and, in some cases, the magnitude of the coupling based on such measured time
series.

###  II.1 Conditional mutual information

In the context of information theory, the amount of information being
transferred from one system to another is quantified by the Shannon entropy
Thomas M. Cover (2006). The Shannon entropy of a random variable $X$ with a
probability distribution $p_{X}$ and its support $R_{X}$ (subset of the range
of $X$) is defined as

$H(X)=-\sum_{x\in R_{X}}p_{X}(x)\log p_{X}(x).$ (3)

The joint entropy of two random variables $X$ and $Y$ with a joint probability
distribution $p_{X,Y}$ is defined as

$H(X,Y)=-\sum_{x\in R_{X},y\in R_{Y}}p_{X,Y}(x,y)\log p_{X,Y}(x,y).$ (4)

Their mutual information $I(X,Y)$ , which is a measure of the amount of common
information contained in $X$ and $Y$, is defined as

$I(X,Y)=H(X)+H(Y)-H(X,Y)$ (5)

and is maximal if $X=Y$ and equal to $0$ if $X$ and $Y$ are independent. The
conditional entropy of $X$ and $Y$ is

$H(X|Y)=-\sum_{x\in R_{X},y\in R_{Y}}p_{X,Y}(x,y)\log p_{X|Y}(x|y),$ (6)

where $p_{X|Y}$ is the conditional probability distribution. Conditional
entropy measures the uncertainty in $X$, provided that everything about $Y$ is
known. Conditional entropy $H(X,Y|Z)$ is defined similarly, which leads to
conditional mutual information (CMI)

$I(X,Y|Z)=H(X|Z)+H(Y|Z)-H(X,Y|Z).$ (7)

It quantifies the average amount of common information contained in $X$ and
$Y$ given the value of $Z$.

While a measurement almost certainly contains uncorrelated noise that is
stochastic in nature, the underlying dynamics of the measured system can still
be deterministic.

Let us consider two possibly coupled systems with trajectories $\mathbf{x}(t)$
and $\mathbf{y}(t)$. Their corresponding time series are
$(\mathbf{x}(t_{i}):i\in\mathbb{Z})$ and $(\mathbf{y}(t_{i}):i\in\mathbb{Z})$.
Their values are distributed according to the probability distributions
$p_{\mathbf{X}}$ and $p_{\mathbf{Y}}$. We think of the values in the series
$\mathbf{x}(t_{i})$ and $\mathbf{y}(t_{i})$ as a realization of the random
variables denoted by $\mathbf{X}$ and $\mathbf{Y}$ at time $t_{i}$. By just
discussing their distributions, we lose the information about time. To restore
the notion of time, we need to introduce the time-lagged random variables
$\mathbf{X}_{\tau}$ and $\mathbf{Y}_{\tau}$ with realizations ${\bf
x}(t_{i+\tau})$ and ${\bf y}(t_{i+\tau})$ at times $t_{i}$ and discuss them in
parallel with their non time-lagged versions.

In order to infer a possible influence of the trajectory $\mathbf{x}(t)$ on
$\mathbf{y}(t)$ from their time series, we need to compute the mutual
information between the random variable of the first system $\mathbf{X}$ and
the time-lagged random variable of the second system $\mathbf{Y}_{\tau}$ given
by $I(\mathbf{X},\mathbf{Y}_{\tau})$. However, the mutual information captures
the connection between $\mathbf{X}$ and $\mathbf{Y}$ as well as between
$\mathbf{Y}$ and $\mathbf{Y}_{\tau}$ when the two systems are not independent.
Therefore, Paluš  _et al._ (2001) proposed to observe the conditional mutual
information $I(\mathbf{X},\mathbf{Y}_{\tau}|\mathbf{Y})$ to infer causal
relations between trajectories $\mathbf{x}(t)$ and $\mathbf{y}(t)$. This
quantity only captures the net information flow from $\mathbf{x}(t)$ to
$\mathbf{y}(t)$. A tweak to increase precision is to instead calculate

$I(\mathbf{X},\Delta_{\tau}\mathbf{Y}|\mathbf{Y})=H(\mathbf{X}|\mathbf{Y})+H(\Delta_{\tau}\mathbf{Y}|\mathbf{Y})-H(\mathbf{X},\Delta_{\tau}\mathbf{Y}|\mathbf{Y}),$
(8)

where $\Delta_{\tau}\mathbf{Y}=\mathbf{Y}_{\tau}-\mathbf{Y}$. This allows us
to look at small local changes in $\mathbf{y}(t)$ rather than the entire range
of data. Information flow from $\mathbf{y}(t)$ to $\mathbf{x}(t)$ is
calculated as $I(\mathbf{Y},\Delta_{\tau}\mathbf{X}|\mathbf{X})$ by swapping
$\mathbf{X}$ and $\mathbf{Y}$.

###  II.2 Information flow

Both $I(\mathbf{X},\mathbf{Y}_{\tau}|\mathbf{Y})$ and
$I(\mathbf{X},\Delta_{\tau}\mathbf{Y}|\mathbf{Y})$ are zero for $\tau=0$.
However, $\lim_{\tau\rightarrow
0^{+}}I(\mathbf{X},\mathbf{Y}_{\tau}|\mathbf{Y})=0$ and $\lim_{\tau\rightarrow
0^{+}}I(\mathbf{X},\Delta_{\tau}\mathbf{Y}|\mathbf{Y})\neq 0$. Since we are
interested in small time lags in Section IV and V, we use the second
definition.

To obtain a numerically more robust measure, we can average the CMIs
$I(\mathbf{X},\Delta_{\tau}\mathbf{Y}|\mathbf{Y})$ over some time lags $\tau$
in order to better approximate the information flow between two systems Paluš
_et al._ (2001).

$I(\mathbf{x}\rightarrow\mathbf{y})=\frac{1}{l}\sum_{i=1}^{l}I(\mathbf{X},\Delta_{i}\mathbf{Y}|\mathbf{Y})$
(9)

The appropriate values of $l$ are discussed in Section IV.1.

When calculating information flow from one subsystem to another, we use the
values of trajectories. The trajectories capture the full extend of
information about the dynamics of a particular system and can be
multidimensional. In realistic situations, however, the measurement of a
system or a subsystem is usually one-dimensional. A possible solution to
account for a multidimensional phase space is the use of time-delay embedding
vectors Rand (1981) or instantaneous phase representation in the case of
oscillatory systems Rosenblum, Pikovsky, and Kurths (1996); Paluš and
Stefanovska (2003).

##  III Numerical analysis

Information theory is thoroughly researched and well understood. However,
numerical computation of entropies is challenging due to limitations and
necessary approximations. It is therefore important to be aware of them and
study their effect on known systems in order to be able to analyze unknown
ones.

###  III.1 Numerical implementation

A convenient way to numerically implement the calculation of CMI is to use the
chain rule for entropies $H(X|Y)=H(X,Y)-H(Y)$ and rewrite (7) as

$I(\mathbf{X},\Delta_{\tau}\mathbf{Y}|\mathbf{Y})=H(\mathbf{X},\mathbf{Y})+H(\Delta_{\tau}\mathbf{Y},\mathbf{Y})-H(\mathbf{Y})-H(\mathbf{X},\Delta_{\tau}\mathbf{Y},\mathbf{Y})\>.$
(10)

The biggest issue in calculating CMI (10) from the involved entropies is the
estimation of the underlying probability distributions by the corresponding
histograms. The latter are computed from long enough series of trajectories:
having a time series of length $n$, we assign its values in each dimension
into $N$ bins, where $N\ll n$, and construct a histogram. Multidimensional
histograms are used to estimate multidimensional probability distributions.
Consider a system with two $k$-dimensional subsystems. When computing the
joint entropy in CMI (10), the histogram used in its calculation has $N^{3k}$
bins. Individual bins must contain a significant number of data points to
sufficiently approximate the distribution. Therefore, the total number of bins
cannot be too large.

The simplest way to construct a histogram is to use equidistant binning, which
means that the widths of all histogram bins are equal. We found this binning
procedure to be sufficient for our problem. Other binning procedures can also
be used, such as marginal equiquantization, where the marginal bins are chosen
so that they each contain the same number of data points Paluš and Vejmelka
(2007).

###  III.2 Duffing oscillator systems

The model we have chosen for numerical analysis is the Duffing oscillator
system, which is studied in detail in Kovacic and Brennan (2011). Its
properties make it very convenient, since it is a nonautonomous oscillator
that can be in a regular or in a chaotic dynamical regime. It is defined by

$\frac{d^{2}x}{dt^{2}}+\delta\frac{dx}{dt}+\alpha x+\beta
x^{3}=\gamma\cos(\omega t).$ (11)

It is a nonlinear periodically forced oscillator with damping. Its trajectory
is denoted by $(x(t),\frac{dx}{dt}(t))$ and the associated phase space is two-
dimensional. The choice of the model parameters
$(\delta,\alpha,\beta,\gamma,\omega)$ yields different systems which can be in
a regular or in a chaotic dynamical regime Ueda (1979). This allows us to test
the behavior of CMI on systems in different regimes.

For regular behavior, a high damping $\delta=0.2$ and a linear potential
$\beta=0$ were chosen. The orbits of the chosen regular systems are circles,
same as in normal driven harmonic oscillators, while chaotic systems exhibit
richer dynamics. For chaotic behavior, it is beneficial to have low damping.
In fact, increasing the damping parameter $\delta$ of a chaotic system makes
it regular. Four different chaotic systems denoted $C_{1}$ to $C_{4}$ and two
regular ones denoted $R_{1}$ and $R_{2}$ were chosen. The parameters and the
maximal Lyapunov exponents of the chosen systems are presented in Table 1. The
maximal Lyapunov exponents were calculated from the stroboscopic map at
$\cos{\omega t}=1$ using the Benettin algorithm Benettin  _et al._ (1980). The
negative maximal Lyapunov exponents of regular systems were obtained due to
orbits converging towards a limit cycle. Two examples of chaotic systems are
shown in Figure 1.

Table 1:  Parameters of the six Duffing systems under consideration and their maximal Lyapunov exponents $\lambda_{\rm max}$. The systems $R_{1},R_{2}$ are in a regular regime. The chaotic systems $C_{1},C_{2}$ have a high nonlinear term and $C_{3},C_{4}$ have a low nonlinear term. System | $\delta$ | $\alpha$ | $\beta$ | $\gamma$ | $\omega$ | $\lambda_{\rm max}$
---|---|---|---|---|---|---
$R_{1}$ | 0.2 | 1 | 0 | 1 | 0.51 | -1.23(2)
$R_{2}$ | 0.2 | 1 | 0 | 1 | 1.83 | -0.34(3)
$C_{1}$ | 0 | 1 | 400 | 1 | 1.1 | 0.6(2)
$C_{2}$ | 0 | 1 | 80 | 1 | 3.1 | 0.2(3)
$C_{3}$ | 0.3 | -1 | 1 | 0.5 | 1.2 | 0.5(7)
$C_{4}$ | 0.2 | -1 | 1 | 0.3 | 1 | 0.9(9)
Figure 1: Chaotic Duffing oscillator systems $C_{1}$ (left) and $C_{4}$
(right). The orange dots represent a stroboscopic map at $\cos(\omega t)=1$
obtained for long times $t_{max}=10^{5}$.

Systems $C_{1}$ and $C_{2}$ are similar to each other, and so are $C_{3}$ and
$C_{4}$, as can be seen from the table of parameters. Systems $C_{1}$ and
$C_{2}$ have a highly nonlinear potential (high $\beta$) and no damping
($\delta=0$), while $C_{3}$ and $C_{4}$ have a small $\beta$ and a high
$\delta$, making their dynamics less "wild" but still chaotic. In fact, the
system $C_{4}$ has a fractal attractor with dimension of around $1.4$
Tarnopolski (2013). The driving frequencies $\omega$ of the different systems
were chosen so that one is not a multiple of another.

###  III.3 Duffing system’s CMIs

Two Duffing systems can be coupled in different ways. We have chosen a
unidirectional linear coupling of the form $+\epsilon_{1}(x_{2}-x_{1})$ due to
its intuitive interpretation. Such a coupling forces the driven system towards
the driving system and can cause synchronization. If the coupling is very
strong, the systems become identical. Coupled Duffing systems are thus defined
as:

$\displaystyle\frac{d^{2}x_{1}}{dt^{2}}+\delta_{1}\frac{dx_{1}}{dt}+\alpha_{1}x_{1}+\beta_{1}x_{1}^{3}=\gamma_{1}\cos(\omega_{1}t)+\epsilon_{1}(x_{2}-x_{1})$
(12)
$\displaystyle\frac{d^{2}x_{2}}{dt^{2}}+\delta_{2}\frac{dx_{2}}{dt}+\alpha_{2}x_{2}+\beta_{2}x_{2}^{3}=\gamma_{2}\cos(\omega_{2}t).$
(13)

While the term driven can also be used to describe an individual Duffing
oscillator, i.e., it is driven by an external periodic force, we use the terms
driving and driven to distinguish between $x_{1}$ and $x_{2}$. We say that the
driven system $x_{1}$ depends on the driving system $x_{2}$.

We write a trajectory of the coupled Duffing system as
$(x_{1}(t),\frac{dx_{1}}{dt}(t),x_{2}(t),\frac{dx_{2}}{dt}(t))$ and
trajectories of its subsystems as ${\bf
x}_{i}(t)=(x_{i}(t),\frac{dx_{i}}{dt}(t))$ for $i=1,2$. The phase space of the
whole system is thus four-dimensional. In our calculations we ensured that our
system is on the attractor, so that the results do not depend on the initial
conditions. This was accomplished by taking the initial conditions
$(x_{1}(0),\frac{dx_{1}}{dt}(0),x_{2}(0),\frac{dx_{2}}{dt}(0))=(0,0,0,0)$ and
only considering the data from the trajectories from time $t=100$ onward. This
way, the energy of the system is gained entirely from the nonautonomous parts
of subsystems (12) and (13).

The trajectory of the whole system, or more precisely, the corresponding time
series, was obtained by first transforming the equations of motion (12) and
(13) to a system of four first order differential equations and then
integrating them using the Runge-Kutta 4 integrator with a fixed integration
step equal to $\Delta t=0.002$ and the time series length of $n=10^{7}$
throughout our study. The use of a constant integration step is convenient for
calculating probability distributions of points on the trajectory. In the
calculation of the histogram of points pertaining to the evaluation of CMI, we
use equidistant binning for each two-dimensional subsystem, as explained in
Section III.1. In this study we use $4$ bins for $x$ direction and $4$ for
$dx/dt$ direction unless otherwise specified. This effectively gives $16$ bins
per subsystem.

The convergence of CMI with the length of the trajectory $t$, more precisely
the number of steps in the sample $n(\approx t/\Delta t)$, is shown in Figure
2. In the direction without coupling (a), the CMI values decrease with the
length of the time series $n$ and stabilise at a nonzero value. We attribute
this to approximation errors of the discrete probability distributions. In the
direction of coupling (b), the values converge to a constant. The convergence
is slower in the presence of chaotic systems, especially for $C_{1}$ and
$C_{2}$. In a certain range of $\epsilon_{1}$ we see that the dependence is
linear in the log-log scale. Fitting a function of the form $f(t)=At^{-C}$ in
this range, we obtain the coefficient $C=1$ for convergence, even for chaotic
systems.

This is a bit surprising, since CMI is essentially a difference of entropies
determined by the empirical distributions which converge for periodic and
stochastic (chaotic) signals as $O(t^{-1})$ and $O(t^{-1/2})$, respectively.
The latter was checked in our study. The $O(t^{-1})$ convergence of CMI is
therefore faster than of its ingredients. This seems to be a consequence of
the fact that the stochastic corrections to entropies of order of magnitude
$O(t^{-1/2})$ are correlated and cancel each other out in the calculation of
CMI, and only the next possible order $O(1/t)$ remains. This is good news for
the usability of CMI, as even shorter samples can deliver relatively precise
results. Further discussion of the convergence is given in Appendix A.

Figure 2: The convergence of CMIs (8) for time lag $\tau=1$ between different
systems in the direction without coupling (a) and in the direction of coupling
(b). The coupling parameter is $\epsilon_{1}=0.5$, the integration step is
$\Delta t=0.002$ and $4$ bins per scalar variable were used in the histograms.

##  IV Analysis of CMI and information flow

Figure 3: The dependence of CMI in both directions on the time lag $\tau$.
Four different systems are examined from (a) to (d).

We have defined CMI between two subsystems and information flow as the average
of CMIs over time lags. In this section, we study the behavior of these
quantities for various combinations of coupled Duffing systems. Mainly, we
examine the dependence of CMI on the time lag and the dependence of
information flow on the coupling parameter. We denote the direction of
coupling by $(C)$ and the direction without coupling by $(N)$, e.g.,
information flow in the direction of coupling is denoted by $I(C)$. Here, $C$
stands for coupled and $N$ stands for not coupled.

###  IV.1 Dependence of CMI on the time lag

First, it is important to observe the dependence of CMIs (8) on the time lag
$\tau$ in order to determine the averaging in the definition of information
flow (9). It is shown for four examples of coupled Duffing systems in Figure
3.

Three main conclusions can be drawn from these dependences. First, we note
that the values in the direction of coupling CMI $(C)$ (upper row) are much
larger than CMI $(N)$ (lower row), as can be seen in all examples from (a) to
(d). Therefore, the CMI values can be used to determine the direction of
coupling. Second, we see that CMI $(C)$ increases with increased coupling
$\epsilon_{1}$. This means that our measure of coupling between systems
increases with coupling. Finally, CMI $(C)$ decreases with increased time lag.
It falls with the correlation time of the time-lagged system. On the contrary,
CMI $(N)$ increases with increased time lag. An exception is the system
$C_{2}\rightarrow C_{1}$, where CMI $(C)$ increases again to an even higher
value at $\tau\Delta t\approx 1$ and then proceeds to decrease.

This shows that the appropriate averaging in the definition (9) is over small
time lags $\tau$, where the values of CMI are large and information about the
driving system is still present in the driven one. The averaging decreases
errors due to dynamic fluctuations and possible noise. The simulation of the
system is noiseless and the errors of dynamical fluctuations are decreased
with the length of the time series. Therefore, taking only one value of CMI is
shown to be sufficient. Since CMI $(C)$ generally falls and CMI $(N)$
generally increases for small time lags, the values at the smallest possible
time lag $\tau=1$ are taken into account. Our measure of information flow in
either direction is therefore the very first point in Figure 3.

Since oscillator systems are similar to themselves when they are time-lagged,
we also detect oscillations at the frequency of the time-lagged system. They
are most expressed for the system $R_{1}\rightarrow C_{2}$. It is important to
mention that while CMI generally decreases with $\tau$ for small coupling, the
dependence can be different for bigger coupling.

###  IV.2 Dependence of information flow on the coupling parameter

Figure 4: The dependence of the information flow (9) on the coupling parameter
$\epsilon_{1}$. The orange lines represent the flow in the direction of
coupling $I(C)$ and the blue lines in the direction without coupling $I(N)$.
The dependences for weak coupling are shown for various systems in (a) to (g).
For the system $R_{1}\rightarrow C_{3}$, strong coupling is also shown in (h).
The calculations were done for $1000$ evenly spaced values of $\epsilon_{1}$
for each case.

We can now study the information flow between different systems to see how it
characterises them. The focus here is on observing the dependence of the
information flow on the coupling parameter. This dependence is shown in Figure
4 for various coupled Duffing systems.

The general behavior is similar for any two coupled systems. The information
flow $I(C)$ is very close to $0$ in the absence of coupling, since the systems
are independent. It then generally increases as the coupling parameter
increases. This can be seen in the examples from (a) to (f) in Figure 4.

For intermediate values, this dependence can have many local extremes until it
begins to decrease and reaches zero in the limit
$\epsilon_{1}\rightarrow\infty$, as indicated by Figure 4 (h). In this limit,
the systems are the same because the orbit of one system is pushed towards
that of the other, and there is no information transfer between identical
systems.

It can also be seen that $I(C)$ is much greater than $I(N)$. This means that
the direction of coupling can be easily inferred from the information flow.
There are some exceptions when it comes to bigger coupling, one of which can
be seen in Figure 4 (g). We have not found a general rule for which systems
exhibit this property. However, it seems to occur when similar pairs of
systems are coupled, for example for all four cases $C_{i}\rightarrow R_{j}$,
for $i,j\in\\{1,2\\}$.

In the limits of weak and strong coupling, the dependence can be monotonic,
which will be discussed in detail in the next section. A difference is
observed between chaotic systems with high and with low nonlinear terms
$\beta$. If the driven system is one with a high nonlinear term $C_{1}$ or
$C_{2}$, the dependence appears much more noisy, especially when two of these
systems are coupled. Monotonicity is harder to show or find in these cases,
which makes solving the inverse problem harder.

##  V The inverse problem

In this section, we are interested in whether the coupling parameter is
uniquely determined by the information flow $I(C)$ in some range of values,
i.e., whether this dependence is monotonic. If this is the case, the coupling
parameter can be inferred from the time series of two known systems by
calculating the information flow.

For most pairs of systems, there are two cases that exhibit monotonic
dependence of the information flow on the coupling parameter, as indicated in
Section IV.2. The first is the limit of strong coupling, i.e., for large
values of the coupling parameter. The second is the limit of weak coupling,
which is most often observed in nature. These limiting cases will be studied
in detail. For intermediate values, the inverse problem is generally not
solvable.

###  V.1 Limit of strong coupling

First, we consider the limiting case of large coupling parameters
$\epsilon_{1}$. As explained in Section IV.2, the information flow $I(C)$
decreases towards zero in this limit. The dependence for large $\epsilon_{1}$
is shown in Figure 5 in double logarithmic scale. Since the dependence in this
scale is linear for sufficiently large $\epsilon_{1}$, we know its form is

$f(\epsilon_{1})=A\epsilon_{1}^{-B}.$ (14)

Using least squares regression on points in the marked areas, the value $B=1$
was obtained for all cases. The obtained formulas hold for a wide range of
values of $\epsilon_{1}$. The value $B=1$ holds regardless of the number of
bins chosen, as does the range of $\epsilon_{1}$ with this dependence (checked
up to $24$ bins per variable). Thus, the inverse problem is solvable in the
limit of large coupling.

Figure 5: Information flow in the direction of coupling for large
$\epsilon_{1}$ for different coupled systems. Vertical lines mark the lower
boundaries of areas used for regression.

###  V.2 Limit of weak coupling

The case of weak coupling is much more interesting, since it is usually
encountered in physical systems. All numerical tests suggest that in most
cases the information flow grows monotonically with increasing $\epsilon_{1}$,
up to a value $\epsilon_{max}$ where the monotonicity ends. It may end due to
a distinct peak in the dependence, as in Figure 4 (f),(g), due to a discrete
jump as in Figure 4 (e), or due to the appearence of noise as in Figure 4 (c).
The analysis of different coupled Duffing systems gives different answers to
the solvability of the inverse problem.

In most cases where the monotonicity ends with a distinct peak in the
dependence on the coupling parameter, it is due to the synchronization of
subsystems. In a mathematical model, the appearance of synchronization can be
observed with conditional Lyapunov exponents (CLE), i.e., the Lyapunov
exponents of the driven subsystem Pecora and Carroll (2015). The dependence of
the information flow on $\epsilon_{1}$ and the dependence of CLE on
$\epsilon_{1}$ is shown in Figure 6 for one of these cases. Synchronization,
which is marked by the largest exponent becoming negative, occurs close to the
point where the dependence of the information flow suddenly changes, which is
similar to observations made in Paluš and Vejmelka (2007); Palus (2018). A
thorough survey on synchronization is given in Pikovsky, Rosenblum, and Kurths
(2001).

In cases of the discrete jumps in the information flow dependence, the
observation of phase state portraits and conditional Lyapunov exponents reveal
bifurcations of the driven subsystem. These cases are examined in detail in
Appendix B.

Figure 6: Information flow in both directions and the largest Lyapunov
exponent of the driven subsystem (conditional Lyapunov exponent) for system
$C_{3}\rightarrow C_{4}$. Gray lines represent the point where the exponent
reaches zero, i.e., the point of synchronization.

Let us now consider two coupled regular systems $R_{2}\rightarrow R_{1}$. The
dependence in Figure 4 (a) indicates a monotonic trend up to
$\epsilon_{max}\approx 0.7$. This value should not strongly depend on the
details of numerical calculation in order to determine the value of
$\epsilon_{max}$ as a property of a system.

By increasing the number of marginal bins $N$, the results approach the exact
analytic dependence for smooth variables Härdle  _et al._ (2004). However,
$\epsilon_{max}$ moves towards $0$ as the number of marginal bins of the
histograms increases, as can be seen in Figure 7. This is observed in all
cases when the driven system is regular. We have done the same test with two
coupled harmonic oscillators and obtained very similar results. Therefore, we
suspect that this is due to the singularity of the probability distribution
$p(x,dx/dt)$ of regular systems, since their orbits are circles in 2D space
and even a very weak coupling changes the probability distributions
significantly so that it cannot be explained by linear perturbation theory.

Figure 7: Information flow in the direction of coupling for unidirectionally
coupled systems $R_{2}\rightarrow R_{1}$, calculated with different numbers of
bins per variable $N$. The dependences are normalized to be easily comparable.

For driven chaotic systems, the value of $\epsilon_{max}$ did not change for
different numbers of histogram bins. The numerical analysis of
$R_{1}\rightarrow C_{2}$ is done in Figure 8. Increasing the number of
marginal bins $N$ while keeping the length of the time series $n$ constant
does not affect the general shape, as seen in Figure 8 (a). However, the
values increase at $\epsilon_{1}=0$ because the bins become poorly populated.
With increasing length $n$ at a high $N=16$, the value decreases again towards
$0$, as seen in Figure 8 (b). This shows that while a finer partition is
generally beneficial, taking too many bins will cause them to be poorly
populated and distorts the results.

Moreover, the dependence is found to converge towards a monotonic dependence
as $n$ is further increased. Due to numerical imprecisions, the dependence
will never be exactly monotonic when considering values of $\epsilon_{1}$ that
are very close to each other. Instead, we show monotonicity empirically by
doing ever more precise calculations. These can be done by increasing $n$, and
perhaps $N$. Decreasing the time step of the integration $\Delta t$ changes
the results only minimally.

Figure 8: Information flow in the direction of coupling for unidirectionally
coupled systems $R_{1}\rightarrow C_{2}$, calculated with data length
$n=10^{6}$ at different number of bins $N$ (a) and with $N=16$ at different
$n$ (b). In (a), the dependences are normalized to be easily comparable.

For most cases of driven chaotic systems, monotonicity has been shown in the
regime of weak coupling strength $\epsilon_{1}$. For the coupling of two
systems with high nonlinear term $C_{i}\rightarrow C_{j}$, for
$i,j\in\\{1,2\\}$, seen in Figure 4 (d), the dependence always contains a lot
of noise that does not decrease with longer time series. This is due to
bifurcations already at very small values of the coupling parameter.
Bifurcations of the driven system destroy the monotonicity, as explained in
the next section.

###  V.3 Information susceptibility

Table 2: Values of $\epsilon_{max}$ for various driven and driving systems. Cases where the monotonicity of the information flow in the low coupling regime could not be shown are marked with a question mark. Systems that exhibit non-continuous dependencies due to bifurcations are marked from B.1 to B.6 as a reference to the appendix where they are examined. |  | Driving system |
---|---|---|---
|  | $R_{1}$ | $R_{2}$ | $C_{1}$ | $C_{2}$ | $C_{3}$ | $C_{4}$ |
Driven system | $C_{1}$ | $\epsilon_{max}=1.(0)$ | $\epsilon_{max}=0.3(8)^{\ref{appendixR2vC1}}$ |  | $?$ | $\epsilon_{max}=1.(9)$ | $\epsilon_{max}=1.(8)$ |
$C_{2}$ | $\epsilon_{max}=3.(9)$ | $\epsilon_{max}=5.(4)$ | $?$ |  | $\epsilon_{max}=0.(9)$ | $\epsilon_{max}=0.(7)$ |
$C_{3}$ | $\epsilon_{max}=0.0(4)$ | $\epsilon_{max}=0.02(1)^{\ref{appendixR2vC3}}$ | $\epsilon_{max}=0.02(3)^{\ref{appendixC1vC3}}$ | $\epsilon_{max}=0.02(1)^{\ref{appendixC2vC3}}$ |  | $\epsilon_{max}=0.7(4)$ |
$C_{4}$ | $\epsilon_{max}=0.3(2)^{\ref{appendixR1vC4}}$ | $\epsilon_{max}=0.5(1)^{\ref{appendixR2vC4}}$ | $\epsilon_{max}=0.4(7)$ | $\epsilon_{max}=1.0(4)$ | $\epsilon_{max}=0.5(8)$ |  |

To determine the values of $\epsilon_{max}$, we introduce a smooth function
$I_{s}$ obtained by fitting splines to the dependence of the information flow
$I(C)$ on the coupling constant. This allows us to calculate the information
susceptibility, which we define as

$\chi_{I}(\epsilon_{1})=\frac{dI_{s}}{d\epsilon_{1}}(\epsilon_{1}).$ (15)

The smooth dependence allows a simple calculation of the derivative, which is
not possible for discrete points of the information flow. The information
susceptibility $\chi_{I}(\epsilon_{1})$ is positive for small $\epsilon_{1}$.
While it is positive, the dependence of CMI on $\epsilon_{1}$ is monotonic and
the inverse problem is solvable. The value of $\epsilon_{max}$ is then marked
by the change of sign of $\chi_{I}(\epsilon_{1})$, i.e., the first point at
which it is equal to $0$.

This was done for systems where a monotonic dependence on the coupling
parameter was found. The function UnivariateSpline from the package
scipy.interpolate in Python 3 was used to obtain $I_{s}$. Parameters chosen
were $k=3$ for cubic splines and a smoothing factor between $0.005\leq s\leq
0.1$.

For the cases where the driven system is either $C_{3}$ or $C_{4}$, this is
straightforward most of the times, as can be seen in Figure 9 (a). An odd
example $R_{2}\rightarrow C_{4}$ can also be seen. It contains a nonsmooth
dependence at $\epsilon_{1}=0.51$. Observing the phase space portraits of the
driven system has shown that this is due to bifurcations. For
$0.51<\epsilon_{1}<0.61$ and for $\epsilon_{1}>0.8$ the driven subsystem
$C_{4}$ is in a regular periodic regime. If the driven system’s dynamics
change significantly at some value of $\epsilon_{1}$ before a distinct peak in
the information flow is observed (in this case at $\epsilon_{1}=0.51$), this
value marks $\epsilon_{max}$. Such cases are shown in Appendix B. A peak in
the information flow can still be seen in such cases, but does not mark
$\epsilon_{max}$.

For the cases seen in Figure 9 (b) where either $C_{1}$ or $C_{2}$ is the
driven system, the exact value of $\epsilon_{max}$ could not be determined due
to a noisy dependence already seen in Figure 4 (c). When these two systems are
driven, their orbits can get stuck in a region of phase space for some time.
This can have a seemingly random effect on the information flow, which results
in the noisiness. The monotonicity could only be reliably shown up to a
certain point where this noise started. In these cases, this point was taken
as $\epsilon_{max}$.

Once a smooth function is obtained, the derivative in (15) is simple to
calculate. Some examples of the information susceptibility and the obtained
values of $\epsilon_{max}$ are shown in Figure 10. The values of $\chi_{I}$
are always smaller in the cases of driven $C_{1},C_{2}$ than for driven
$C_{3},C_{4}$, implying that these systems are less adaptable to the driving
system.

Figure 9: Information flow data and their smoothed fit as a function of the
coupling parameter $\epsilon_{1}$ for different pairs of systems, driven
system having a small (a) and a large (b) nonlinear term. Figure 10:
Information susceptibility (15) for different systems. Dashed vertical lines
mark the obtained values of $\epsilon_{max}$.

Finally, we show the results for the solvability of the inverse problem in the
limit of weak coupling in Table 2. Only driven chaotic systems are considered,
since the results for driven regular systems depend strongly on the partition
of the phase space. Some properties of the individual systems and the
information flow can be derived from the table.

Most notably, a difference between chaotic systems with a large and with a
small nonlinear term is seen. For driven chaotic systems with a large
nonlinear term, $\epsilon_{max}$ was mostly marked by the appearence of noise.
Nevertheless, the values of $\epsilon_{max}$ are generally larger than for
driven $C_{3}$, $C_{4}$. This means that $C_{1}$ and $C_{2}$ are harder to
adapt to the driving system, which is also implied by the smaller values of
information flow and information susceptibility. For driven $C_{3}$ or
$C_{4}$, the noisy dependence was never observed and the values converged
faster.

There are six cases in which $\epsilon_{max}$ is marked by bifurcations in the
driven subsystem. In all these cases except $R_{1}\rightarrow C_{4}$, there is
a discrete jump in the information flow at the point of system change.

In summary, we have shown that the coupling parameter of most systems is
uniquely determined by the information flow in the limit of weak coupling. In
measured time series, the models of the measured systems are usually unknown
and surrogate data are often needed to determine the direction of coupling. In
coupled Duffing systems, we can determine both the direction and the magnitude
of the coupling in the limit of weak coupling. This is possible because the
model of the systems is known, or rather known not to change. Only the
coupling parameter changes, while the other system parameters are kept
constant. If this was known to be the case in a realistic system, inferring
the direction as well as the magnitude of coupling might be feasible solely by
calculating the information flow.

##  VI Conclusion

One of the challenges in studying complex systems is to infer the coupling
between individual systems from their time series. We studied the behavior of
information-based coupling detection methods on coupled Duffing oscillator
systems. The information flow successfully determined the direction of
coupling between chaotic or regular systems. This was always possible for
small coupling, and in most cases also for large coupling.

To find the numerical properties, the convergence of CMI with the length of
the time series and the dependence on different partitions were observed. The
convergence was found to be faster than expected, and a finer partition to be
generally beneficial as long as the individual bins are well populated.

While the absolute value of CMI is considered of little value in experimental
data, we can analyze it thoroughly in mathematical models. Thus, we studied
the dependence of CMI on the coupling parameter. We evaluated the usefulness
of CMI in solving the inverse problem of inferring the coupling parameter. By
defining the information susceptibility, we found that in some cases the
coupling parameter is uniquely determined by the information flow. This is the
case in the limit of strong coupling, and at least in some cases in the limit
of weak coupling. Although this knowledge cannot be applied to arbitrary
experimental data, it might be possible if accurate models of the systems were
known.

The inverse problem proved to be complex because bifurcations can occur when
systems are coupled. The exact classification of systems for which it is
solvable remains to be seen.

###### Acknowledgements.

The authors MB and PB acknowledge the financial support from the Slovenian
Research Agency (research program No. P2–0001). MH acknowledges the financial
support from the Slovenian Research Agency (research program No. P1–0402 and
research projects No. J1–1698).

##  Appendix A Convergence of CMI

The expected convergence of histograms with time $p(t)=p_{0}+At^{-C}$ is
$C=1/2$ for chaotic and $C=1$ for regular systems. In Figure 2 it was found
that the convergence of the information flow is $C=1$ for all systems. To find
out why this is the case, we tested the convergence of different variables for
the system $C_{2}\rightarrow C_{1}$. For this test, we chose $\epsilon_{1}=0$,
since we expect CMI to converge close to $0$.

We looked at the convergence of the conditional mutual information (8) towards
$0$, the four entropies involved in the calculation of CMI (10) and the
probability distributions used to calculate these entropies. In Figure 11, the
probability distributions are denoted as $p_{1}$ to $p_{4}$ and the entropies
as $H_{1}$ to $H_{4}$, according to the order in (10). The entropy $H_{3}$ has
two arguments, the entropies $H_{1}$ and $H_{2}$ have four arguments, and the
entropy $H_{4}$ has six arguments. The same holds for the corresponding
distributions.

As expected, the empirical probability distributions $p_{i}$ converge with
$C=1/2$. This was verified by treating the distributions as vectors and
observing their convergence with increasing length of the series $n$ in 2-norm
$|\cdot|_{2}$. The converged values are denoted by $p_{\infty}$ and
$H_{\infty}$. The entropies $H_{i}$ also converge slower than with $C=1$.
However, taking the sum and difference of the entropies defining the
information flow CMI, the convergence is evidently faster than the
constituting terms with $C=1$.

Figure 11: The convergence of probabilities (a), entropies (b) and conditional
mutual information (time-lagging $C_{2}$ with $\tau=1$) (c) with the length of
the time series $n$ for the system $C_{2}\rightarrow C_{1}$ with
$\epsilon_{1}=0$. The black lines are given to indicate the type of
convergence.

##  Appendix B Systems with bifurcations

Some systems exhibit bifurcations, i.e., sudden qualitative changes in the
system behavior caused by a small change in the coupling parameter
$\epsilon_{1}$. In such cases, observing the phase space and conditional
Lyapunov exponents helps in understanding the dependence of the information
flow and the solvability of the inverse problem.

###  B.1 System $R_{2}\rightarrow C_{1}$

The information flow dependence, the CLE dependence, and the phase space
portraits for the system $R_{2}\rightarrow C_{1}$ are shown in Figure 12. The
inverse problem is solvable for $\epsilon_{1}\leq 0.38$, where the dynamics of
the system change. The change can also be seen from smaller values of the CLE.

Figure 12: The information flow and CLE dependences on the coupling parameter
and the corresponding phase space portraits of the driven subsystem for
$R_{2}\rightarrow C_{1}$. The orange dots represent a stroboscopic map at
$\cos(\omega_{1}t)=1$. Figure 13: The information flow and CLE dependences on
the coupling parameter and the corresponding phase space portraits of the
driven subsystem for $R_{2}\rightarrow C_{4}$. The orange dots represent a
stroboscopic map at $\cos(\omega_{1}t)=1$.

###  B.2 System $R_{2}\rightarrow C_{4}$

The information flow dependence, the CLE dependence, and the phase space
portraits for the system $R_{2}\rightarrow C_{4}$ are shown in Figure 13. In
the regions $0.51<\epsilon_{1}<0.61$ and $\epsilon_{1}<0.8$, the driven
subsystem exhibits regular dynamics and the information flow takes higher
values, while the CLE takes negative values. The inverse problem is solvable
until the first appearance of regular dynamics at $\epsilon_{max}=0.51$.

###  B.3 System $C_{2}\rightarrow C_{3}$

The information flow dependence, the CLE dependence, and the phase space
portraits for the system $C_{2}\rightarrow C_{3}$ are shown in Figure 14. The
system dynamics change several times, the first time at
$\epsilon_{max}=0.021$. This is seen in an intense behaviour of the
information flow and the CLE values.

Figure 14: The information flow and CLE dependences on the coupling parameter
and the corresponding phase space portraits of the driven subsystem for
$C_{2}\rightarrow C_{3}$. The orange dots represent a stroboscopic map at
$\cos(\omega_{1}t)=1$.

###  B.4 System $R_{1}\rightarrow C_{4}$

The information flow dependence, the CLE dependence, and the phase space
portraits for the system $R_{1}\rightarrow C_{4}$ are shown in Figure 15. The
inverse problem is solvable up to $\epsilon_{max}=0.32$, where the dynamics of
the driven subsystem become regular, which is reflected in the negative values
of the CLE. The driven subsystem becomes periodic with a period of exactly
$100$, i.e. $t_{period}=100\cdot 2\pi/\omega_{1}$.

Figure 15: The information flow and CLE dependences on the coupling parameter
and the corresponding phase space portraits of the driven subsystem for
$R_{1}\rightarrow C_{4}$. The orange dots represent a stroboscopic map at
$\cos(\omega_{1}t)=1$.

###  B.5 System $R_{2}\rightarrow C_{3}$

The information flow dependence, the CLE dependence, and the phase space
portraits for the system $R_{2}\rightarrow C_{3}$ are shown in Figure 16. The
dynamics of the driven subsystem become regular at $\epsilon_{max}=0.022$,
which is reflected in the negative values of the CLE.

Figure 16: The information flow and CLE dependences on the coupling parameter
and the corresponding phase space portraits of the driven subsystem for
$R_{2}\rightarrow C_{3}$. The orange dots represent a stroboscopic map at
$\cos(\omega_{1}t)=1$. Figure 17: The information flow and CLE dependences on
the coupling parameter and the corresponding phase space portraits of the
driven subsystem for $C_{1}\rightarrow C_{3}$. The orange dots represent a
stroboscopic map at $\cos(\omega_{1}t)=1$.

###  B.6 System $C_{1}\rightarrow C_{3}$

The information flow dependence, the CLE dependence, and the phase space
portraits for the system $C_{1}\rightarrow C_{3}$ are shown in Figure 17. The
dynamics of the driven subsystem become regular at $\epsilon_{max}=0.023$,
which is reflected in the negative values of the CLE.

## Data availability

The data that support the analysis of this article have been generated by the
authors and can be fully reproduced from the repository
https://repo.ijs.si/mbresar/duffing-cmi Python code for studying coupling
between Duffing systems with conditional mutual information .

## References

  * Thomas M. Cover (2006) J. A. T. Thomas M. Cover,  _Elements of Information Theory_ (Wiley John + Sons, 2006).
  * Holmes (1979) P. Holmes, “A nonlinear oscillator with a strange attractor,” Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences 292, 419–448 (1979).
  * Paluš  _et al._ (2001) M. Paluš, V. Komárek, Z. Hrnčíř, and K. Štěrbová, “Synchronization as adjustment of information rates: Detection from bivariate time series,” Physical Review E 63 (2001), 10.1103/physreve.63.046211.
  * Musizza  _et al._ (2007) B. Musizza, A. Stefanovska, P. V. E. McClintock, M. Paluš, J. Petrovčič, S. Ribarič, and F. F. Bajrović, “Interactions between cardiac, respiratory and EEG-$\delta$ oscillations in rats during anaesthesia,” The Journal of Physiology 580, 315–326 (2007).
  * Paluš and Vejmelka (2007) M. Paluš and M. Vejmelka, “Directionality of coupling from bivariate time series: How to avoid false causalities and missed connections,” Physical Review E 75 (2007), 10.1103/physreve.75.056211.
  * Paluš and Stefanovska (2003) M. Paluš and A. Stefanovska, “Direction of coupling from phases of interacting oscillators: An information-theoretic approach,” Physical Review E 67 (2003), 10.1103/physreve.67.055201.
  * Clemson and Stefanovska (2014) P. T. Clemson and A. Stefanovska, “Discerning non-autonomous dynamics,” Physics Reports 542, 297 – 368 (2014).
  * Schreiber (2000) T. Schreiber, “Measuring information transfer,” Phys. Rev. Lett. 85, 461–464 (2000).
  * Rand (1981) D. A. Rand,  _Dynamical systems and turbulence, Warwick 1980 : proceedings of a symposium held at the University of Warwick 1979-80_ (Springer-Verlag, Berlin New York, 1981).
  * Rosenblum, Pikovsky, and Kurths (1996) M. G. Rosenblum, A. S. Pikovsky, and J. Kurths, “Phase synchronization of chaotic oscillators,” Phys. Rev. Lett. 76, 1804–1807 (1996).
  * Kovacic and Brennan (2011) I. Kovacic and M. J. Brennan, eds.,  _The Duffing Equation_ (John Wiley & Sons, Ltd, 2011).
  * Ueda (1979) Ueda, “Randomly transitional phenomena in the system governed by duffing’s equation,” Journal of Statistical Physics (1979).
  * Benettin  _et al._ (1980) G. Benettin, L. Galgani, A. Giorgilli, and J.-M. Strelcyn, “Lyapunov characteristic exponents for smooth dynamical systems and for hamiltonian systems: A method for computing all of them. part 2: Numerical application,” Meccanica 15, 21–30 (1980).
  * Tarnopolski (2013) M. Tarnopolski, “On the fractal dimension of the duffing attractor,” Romanian Reports in Physics 66 (2013).
  * Pecora and Carroll (2015) L. M. Pecora and T. L. Carroll, “Synchronization of chaotic systems,” Chaos: An Interdisciplinary Journal of Nonlinear Science 25, 097611 (2015).
  * Palus (2018) M. Palus, “Linked by dynamics: Wavelet-based mutual information rate as a connectivity measure and scale-specific networks,” in  _Advances in Nonlinear Geosciences_, edited by A. A. Tsonis (Springer, Cham, 2018) pp. 427–463.
  * Pikovsky, Rosenblum, and Kurths (2001) A. Pikovsky, M. Rosenblum, and J. Kurths,  _Synchronization: A Universal Concept in Nonlinear Sciences_, Cambridge Nonlinear Science Series (Cambridge University Press, 2001).
  * Härdle  _et al._ (2004) W. Härdle, M. Müller, S. A. Sperlich, and A. Werwatz,  _Nonparametric and Semiparametric Models_ (Springer, 2004).
  * (19) Python code for studying coupling between Duffing systems with conditional mutual information, https://repo.ijs.si/mbresar/duffing-cmi.