# TREET: TRansfer Entropy Estimation via Transformer

Omer Luxembourg, , Dor Tsur, , Haim Permuter All authors are with the
Department of Electrical and Computer Engineering of Ben-Gurion University of
the Negev.

###### Abstract

Transfer entropy (TE) is a measurement in information theory that reveals the
directional flow of information between processes, providing valuable insights
for a wide range of real-world applications. This work proposes Transfer
Entropy Estimation via Transformers (TREET), a novel transformer-based
approach for estimating the TE for stationary processes. The proposed approach
employs Donsker-Vardhan (DV) representation to TE and leverages the attention
mechanism for the task of neural estimation. We propose a detailed theoretical
and empirical study of the TREET, comparing it to existing methods. To
increase its applicability, we design an estimated TE optimization scheme that
is motivated by the functional representation lemma. Afterwards, we take
advantage of the joint optimization scheme to optimize the capacity of
communication channels with memory, which is a canonical optimization problem
in information theory, and show the memory capabilities of our estimator.
Finally, we apply TREET to real-world feature analysis. Our work, applied with
state-of-the-art deep learning methods, opens a new door for communication
problems which are yet to be solved.

###### Index Terms:

Deep Learning, Information Theory, Transfer Entropy, Transformers, Neural
Estimation, Communication Channels

##  I Introduction

Transfer entropy (TE), introduced by Schreiber [1], stands as a pivotal
information-theoretic measurement that captures the coupling dynamics within
temporally evolving systems [2]. Derived from the principles of mutual
information (MI), TE is distinctive for its inherent asymmetry, strategically
employed in diverse applications for causal analysis [3]. TE serves as a
robust measure of the directed, asymmetric information flow between two
stochastic processes. Specifically, TE quantifies the reduction in uncertainty
about prospective values, incorporating the past values of one process to
predict the future values of another [4].

Transformers [5] are a neural network (NN) architecture that in recent years
became prevalent in many fields of research, including time series forecasting
[6, 7, 8]. With the attention mechanism, transformers are avoiding traditional
recurrent flow of states, such in recurrent NNs (RNNs), with a bounded length
of input sequence, yet manage to achieve great results [9].

The versatility of TE is evident in various domains. In neuroscience, it has
proven to be effective in deciphering functional connectivity among neurons
and between neurons to various physical tasks [10, 11, 12]. Moreover, analysis
between visual sensors and movement actuators in embodied cognitive systems
can be one via TE [13]. TE proves indispensable as it empowers social networks
to delve into the intricate dynamics of social influence and helps to
encompass the spread of misinformation and scrutinize statistical causality in
various text-related events, including the spread of misinformation [14]. [15]
utilizes a variant of TE in conjunction with an RNN to predict the direction
of the US stock market, incorporating TE as input feature. Lately [16]
suggested a greedy algorithm for feature selection while leveraging the
connection between each feature and the target with TE.

###  I-A Estimation and Optimization of Transfer Entropy

Estimating information theoretic measures such as TE is challenging, due to
lack of samples and unknown underlying distribution of the data, that may lead
to larger errors in estimation [4]. Among the estimation methods, Kernel
Density Estimation (KDE) [17], $k$-nearest neighbors (KNN) [18, 19], and the
Kraskov-Stögbauer-Grassberger (KSG) technique are noteworthy, despite their
struggle with the bias-variance trade-off. Notably, Granger causality [20]
also serves as a foundational method for TE estimation, particularly in
Gaussian joint processes where it is equivalent to TE, highlighting its
effectiveness in quantifying information flow between sequences [21, 22].

Recent advancements have embraced neural estimation methods for information
theory measures, leveraging the Donsker-Varadhan (DV) variational formula for
Kullback-Leibler (KL) divergence to accurately estimate metrics like MI,
directed information (DI) rate, and TE [23, 24, 25, 26]. These methods,
including MI neural estimation (MINE), DI neural estimation (DINE), and the
intrinsic TE neural estimator (ITENE) for MI, DI rate and variant of TE
respectively, utilize NNs to solve optimization problems, facilitating the
estimation process even in complex scenarios, such as estimating the channel
capacity of continuous channels in the presence of memory. However, challenges
arise with larger context windows, necessitating adjustments for accurate
estimation in broader historical contexts.

NNs have revolutionized fields such as computer vision and natural language
processing, with transformers now leading time series analysis, overtaking
RNNs due to their superior performance [27, 7, 6, 8]. Despite the widespread
use of MINE for information measure estimation, it subjects to certain
limitations that undergone comprehensive scrutiny [28, 29, 30]. To address
these issues, our novel approach, inspired by DINE, utilizes the capabilities
of transformers. This method specifically tailors TE estimation to the
sequential dynamics of data, offering a significant improvement over
conventional techniques that overlook temporal relationships.

###  I-B Contributions

In this work, we introduce TREET, a novel approach for estimating TE using
transformers. Our method is grounded in the DV representation, leveraging
attention-based NNs adapted to meet the structural constraints of DV
optimization. Theoretically, we establish the consistency of TREET and develop
an auxiliary neural distribution generator (NDG) module to facilitate TE
optimization, utilizing transformers.

Empirically, we showcase the versatility of TREET across various applications.
We present our findings on channel capacity estimation, illustrating the
efficacy of the TREET and NDG in a joint optimization and estimation process,
supported by theoretical validations. Additionally, we highlight the memory
capabilities of TREET through a long memory channel example. A key empirical
contribution is our feature analysis of the Apnea dataset, illustrating the
potential of TREET for comprehensive real-world data analysis, paving the way
for broader applications in future research.

###  I-C Organizations

The paper is structured as follows, Section II provides background on
information theory and NNs. Section III presents the TREET, it’s theoretical
guarantees and practical implementations, while the optimization of TREET is
described in Section IV. Experimental results for channel capacity estimation,
memory analysis and features analysis on physiological data are shown in
Section V. Section VI concludes the paper and discuss future research
potentials.

##  II Background

In this section, we elaborate on the preliminaries necessary to present our
method. We familiarize the reader with our notation and provide the formal
definition of TE, then relate it to DI. Subsequently, we introduce the concept
of transformer NN, define them, and discuss the theorem of universal
approximation. Lastly, we present the use of NN as estimators for information
theory measurements.

###  II-A Notation

Calligraphic letters, such as $\mathcal{X}$, denote subsets of the
$d$-dimensional Euclidean space, $\mathbb{R}^{d}$. Expectations are
represented by $\mathbb{E}$, with all random variables defined in the
probability space $(\Omega,\mathcal{F},\mathbb{P})$. The collection of Borel
probability measures on $\mathcal{X}$ is indicated as
$\mathcal{P}(\mathcal{X})$, and $\mathcal{P}_{\mathsf{ac}}(\mathcal{X})$
specifically refers to those measures that are absolutely continuous with
respect to Lebesgue measure, with their densities denoted by lowercase $p$.
Random variables and vectors are uppercase, e.g., $X$, and stochastic
processes are in blackboard bold, e.g., $\mathbb{X}:=(X_{t})_{t\in\mathbb{Z}}$
for discrete time $t$.

The sequence of $l$ samples from time $t$ in process $\mathbb{X}$ is
$X_{t}^{t+l}:=[X_{t},\ldots,X_{t+l}]^{\top}$, and for stationary processes,
$X^{l}:=[X_{1},\ldots,X_{l}]^{\top}$. For a measure
$\mathcal{Q}\in\mathcal{P}_{\mathsf{ac}}(\mathcal{X})$ with PDF $q$, cross
entropy between $P$ and $Q$ is
$\mathsf{h}_{\mathsf{CE}}(P,Q):=-\mathbb{E}_{P}[\log q]$. Differential entropy
of $X\sim P$ is $\mathsf{h}(X):=\mathsf{h}_{\mathsf{CE}}(P,P)$. If $Q\ll P$,
the KL divergence
$\mathsf{D}_{\mathsf{KL}}(P\|Q):=\mathbb{E}_{P}[\log\frac{\mathrm{d}P}{\mathrm{d}Q}]$.
MI between $(X,Y)\sim P_{XY}$ is
$\mathsf{I}(X;Y):=\mathsf{D}_{\mathsf{KL}}(P_{XY}\|P_{X}\otimes P_{Y})$.
Conditional KL divergence for $P_{Y|X},Q_{Y|X}$ given $X\sim P_{X}$ is
$\mathsf{D}_{\mathsf{KL}}(P_{Y|X}\|Q_{Y|X}|P_{X})$, and conditional MI for
$(X,Y,Z)\sim P_{XYZ}$ is
$\mathsf{I}(X;Y|Z):=\mathsf{D}_{\mathsf{KL}}(P_{XY|Z}\|P_{X|Z}\otimes
P_{Y|Z}|P_{Z})$.

###  II-B Transfer Entropy

TE quantifies causal influence from the past of one sequence on the present of
another, formally given by

######  Definition 1 (Transfer Entropy)

For jointly distibuted processes $\mathbb{X}$ and $\mathbb{Y}$ and
$k,l<\infty$, the TE, with parameters $(k,l)$, is given by

$\mathsf{TE}_{X\to
Y}(t;k,l):=\mathsf{I}\left(X^{t}_{t-l};Y_{t}|Y^{t-1}_{t-k}\right).$ (1)

When the considered processes are jointly stationary, we can omit the temporal
index $t$ in(1) and write TE as

$\mathsf{TE}_{X\to Y}(k,l):=\mathsf{I}\left(X^{l};Y_{l}|Y^{l-1}_{l-k}\right).$
(2)

Note that our definition of TE includes $X_{t}$, whereas other definitions [1,
31, 21, 26] define TE as $\mathsf{I}(X^{l-1};Y_{l}|Y^{l-1}_{l-k})$. The main
difference of this change is that our definition reduces to MI when the
processes are jointly independent and identically distributed (i.i.d.), i.e.,
$\mathsf{TE}_{X\to Y}(0,0)=I(X;Y)$. When $k=l$, we use the abbreviated
notation $\mathsf{TE}_{X\to Y}(l)$.

###  II-C Relation to Directed Information

DI [32, 33] is given by

$\displaystyle\mathsf{I}\left(X^{n}\rightarrow{}Y^{n}\right)$
$\displaystyle:=\sum_{i=1}^{n}\mathsf{I}\left(X^{i};Y_{i}|Y^{i-1}\right)$
$\displaystyle=\sum_{i=1}^{n}\mathsf{TE}_{X\to Y}(i).$ (3)

When $(X^{n},Y^{n})$ are $n$-fold projections of some stochastic processes
$\mathbb{X}$ and $\mathbb{Y}$ and from (2) we can observe that DI is the sum
of TE. We can define the DI rate, which is given by

$\displaystyle\mathsf{I}(\mathbb{X}\to\mathbb{Y})$
$\displaystyle:=\lim_{n\to\infty}\frac{1}{n}\mathsf{I}\left(X^{n}\rightarrow{}Y^{n}\right)$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{=}}\lim_{n\to\infty}\mathsf{I}(X^{n};Y_{n}|Y^{n-1}),$
(4)

and $(a)$ holds due to stationarity [25]. The DI rate can be therefore
observed as a limit case of TE, where the limit of TE exists since
$\mathsf{TE}_{X\to
Y}(l)=\mathsf{h}(Y_{l}|Y^{l-1})-\mathsf{h}(Y_{l}|X^{l},Y^{l-1})$ and each
conditional entropy is non-negative, decreasing for increasing $l$
(conditioning reduces entropy). Moreover, under appropriate Markov
assumptions,
$Y_{l}-Y^{l-1}-Y^{-1}_{-\infty},Y_{l}-(X^{l},Y^{l-1})-(X^{-1}_{-\infty},Y^{-1}_{-\infty})$,
we can form an equality between TE and DI (refer to Appendix VII-A for further
details).

###  II-D Attention Mechanism and Transformers

NNs serve as a universal function approximators [34], excelling in capturing
complex relationships within the data. NNs capabilities can be enhanced by the
attention111This paper considers the dot-product attention. mechanism [27],
which selects various combinations of inputs according to their significance
to the predictions. The attention can address time series datasets, while
weighting over each time input. Attention comprises queries, which represent
the current temporal focus, keys, which match against the queries to determine
relevance, and values, that contain the NN inputs to be weighted and act as a
memory function in time-series. Attention is formally given by

######  Definition 2 (Attention)

Let $Q=XW_{Q},K=XW_{K}$ and $V=XW_{V}$
($W_{Q},W_{K},W_{V}\in\mathbb{R}^{x\times d}$) be the queries, keys and
values, which are linear projections of the input
$X=(x_{1},x_{2},\ldots,x_{n})$, $\forall i:x_{i}\in\mathbb{R}^{x}$. Then the
attention is as follows,

$\mathsf{Attn}(X)=\mathsf{softmax}\left(QK^{\top}\right)V,$

where $\mathsf{softmax}(Z)_{j}={\exp{[Z_{j}]}}/{\sum_{m=1}^{n}\exp{[Z_{m}]}}$
where $Z\in\mathbb{R}^{n}$, is performed for each column of the dot-product
between queries and keys separately. The weights are calculated from the
softmax operation and the output is weighted from the values $V$.

Transformers utilize positional encoding (PE), which is a maps each input to
the attention according to its ordinal index, and is usually applied before
the attention for time series applications to deformation of the sequence
structure. In addition, a generalization of the transformer uses multi-head
attention,

$\mathsf{MultiHead}\text{-}\mathsf{Attn}(X)=[H_{1},\ldots,H_{h}]W_{0},$ (5)

where

$H_{i}=\mathsf{softmax}\left(Q^{(i)}{K^{(i)}}^{\top}\right)V^{(i)},\quad
i=1,\ldots,n,$ (6)

where $W_{0}\in\mathbb{R}^{dh\times d}$ is learnable parameter, and
$Q^{(i)},K^{(i)},V^{(i)}$ are the $i^{\text{th}}$ learnable projections of
queries keys and values [27]. A transformer block is a sequence-to-sequence
function, which maps input sequence of to an output sequence. Each transformer
block consists of attention layer and time-wise feed-forward layer.

Formally, a transformer is given as follows.

######  Definition 3 (Transformer function class)

Let $d_{i},d_{o},l,v\in\mathbb{N}$. The class of transformers with $v$
neurons, denoted
$\mathcal{G}_{\mathsf{tf}}^{(d_{i},d_{o},l,v)}:\mathbb{R}^{d_{i}\times
l}\to\mathbb{R}^{d_{o}\times l}$, is the set of discrete-time with the
following structure:

$\displaystyle X_{\mathsf{pe}}=W_{1}\cdot X+E,$ (7a)
$\displaystyle\mathsf{Attn}(X_{\mathsf{pe}})=X_{\mathsf{pe}}+\sum_{i=1}^{h}W_{O}^{i}W_{V}^{i}X_{\mathsf{pe}}$
$\displaystyle\qquad\qquad\cdot\mathsf{softmax}\left[(W_{K}^{i}X_{\mathsf{pe}})^{\top}(W_{Q}^{i}X_{\mathsf{pe}})\right],$
(7b) $\displaystyle
Y=\mathsf{FF}(X_{\mathsf{pe}})=\mathsf{Attn}(X_{\mathsf{pe}})+W_{3}\cdot\sigma_{\mathsf{R}}$
$\displaystyle\qquad\qquad\left(W_{2}\cdot\mathsf{Attn}(X_{\mathsf{pe}})+b_{1}\mathbf{1}_{l}^{\top}\right)+b_{2}\mathbf{1}_{l}^{\top},$
(7c)

where $X\in\mathbb{R}^{d_{i}\times l}$ is the input sequence of $l$ samples,
$Y\in\mathbb{R}^{d_{o}\times l}$ is the transformer output,
$W_{1}\in\mathbb{R}^{d_{e}\times d_{i}},W_{O}^{i}\in\mathbb{R}^{d_{e}\times
d_{m}},W_{Q}^{i},W_{K}^{i},W_{V}^{i}\in\mathbb{R}^{d_{m}\times
d_{e}},W_{2}\in\mathbb{R}^{d_{r}\times d_{e}},W_{3}\in\mathbb{R}^{d_{e}\times
d_{r}},b_{1}\in\mathbb{R}^{d_{r}},b_{2}\in\mathbb{R}^{d_{e}}$ are the weights
and biases of the network, $E\in\mathbb{R}^{d_{e}\times l}$ is the PE of the
input. The number of heads $h$ and the head size $d_{m}$ are the parameters of
the multi-head attention while $d_{r}$ is the hidden dimension of the feed-
forward (FF) layer. Assuming that the input is an additive product with it’s
positional encoding product, before the transformer blocks. The class of
transformers with dimensions $(d_{i},d_{o},l)$ is thus given by

$\mathcal{G}_{\mathsf{tf}}^{(d_{i},d_{o},l)}:=\bigcup_{v\in\mathbb{N}}\mathcal{G}_{\mathsf{tf}}^{(d_{i},d_{o},l,v)}.$
(8)

Transformers are a universal approximation class of sqeuence to sequence
mappings:

######  Theorem 1 (Universal approx. for transformers)

Let $\epsilon>0$, $l\in\mathbb{N}$. $\mathcal{U}\subset\mathbb{R}^{T\times
d_{i}},\mathcal{Z}\subset\mathbb{R}^{l\times d_{o}}$ be open sets, and
$f:\mathcal{U}\to\mathcal{Z}$ be a continuous vector-valued function. Then,
there exist $v\in\mathbb{N}$ and a $v$-neuron transformer
$g\in\mathcal{G}_{\mathsf{tf}}^{(d_{i},d_{o},l,v)}$ (as in Definition 3, such
that for any sequence of inputs $\\{u^{l}\\}\in\mathcal{U}$ and sequence of
outputs $\\{z^{l}\\}\in\mathcal{Z}$, we have

$\left\|f(u^{l})-g(u^{l})\right\|_{1}\leq\epsilon,$ (9)

This paper utilizes the class of causal transformers,
$\mathcal{G}_{\mathsf{ctf}}$, which is built upon $\mathcal{G}_{\mathsf{tf}}$.
To this end, we use the notion of causal functions

######  Definition 4 (Causal Function)

Let $\mathcal{F}:\mathbb{R}^{d_{u}\times L}\to\mathbb{R}^{d_{z}\times L}$ be a
function for $d_{u},d_{z},L\in\mathbb{N}$. For a series of inputs
$U=\\{u_{t}\\}_{t=1}^{L}$ and function outputs $Z=\\{z_{t}\\}_{t=1}^{L}$, the
function $\mathcal{F}$ is defined as a causal function iff, for any
$t_{0},t\in\\{1,\ldots,L\\},t\leq t_{0}$, the output $z_{t_{0}}$ depends only
on the inputs $u_{t}$.

In order to achieve a causal mapping function, we apply causal mask on the
attention scores, to result with a mapping that overlooks dependence on future
elements as in Definition 4 (note that only the attention perform time mixing,
thus leaving the rest of the transformers intact is valid). The causal mask,
denoted as $M\in\mathbb{R}^{l\times l}$, multiplied element-wise with the dot-
product of keys and queries, before the softmax operation, and is given by

$M_{[i,j]}=\begin{cases}1&\text{if}j\leq i\\\
-\infty&\text{otherwise}\end{cases}$ (10)

where $-\infty$ nullifies the corresponding entries after the softmax
operation. Hence, (7b) is written as,

$\displaystyle\mathsf{Attn}(X_{\mathsf{pe}})=X_{\mathsf{pe}}+\sum_{i=1}^{h}W_{O}^{i}W_{V}^{i}X_{\mathsf{pe}}$
$\displaystyle\qquad\cdot\mathsf{softmax}\left[(W_{K}^{i}X_{\mathsf{pe}})^{\top}(W_{Q}^{i}X_{\mathsf{pe}})\odot
M\right].$ (11)

Although changing the dot-product operation of the attention, the model
remains consistent with the universality framework, which is grounded in the
architecture’s capacity for pair-wise operations rather than being constrained
by the specifics of the attention mechanism [5]. Thus, zeroing out elements in
the input sequence aligns with the transformers function class.

###  II-E Neural Estimation

Neural estimation is a methodology that utilizes NN optimization for the
optimization of various functionals. First introduced by the MINE, neural
estimators often utilize the DV representation [23, Theorem 3.2].

######  Theorem 2 (DV representation)

For any, $P,Q\in\mathcal{P}(\mathcal{X})$, we have

$\mathsf{D}_{\mathsf{KL}}(P\|Q)=\sup_{f:\mathcal{X}\rightarrow{}\mathbb{R}}\mathbb{E}_{P}[f]-\log\left(\mathbb{E}_{Q}[e^{f}]\right),$
(12)

where the supremum is taken over all measurable functions $f$ with finite
expectations.

The MINE is then obtained by approximating the DV optimization function class
with the class of NNs, and estimating expectations with samples means. The
MINE was generalized to additional information theoretic quantities, such as
conditional MI [35], total correlation [36], and TE with its variations [26].
Furthermore, a generalization of the MINE to DI was proposed in [25]. Drawing
inspiration from the MI and DI neural estimators, we develop a provably
consistent estimator of TE that utilizes the power of the attention mechanism.

##  III Estimation of Transfer Entropy

In this paper we harness to computational power of transformers architecture
and modify the attention mechanism to result with TREET, a new estimator of TE
for high dimensional continuous data. We begin by deriving the estimator, then
account for its theoretical guarantees. Finally, we provide implementation
details, outlining the modifications of attention to neural estimation.

###  III-A Estimator Derivation

The TREET provides an estimate of the TE (1) from a set of samples
$D_{n,l}=(X^{n+l},Y^{n+l})\sim P_{X^{n+l},Y^{n+l}}$. We assume that $l\ll n$
and this omit the dependence on $l$ in the dataset notation. The TREET
decomposes TE into two KL terms, each of which taken between the the data
distribution and a corresponding reference distribution $\widetilde{P}_{Y}$,
such that $\widetilde{P}_{Y}$ can be sampled freely. To derive TREET, first we
represent TE as subtraction of KL divergences w.r.t. and absolutely continuous
$\widetilde{P}_{Y}$. We propose the following.

######  Lemma 1 (TE as KL Divergences)

TE decomposes as

$\mathsf{TE}_{X\to
Y}(l)=\mathsf{D}_{Y_{l}|Y^{l-1}X^{l}\|\widetilde{Y}_{l}}-\mathsf{D}_{Y_{l}|Y^{l-1}\|\widetilde{Y}_{l}},$
(13)

where

$\displaystyle\mathsf{D}_{Y_{l}|Y^{l-1}\|\widetilde{Y}_{l}}:=\mathsf{D}_{\mathsf{KL}}\left(P_{Y_{l}|Y^{l-1}}\|\widetilde{P}_{Y_{l}}\Big{|}P_{Y^{l-1}}\right),$
(14a) $\displaystyle\mathsf{D}_{Y_{l}|Y^{l-1}X^{l}\|\widetilde{Y}_{l}}:=$
$\displaystyle\qquad\qquad\mathsf{D}_{\mathsf{KL}}\left(P_{Y_{l}|Y^{l-1}X^{l}}\|\widetilde{P}_{Y_{l}}\Big{|}P_{Y^{l-1}X^{l}}\right),$
(14b)

and the conditional KL divergence is
$\mathsf{D}_{\mathsf{KL}}(P_{X|Z}\|P_{Y|Z}|P_{Z}):=\mathbb{E}_{Z}[\mathsf{D}_{\mathsf{KL}}(P_{X|Z}\|P_{Y|Z})]$.

Lemma 1 is proved in Section VII-C, and follows basic information theoretic
properties of TE and KL divergences. Utilizing the DV variational
representation (12) and performing the optimization over a set of causal
transformer architectures
$\mathcal{G}_{\mathsf{ctf}}^{Y}:=\mathcal{G}_{\mathsf{ctf}}^{(d_{y},1,l,v_{y})},\mathcal{G}_{\mathsf{ctf}}^{XY}:=\mathcal{G}_{\mathsf{ctf}}^{(d_{x}+d_{y},1,l,v_{xy})}$
for given $l,d_{y},d_{x},v_{y},v_{xy}\in\mathbb{N}$, and let
$g_{y}\in\mathcal{G}_{\mathsf{ctf}}^{Y},g_{xy}\in\mathcal{G}_{\mathsf{ctf}}^{XY}$.
Each KL divergence can be approximated with transformers by

$\displaystyle\mathsf{D}_{Y_{l}|Y^{l-1}\|\widetilde{Y}_{l}}=$
$\displaystyle\sup_{g_{y}\in\mathcal{G}_{\mathsf{ctf}}^{Y}}\mathbb{E}\left[g_{y}\left(Y^{l}\right)\right]$
$\displaystyle-\log\left(\mathbb{E}\left[e^{g_{y}\left(\widetilde{Y}_{l},Y^{l-1}\right)}\right]\right),$
(15a) $\displaystyle\mathsf{D}_{Y_{l}|Y^{l-1}X^{l}\|\widetilde{Y}_{l}}=$
$\displaystyle\sup_{g_{xy}\in\mathcal{G}_{\mathsf{ctf}}^{XY}}\mathbb{E}\left[g_{xy}\left(Y^{l},X^{l}\right)\right]$
$\displaystyle-\log\left(\mathbb{E}\left[e^{g_{xy}\left(\widetilde{Y}_{l},Y^{l-1},X^{l}\right)}\right]\right).$
(15b)

Finally, replacing expectations with sample means in (15), we result with the
TREET, given by

$\displaystyle\widehat{\mathsf{TE}}_{X\to Y}(D_{n};l)$
$\displaystyle:=\sup_{g_{xy}\in\mathcal{G}_{\mathsf{ctf}}^{XY}}\sup_{g_{y}\in\mathcal{G}_{\mathsf{ctf}}^{Y}}\widehat{\mathsf{TE}}_{X\to
Y}(D_{n},g_{y},g_{xy};l)$ (16a)
$\displaystyle=\sup_{g_{xy}\in\mathcal{G}_{\mathsf{ctf}}^{XY}}\widehat{\mathsf{D}}_{Y_{l}|Y^{l-1}X^{l}\|\widetilde{Y}_{l}}(D_{n},g_{xy})$
$\displaystyle\quad-\sup_{g_{y}\in\mathcal{G}_{\mathsf{ctf}}^{Y}}\widehat{\mathsf{D}}_{Y_{l}|Y^{l-1}\|\widetilde{Y}_{l}}(D_{n},g_{y}),$
(16b)

where

$\displaystyle\widehat{\mathsf{D}}_{Y_{l}|Y^{l-1}\|\widetilde{Y}_{l}}(D_{n},g_{y}):=\frac{1}{n}\sum_{i=1}^{n}g_{y}\left(Y_{i}^{i+l}\right)$
$\displaystyle\quad\qquad\qquad\qquad-\log\left(\frac{1}{n}\sum_{i=1}^{n}e^{g_{y}\left(\widetilde{Y}_{i+l},Y_{i}^{i+l-1}\right)}\right),$
(17a)
$\displaystyle\widehat{\mathsf{D}}_{Y_{l}|Y^{l-1}X^{l}\|\widetilde{Y}_{l}}(D_{n},g_{xy}):=\frac{1}{n}\sum_{i=1}^{n}g_{xy}\left(Y_{i}^{i+l},X_{i}^{i+l}\right)$
$\displaystyle\quad\qquad\qquad-\log\left(\frac{1}{n}\sum_{i=1}^{n}e^{g_{xy}\left(\widetilde{Y}_{i+l},Y_{i}^{i+l-1},X_{i}^{i+l}\right)}\right).$
(17b)

where $\widetilde{Y}$ is i.i.d. under absolutely continuous reference
measurement under the alphabet $\mathcal{Y}$, and
$\sup_{g_{xy}\in\mathcal{G}_{\mathsf{ctf}}^{XY}}\widehat{\mathsf{D}}_{Y_{l}|Y^{l-1}X^{l}\|\widetilde{Y}_{l}}(D_{n},g_{xy})=\mathsf{D}_{Y_{l}|Y^{l-1}X^{l}\|\widetilde{Y}_{l}}$,
$\sup_{g_{y}\in\mathcal{G}_{\mathsf{ctf}}^{Y}}\widehat{\mathsf{D}}_{Y_{l}|Y^{l-1}\|\widetilde{Y}_{l}}(D_{n},g_{y})=\mathsf{D}_{Y_{l}|Y^{l-1}\|\widetilde{Y}_{l}}$.
The TREET (16b) is consequently given as the subtraction of the solutions of
two optimization problems (17), each optimizing its own NN. The optimization
is then performed via gradient ascent with the corresponding model
$g_{y},g_{xy}$ and the dataset $D_{n}$. Note that Definition 15 states that
each KL divergence we estimate with DV representation, achieved by function of
$(Y^{l},\widetilde{Y}_{l})$ and $(Y^{l},X^{l},\widetilde{Y}_{l})$, thus we use
a causal transformer to opt $g_{y},g_{xy}$ in order to ignore relative future
values of the process. In the following section, we prove that the TREET
implemented with causal transformers is a consistent estimator of TE.

###  III-B Theoretical Guarantees

To further understand the capabilities of TREET, we explore its theoretical
guarantees. Assuming joint stationarity of $(\mathbb{X},\mathbb{Y})$ and using
the causal transformers function class, we introduce the consistency of the
proposed TREET.

######  Theorem 3 (TREET consistency)

Let $\mathbb{X}$ and $\mathbb{Y}$ be jointly stationary, ergodic stochastic
processes. TREET is strongly consistent estimator of $\mathsf{TE}_{X\to Y}(l)$
for $l\in\mathbb{N}$, i.e. $\mathbb{P}-a.s.$ for every $\epsilon>0$ there
exists and $N\in\mathbb{N}$ such that for every $n>N$ we have

$\left|\widehat{\mathsf{TE}}_{X\to Y}(D_{n};l)-\mathsf{TE}_{X\to
Y}(l)\right|\leq\epsilon$ (18)

where $l$ is the memory parameter of the TE.

The proof following the steps of representation step - represents TE as a
subtraction of two DV potentials, estimation step - proves that the DV
potentials is achievable by empirical mean of a given set of samples, and
approximation step - shows that the estimator built upon causal transformers
converges to TE with the corresponding memory parameter. The proof is given in
the Appendix VII-B.

ReferenceSamplerTREET$\theta_{y}$DV
Loss$\widehat{\mathsf{D}}_{Y_{l}|Y^{l-1}\|\widetilde{Y}_{l}}(D_{n},g_{\theta_{y}})$Eq.
(17a)$\widetilde{Y}_{t}$$Y^{t}_{t-l}$$g_{\theta_{y}}(Y_{t},Y_{t-l}^{t-1})$$g_{\theta_{y}}(\widetilde{Y}_{t},Y_{t-l}^{t-1})$$\nabla_{\theta_{y}}\widehat{\mathsf{D}}_{Y}(D_{n},g_{\theta_{y}})$

Figure 1: The estimator architecture for the calculation of
$\widehat{\mathsf{D}}_{Y_{l}|Y^{l-1}\|\widetilde{Y}_{l}}(D_{n},g_{\theta_{y}})$.

###  III-C Algorithm and Implementation

This section describes the TREET implementation. We present an overview scheme
for estimation of TE and describe the algorithm. Afterwards we show the
architecture of the TREET and going deeper to its implementation. Formally,
the transformer function class is parametric models whose finitely many
parameters are a subset of a parameter space
$\Theta\subset\mathbb{R}^{d},d\in\mathbb{N}$. For a fixed $v<\infty$ neurons,
denote the functions from the above classes as $g\in\mathcal{G}_{v}$, or their
corresponding parameterized form $g_{\theta},\theta\in\Theta$. We therefore
denote the corresponding networks with $g_{\theta_{y}},g_{\theta_{xy}}$, where
$\theta_{y}$ and $\theta_{xy}$ are the network parameters, respectively. We
now describe the algorithm design and transformer modifications.

####  III-C1 Overview and Algorithm

The TREET algorithm follows an iterative joint optimization of
$\widehat{\mathsf{D}}_{Y_{l}|Y^{l-1}X^{l}\|\widetilde{Y}_{l}}$ and
$\widehat{\mathsf{D}}_{Y_{l}|Y^{l-1}X^{l}\|\widetilde{Y}_{l}}$ through
iterative mini-batch gradient optimization. The algorithm inputs are
$l,D_{n}$, which are the TE parameter and the dataset, respectively. Every
iteration begins with feeding mini-batch sized $m<n$ with sequences length $l$
in each model, followed by the calculation of both DV potentials (17), that
construct $\mathsf{TE}_{X\to Y}(l)$ (16a). The calculated objective is then
used for gradient-based optimization of the NN parameters. The iterative
process continues until a stopping criteria is met, typically defined as the
convergence of $\mathsf{TE}_{X\to Y}(l)$ within a specified tolerance
parameter $\epsilon>0$. For evaluation, we produce the estimated TE again for
many samples as possible, and taking the averaged results. The full pipeline
for estimating
$\widehat{\mathsf{D}}_{Y_{l}|Y^{l-1}X^{l}\|\widetilde{Y}_{l}}(D_{n},g_{\theta_{y}})$
is presented in Figure 1, and the complete list of steps is given in Algorithm
1.

$g_{y}({Y_{l},Y^{l-1}})$ Dense Add & Norm Feed Forward Concat & Norm Fixed-
Past- Causal-Attention Embeddings+ Position Value
$\begin{bmatrix}Y_{1},\cdots,Y_{l}\end{bmatrix}$
KVQ$g_{y}({\widetilde{Y}_{l},Y^{l-1}})$ Dense Add & Norm Feed Forward Concat &
Norm Modified Fixed-Past- Causal-Attention Embeddings+ Position Value
$\begin{bmatrix}\smash{\widetilde{Y}_{1},\cdots,\widetilde{Y}_{l}}\end{bmatrix}$
$\widetilde{\textbf{K}}$$\widetilde{\textbf{V}}$$\smash{\widetilde{\textbf{Q}}}$K,V

Figure 2: The TREET architecture for $g_{y}$, with memory parameter $l$. It
illustrated using a single sequence as an example. However, it is capable of
parallel processing for sequences lengthed $L>l$, and in such cases, the
number of outputs for the function will be $L-l+1$. Both transformer share the
same weights. The FPCA and the modified FPCA are as elaborated in Section
III-C2. Algorithm 1 TREET

Input: Joint process samples $D_{n}$; Observation length $l\in\mathbb{N}$.
Output: $\widehat{\mathsf{TE}}_{X\to Y}(D_{n};l)$ \- TE estimation.

1:NNs initialization $g_{\theta_{y}}$, $g_{\theta_{xy}}$ with corresponding
parameters $\theta_{y},\theta_{xy}$.

2:Step 1 – Optimization:

3:repeat

4: Draw a batch $B_{m}$: $m<n$ sub-sequences, length $L>l$ from $D_{n}$, with
reference samples $P_{\widetilde{Y}}$ for each.

5: Compute both potentials
$\widehat{\mathsf{D}}_{Y_{l}|Y^{l-1}X^{l}\|\widetilde{Y}_{l}}(B_{m},g_{\theta_{xy}})$,
$\widehat{\mathsf{D}}_{Y_{l}|Y^{l-1}\|\widetilde{Y}_{l}}(B_{m},g_{\theta_{y}})$
via (17).

6: Update parameters:

7:
$\theta_{xy}\leftarrow\theta_{xy}+\nabla_{\theta_{xy}}\widehat{\mathsf{D}}_{Y_{l}|Y^{l-1}X^{l}\|\widetilde{Y}_{l}}(B_{m},g_{\theta_{xy}})$

8:
$\theta_{y}\leftarrow\theta_{y}+\nabla_{\theta_{y}}\widehat{\mathsf{D}}_{Y_{l}|Y^{l-1}\|\widetilde{Y}_{l}}(B_{m},g_{\theta_{y}})$

9:until convergence criteria.

10:Step 2 – Evaluation: Evaluate for a sub-sequence (17) and (16a) to obtain
$\widehat{\mathsf{TE}}_{X\to Y}(D_{n};l)$.

$\displaystyle QK^{\top}=\begin{bmatrix}{q}_{(t+l)}\cdot
k_{(t)}&\ldots&{q}_{(t+l)}\cdot{k}_{(t+l)}&-\infty&\ldots&-\infty\\\
-\infty&{q}_{(t+l+1)}\cdot
k_{(t+1)}&\ldots&{q}_{(t+l+1)}\cdot{k}_{(t+l+1)}&-\infty&\vdots\\\
\vdots&-\infty&\ddots&&\ddots&-\infty\\\
-\infty&\ldots&-\infty&{q}_{(t+L)}\cdot
k_{(t+L-l)}&\ldots&{q}_{(t+L)}\cdot{k}_{(t+L)}\\\ \end{bmatrix}$

(19)
$\scalebox{0.9}{\mbox{$\displaystyle\widehat{QK^{\top}}=\begin{bmatrix}\widetilde{q}_{(t+l)}\cdot
k_{(t)}&\ldots&\widetilde{q}_{(t+l)}\cdot\widetilde{k}_{(t+l)}&-\infty&\ldots&-\infty\\\
-\infty&\widetilde{q}_{(t+l+1)}\cdot
k_{(t+1)}&\ldots&\widetilde{q}_{(t+l+1)}\cdot\widetilde{k}_{(t+l+1)}&-\infty&\vdots\\\
\vdots&-\infty&\ddots&&\ddots&-\infty\\\
-\infty&\ldots&-\infty&\widetilde{q}_{(t+L)}\cdot
k_{(t+L-l)}&\ldots&\widetilde{q}_{(t+L)}\cdot\widetilde{k}_{(t+L)}\\\
\end{bmatrix}$}}.$ (20)

The DV representation, (12), suggest that the function $f$ is the same
function for both terms that construct it, i.e. the weights are shared for
both network propagation. Hence, our transformer in TREET, which constructed
by (17), is the same for both terms, for both DV representation. Exemplifying
with $\widehat{\mathsf{D}}_{Y_{l}|Y^{l-1}\|\widetilde{Y}_{l}}$, both
$g_{y}(Y^{l})$ and $g_{y}(\widetilde{Y}^{l},Y^{l-1})$ use the same learning
parameters, from positional encoding layers and attention to FF layers. The
only difference we have between the two terms, is in how we operate the
attention mechanism, which essentially re-use keys and values generated for
the first term, to generate the later one. This model for $g_{y}$ is
visualized in Figure 2. The following sections will elaborate about the
proposed fixed past causal attention (FPCA) that constructs the TREET and its
variation, modified fixed past causal attention which is required for the
reference measurement sampling.

####  III-C2 Fixed Past Causal Attention

In this section we describe the calculation of $g_{y}(Y^{l})$ (17a). For
$L>l$, the causal attention architecture is constructed by queries
$Q\in\mathbb{R}^{L\times d_{o}}$, keys $K\in\mathbb{R}^{L\times d_{o}}$ and
values $V\in\mathbb{R}^{L\times d_{o}}$, and given by

$\displaystyle Q=\begin{bmatrix}q_{1},\ldots,q_{L}\end{bmatrix}$
${}^{\top},\quad
K=\begin{bmatrix}k_{1},\ldots,k_{L}\end{bmatrix}^{\top},\quad$ $\displaystyle
V=$ $\displaystyle\begin{bmatrix}v_{1},\ldots,v_{L}\end{bmatrix}^{\top},\quad$
(21)

where $L\in\mathbb{N}$ is the temporal length of the model’s input and
$l\in\mathbb{N}$ is the memory parameter which determines $\mathsf{TE}_{X\to
Y}(l)$. The $i^{\text{th}}$ step queries, keys and values in the are given by
$q_{i}=W_{Q}x_{i},k_{i}=W_{K}x_{i},v_{i}=W_{V}x_{i}$, and
$M\in\mathbb{R}^{L\times L}$, thus we denote

$\mathsf{Causal}\text{-}\mathsf{Attention}:=\mathsf{softmax}(QK^{\top}\odot
M)V.$ (22)

To result in a valid causal attention mechanism and fixed length input size
$l$, we propose the FPCA, where each query will be multiplied only by its
relative time key and past $l-1$ keys. Thus the dot-product operation is
between $[t+l,t+L]$ queries time steps and $[t,t+L]$ keys time steps. In
addition, we change the causal mask (i.e. for the above case
$M\in\mathbb{R}^{L-l\times L}$) to a Toeplitz-like band matrix mask
$M^{\prime}\in\mathbb{R}^{L-l\times L}$,

$M^{\prime}_{[i,j]}=\begin{cases}1&\text{if}j-i<l\text{and}j\geq i\\\
-\infty&\text{otherwise}\end{cases}$ (23)

The resulting queries and keys dot-product with the given mask yields (19),
and the FPCA is given by

$\mathsf{FPCA}:=\mathsf{softmax}(QK^{\top}\odot M^{\prime})V$ (24)

Note that we use only the last $L-l+1$ results from the transformer outputs
sequence, in order to keep the past information fix to length $l$. With the
given FPCA, the required DV output $g_{y}(Y^{l})$ can be calculated.

####  III-C3 Reference Sampling with FPCA

This section describes the calculation of $g_{y}(\widetilde{Y}_{l},Y^{l-1})$
(17a). As mentioned in Section III-A, the reference distribution,
$\widetilde{Y}$, is independently drawn from some absolute continuous PDF over
the bounds of $\mathcal{Y}$222If $\mathcal{Y}$ is not bounded, we set manually
the maximal bounds regarding to the dataset properties.. Denote the scoring
value after the softmax operation of causal attention for query $q_{i}$ with
key $k_{j}$ as

$\mathcal{R}_{(i,j)}:=\frac{e^{q_{i}\cdot k_{j}}}{\sum_{m=1}^{i}e^{q_{i}\cdot
k_{m}}},\quad j\leq i.$ (25)

Since we use FPCA with memory parameter $l$, the output at time $t$ can be
written as
$\mathsf{FPCA}_{t}=\mathcal{R}_{(t,t-l)}v_{t-l}+\mathcal{R}_{(t,t-l+1)}v_{t-l+1}+\ldots+\mathcal{R}_{(t,t)}v_{t}$.
In order to calculate $g_{y}(\widetilde{Y}_{t},Y^{t-1}_{t-l})$, the output of
the time mixing must contain the original distribution from past information
up to time step $t-1$. This information is stored in the keys and values of
the previous resulted $\mathsf{FPCA}_{t}$, for $g_{y}(Y^{t}_{t-l})$. Thus, the
modified FPCA should be constructed by

$\widetilde{\mathcal{R}}_{(i,j)}:=\begin{cases}\frac{e^{\widetilde{q}_{i}\cdot
k_{j}}}{e^{\widetilde{q}_{i}\cdot\widetilde{k}_{i}}+\sum_{m=1}^{i-1}e^{\widetilde{q}_{i}\cdot
k_{m}}}&,i\neq j\\\\[5.0pt]
\frac{e^{\widetilde{q}_{i}\cdot\widetilde{k}_{i}}}{e^{\widetilde{q}_{i}\cdot\widetilde{k}_{i}}+\sum_{m=1}^{i-1}e^{\widetilde{q}_{i}\cdot
k_{m}}}&,i=j.\end{cases}$ (26)

In this case, the modified FPCA can be written as
$\mathsf{Modified}\text{-}\mathsf{FPCA}_{t}=\widetilde{\mathcal{R}}_{(t,t-l)}v_{t-l}+\widetilde{\mathcal{R}}_{(t,t-l+1)}v_{t-l+1}+\ldots+\widetilde{\mathcal{R}}_{(t,t)}\widetilde{v}_{t}$,
for time $t$. Summarizing, the second term (17a) generated by the modified
FPCA, contains all keys and values of the relative past and query, key and
value for the current present, which is a function of the reference
distribution. The dot-product matrix between queries and keys can be written
as (20), for $l\leq L$, and modified FPCA is written as
$\mathsf{Modified}\text{-}\mathsf{FPCA}=\mathsf{softmax}(\widehat{QK^{\top}}\odot
M^{\prime})V$.

The extension of FPCA and modified FPCA to multi-head version is immediate. In
our implementation, the reference inputs, for the second term in (17a) are
drawn from the uniform measure on the bounding box of the current batch of $Y$
samples, while theoretically, our method allows to draw from any positive
continuous distribution measure; for further details check Appendix VII-B. The
implementation of
$\widehat{\mathsf{D}}_{Y_{l}|Y^{l-1}X^{l}\|\widetilde{Y}_{l}}(D_{n},g_{xy})$
(17b) is obtained by concatenating the $X^{l}$ values with the corresponding
$Y^{l}$ or $\widetilde{Y}^{l}$ values for both the FPCA and the modified FPCA,
respectively.

##  IV Optimization of Estimated Transfer Entropy

Many applications can leverage the optimization of TE, in a given data-driven
setting. The optimization of which can be done via controlling the
distribution of the process $\mathbb{X}$ and $\mathbb{Y}$, independently or
jointly. An important example for such a setting would be communication
channels, whose capacity can be characterized with a certain TE term.

$\begin{bmatrix}X_{t-L}\\\
0\end{bmatrix}$$\begin{matrix}\cdots\end{matrix}$$\cdots$$\begin{bmatrix}X_{t-1}\\\
0\end{bmatrix}$$\begin{bmatrix}0\\\ U_{t}\end{bmatrix}$$X_{t}$NDGTransformer

Figure 3: The recursive process for the NDG with transformers. If feedback
presents, the input includes the past channel output, concatenated with the
corresponding value $X_{i}$ at time $i$.

######  Remark 1 (Channel capacity)

As presented in [25], consider channels with and without feedback links from
the channel output back to the encoder. The feedforward capacity of a channel
sequence $\left\\{P_{Y^{n}\|X^{n}}\right\\}$, for $n\in\mathbb{N}$ is

$C_{\mathsf{FF}}=\lim_{n\to\infty}\sup_{P_{X^{n}}}\frac{1}{n}\mathsf{I}(X^{n};Y^{n}),$
(27)

while the feedback capacity is

$C_{\mathsf{FB}}=\lim_{n\to\infty}\sup_{P_{X^{n}\|Y^{n-1}}}\frac{1}{n}\mathsf{I}(X^{n}\to
Y^{n}).$ (28)

The achievability of the capacities is further discussed in [37], [38]. [32]
showed that for non-feedback scenario, the optimization problem over
$P_{X^{n}\|Y^{n}}$ can be translated to $P_{X^{n}}$, which support the use of
DI rate for both optimization problems. Under some conditions TE can estimate
the capacity as well, as presented in Section II-C.

Focusing on channel capacity, we assume that we can control the input
distribution sampling mechanism and propose an algorithm for the optimization
of estimated TE with respect to the input generator. We refer to this model as
the Neural Distribution Generator (NDG). The estimated TE optimization
methodology is inspired by the proposed methods from [25]. However, the
adaptation to transformer architectures considers a different implementation
of the proposed scheme. As the TREET estimates TE from samples, the NDG is
defined as a generative model of the input distribution samples, and is
optimized with the goal of maximizing the downstream estimated TE.

######  Lemma 2 (Optimal TE)

Let $(\mathbb{X},\mathbb{Y})$ be jointly stationary processes, and the TE with
memory parameter $l\in\mathbb{N}$. Then, the maximal TE, $\mathsf{TE}_{X\to
Y}^{\star}(l)$, by $\mathbb{X}$ is

$\mathsf{TE}_{X\to
Y}^{\star}(l):=\sup_{P_{X^{l}}}\mathsf{I}(X^{l};Y_{l}|Y^{l-1}).$ (29)

The proof of the lemma is given in Appendix VII-D. The lemma suggests that NDG
with input sequence length $l$ is enough to achieve maximum TE with memory
parameter $l$, for independently controlling the distribution of $\mathbb{X}$.

The NDG calculates a sequence of channel input $X^{l}$ through the mapping

$h_{\phi}:(U_{i},Z^{i-1}_{i-l})\rightarrow X_{i}^{\phi},\qquad i=1,\ldots,n,$
(30)

where $U_{i}$ is the random noise drawn from
$P_{U}\in\mathcal{P}_{\mathsf{ac}}(\mathcal{U}),\mathcal{U}\subset\mathbb{R}^{d_{x}}$,
which cause the stochasticity, $Z_{i-l}^{i-1}$ are the past observation of the
generated process created by the model, and the channel corresponding outputs
if feedback exists, and $h_{\phi}$ is the parametric NDG mapping with
parameters $\phi\in\Phi$. By the functional representation lemma [39] and the
restated lemma in [25], we can achieve the distribution of $X$ from an NN
function. After $l$ iterations, with $Z_{i-l}^{i-1}$ storing the relevant
information of previous iterations, the NDG generates the whole sequence.
Transformers need access to the whole sequence at once, in contrast to RNNs
where a single state can theoretically represent the past sequence.

Thus, the input sequence to the NDG with transformer is created via past
outputs of the transformer itself (and corresponding channel outputs if
feedback exists) as depicted in Figure 3, and the past observations are taken
without gradients to prevent backpropagation through iterations. In our
implementation, the input convention for each time step is a concatenated
vector of $[X_{i},U_{i}]^{\top}$ to maintain the input structure of samples
and random noise for every projection in the network, while in relative
history time steps, the noise is replaced with zero, and for the present time,
the input $X_{i}$ is replaced with a zero vector of the same dimension.

NDG Transformer $\phi$Channel$P_{Y_{i}|X^{i}_{i-l},Y^{i-1}_{i-l}}$TREET
Transformer
$g_{\theta_{y}},g_{\theta_{xy}}$$X_{i-l}^{i-1}$$X_{i-l}^{i-1},Y_{i-l}^{i-1}$DV
Loss$U_{i}$$X_{i}$$Y_{i},X_{i}$$\nabla g_{\theta_{y}},\nabla
g_{\theta_{xy}}$$\nabla\phi$$\star$

Figure 4: Complete system for estimating and optimizing TREET with NDG while
Altering between the models to train on. ($\star$) If feedback presents, past
channel output realizations included. Algorithm 2 Continuous TREET
optimization

Input: Continuous sequence-to-sequence system $\mathcal{S}$; Observation
length $l\in\mathbb{N}$.
Output: $\widehat{\mathsf{TE}}_{X\to Y}^{\star}(U^{n};l)$, optimized NDG.

1:NNs initialization $g_{\theta_{y}}$, $g_{\theta_{xy}}$ and $h_{\phi}$ with
corresponding parameters $\theta_{y},\theta_{xy},\phi$.

2:repeat

3: Draw noise $U^{m}$, $m<n$.

4: Compute batch $B_{m}^{\phi}$ sized $m$ using NDG, $\mathcal{S}$

5: if training TREET then

6: Perform TREET optimization - Step 1 in Algorithm 1.

7: else Train NDG

8: Compute $\widehat{\mathsf{TE}}_{X\to
Y}(B_{m}^{\phi},g_{\theta_{y}},g_{\theta_{xy}},h_{\phi};l)$ using (16a).

9: Update NDG parameters:

10: $\quad\phi\leftarrow\phi+\nabla_{\phi}\widehat{\mathsf{TE}}_{X\to
Y}(B_{m}^{\phi},g_{\theta_{y}},g_{\theta_{xy}},h_{\phi};l)$

11:until convergence criteria.

12:Draw $U^{m}$ to produce $l$ length sequence and evaluate
$\widehat{\mathsf{TE}}_{X\to Y}(D_{n}^{\phi};l)$.

13:return $\widehat{\mathsf{TE}}_{X\to Y}^{\star}(U^{n};l)$, optimized NDG.

To estimate and optimize the TE at once, we jointly training the TREET and the
NDG. As described in Algorithm 2 and Figure 4, in each iteration we only
update one of those models by maximizing each DV potentials,
$\widehat{\mathsf{D}}_{Y_{l}|Y^{l-1}\|\widetilde{Y}_{l}},\widehat{\mathsf{D}}_{Y_{l}|Y^{l-1}X^{l}\|\widetilde{Y}_{l}}$,
for TREET model and by maximizing $\widehat{\mathsf{TE}}_{X\to Y}$ for NDG
model. The entire pipeline of one iteration will create a sequence length $l$
from the NDG, by iterating it $l$ times with some initiated zero values -
which creates the dataset, $D_{n}^{\phi}=(X^{\phi,n},Y^{\phi,n})$. Afterwards,
feeding each sample sequence from the dataset to TREET for achieving the
corresponding networks’ outputs as mentioned back in Algorithm 1, which
constructs the loss that will generate gradients backward according to each
model. Next, we demonstrate the power of the proposed method for estimating
the channel capacity.

##  V Experimental Results

In this section we demonstrate the utility and performance of the proposed
algorithms to several tasks, encompasing both estimation and optimization of
TE. Additionally, we perform a comparison with the RNN-based DI rate
estimation scheme from [25] and discusses the results. While the experiments
are about continuous spaces of distributions, discrete spaces can easily be
applied to the algorithm of TREET, but not for the optimization part with NDG
as presented. The TREET consists of transformer with a single fixed-past-
causal-attention layer follows by a single FF layer, and the NDG also consists
of transformer with a single causal-attention layer and a single FF layer.
Further details on the implementation parameters are provided in Appendix
VII-E. All simulations are implemented in PyTorch and the code implementation
can be found at this repository: https://github.com/omerlux/TREET.

$-1$$-0.5$$0$$0.5$$1$$0$$0.2$$0.4$EmpiricalProbabilityInput
Noise$-2$$0$$2$$0$$0.1$$0.2$$0.3$EmpiricalProbabilityOutput $X$

Figure 5: NDG input noise and output $X$, $0$ SNR case.

$-20$$-10$$0$$10$$0$$0.5$$1$$1.5$$\displaystyle
P/\sigma^{2}\mathrm{[dB]}$$\displaystyle\mathrm{Capacity}$TheoreticalDINE
est.TREET est.

(a) Capacity as function of SNR.

$1$$2$$3$$4$$5$$6$$0.5$$1$$1.5$$2$Vector Dimension
$\displaystyle\mathrm{(0dB)}$CapacityTheoreticalDINE est.TREET est.

(b) Capacity as function of vector dimension.

Figure 6: Channel capacity estimation of AWGN channel. 6(a) presents the
estimation of TREET and DINE in comparison to theoretical capacities for
various SNR rates. $P$ represents the power constraint and the noise parameter
is $\sigma^{2}$. 6(b) shows the capacity estimation with variate vector
process dimension for the $0$ dB case. For the case of power constraint $P=1$
theoretically the input $P_{X}$ distributes $\mathcal{N}(0,1)$.

$-20$$-10$$0$$10$$0$$0.5$$1$$1.5$$\displaystyle
P/\sigma^{2}\mathrm{[dB]}$Feedback CapacityTheoreticalDINE est.TREET est.

(a) Feedback Capacity.

$-20$$-10$$0$$10$$0$$0.5$$1$$1.5$$\displaystyle
P/\sigma^{2}\mathrm{[dB]}$Feedforward CapacityTheoreticalDINE est.TREET est.

(b) Feedforward Capacity.

Figure 7: Channel capacity estimation of Gaussian MA(1) channel with variate
SNR.

$-20$$-10$$0$$10$$0$$0.5$$1$$1.5$$\displaystyle
P/\sigma^{2}\mathrm{[dB]}$Feedback CapacityTheoreticalDINE est.TREET est.

(a) Feedback Capacity.

$-20$$-10$$0$$10$$0$$0.5$$1$$1.5$$\displaystyle
P/\sigma^{2}\mathrm{[dB]}$Feedforward CapacityTheoreticalDINE est.TREET est.

(b) Feedforward Capacity.

Figure 8: Channel capacity estimation of Gaussian AR(1) channel with variate
SNR.

###  V-A Empirical Capacity Estimation for Finite Memory Processes

As shown in Section II-C, under some conditions TE is equal to the DI rate and
converges to it. DINE [25] proved that it can estimate a capacity of
stationary channels by optimizing the DI rate estimator and the input
distribution of the channel $P_{X}$, as Remark 1 mentions about channel
capacities. Our experiments has shown that for memory-less channels and for
memory channels with and without feedback, TREET can approximate the capacity
with the joint optimization procedure of the input distribution (NDG) and TE
estimation. Important to note that channel input constraints are essential for
ensuring that the transmitted signals are well-suited to the channel’s
characteristics and limitations. In our experiment we applied the power
constraint on the input signal of the channel, which is implemented via
normalizing a batch of samples to a certain statistics according to the power
constraint.

####  V-A1 AWGN Channel

Consider additive white Gaussian noise (AWGN) channel with i.i.d. noise,

$\displaystyle Z_{i}$ $\displaystyle\sim\mathcal{N}(0,\sigma^{2}),$
$\displaystyle Y_{i}$ $\displaystyle=X_{i}+Z_{i},\qquad i\in\mathbb{Z},$

$X_{i}$ is the channel’s input sequence, coupled with the average power
constraint $\mathbb{E}[X_{i}^{2}]\leq P$. The capacity of this channel is
simple for analytical calculation, and is given by the following formula
$\text{C}=0.5\log\left(1+{P}/{\sigma^{2}}\right)$. Since the process is
memoryless, maximized both DI rate and TE (for any $l\geq 0$) coincide with
the channel capacity. We estimated and optimized the TREET according to
Algorithm 2 and compared the model performance with the DI rate estimation and
optimization scheme from [25]. Results are presented in Figure 6. It can be
seen that both are estimating the right capacities which are the MI, although
their access to multiple past observations, and in addition, the dimension of
the input vector and output vector of a channel changes the capacities, as
expected. Note that larger dimensions cause error of in estimation and still
is an open academic research [28]. To further analyse the learned
distribution, we visualize the optimized NDG mapping in Figure 5. It can be
seen that the optimized NDG maps the uniform ($U\sim[-1,1]$) inputs into
Gassian samples. This observation meets our expectations, as the capacity
achieving AWGN distribution is Gaussian [40].

####  V-A2 Gaussian MA(1) Channel

Given a Moving Average (MA) Gaussian noise channel with order 1

$\displaystyle Z_{i}$ $\displaystyle=N_{i}+\alpha N_{i-1},$ $\displaystyle
Y_{i}$ $\displaystyle=X_{i}+Z_{i},\qquad i\in\mathbb{Z},$

where $N_{i}\sim\mathcal{N}(0,\sigma^{2})$ are i.i.d. and $X_{i}$ is the input
to the channel with power constraint $\mathbb{E}[X_{i}^{2}]\leq P$. We apply
Algorithm 2 to both feedforward and feedback settings, comparing with ground
truth solutions and the DI-based scheme. The feedforward capacity can be
calculated with the water-filling algorithm [34], while the feedback capacity
can be calculated by the root of forth order polynomial equation [41]. As seen
in Figure 7 our method successfully estimated the capacity for a wide range of
SNR values.

$0$$20$$40$$60$$80$$100$$120$$5$$10$$15$$20$$25$Attention Relative Time Step
($t-i$)Series Index$0.02$$0.07$$0.12$Attention weight

(a) 130 input length.

$0$$15$$30$$45$$60$$75$$90$$5$$10$$15$$20$$25$Attention Relative Time Step
($t-i$)Series Index$0.010$$0.015$$0.020$Attention weight

(b) 90 input length.

Figure 9: Attention weights at training convergence of TREET optimized by NDG.
Each row in the matrices represent a different input sequence and the columns
are the weights of past values $i$ from current prediction $t$ (i.e. $i=0$
represent the present prediction $t$). For the GMA(100) process, it can be
observed that giving enough time-steps, TREET easily observes at the related
time $i=100$ where the information needed to be gathered from, while shorter
length lead to instability of training.

####  V-A3 Gaussian AR(1) Channel

The case of autoregressive (AR) Gaussian noise channel of order 1 is similar,

$\displaystyle Z_{i}$ $\displaystyle=N_{i}+\alpha Z_{i-1},$ $\displaystyle
Y_{i}$ $\displaystyle=X_{i}+Z_{i},\qquad i\in\mathbb{Z},$

where $N_{i}\sim\mathcal{N}(0,\sigma^{2})$ are i.i.d. and $X_{i}$ is the input
to the channel with power constraint $\mathbb{E}[X_{i}^{2}]\leq P$. The
capacity is also affected by the existence of feedback. Feedforward capacity
can be solved with the water-filling algorithm [34] and [42] prove how to
achieve the feedback capacity analytically. Figure 8 compares the results of
DINE and TREET estimators.

###  V-B Memory Capabilities Analysis

Previous experiments showed the capability of correctly estimating TE. Delving
deeper into the memory effectiveness of TREET, we tested Gaussian MA channel,
as presented before, but with time delay of 100 steps.

$\displaystyle Z_{i}$ $\displaystyle=N_{i}+\alpha N_{i-100},$ $\displaystyle
Y_{i}$ $\displaystyle=X_{i}+Z_{i},\qquad i\in\mathbb{Z}.$

It is notable that to achieve the correct capacity, the information from
$t-100$ steps must be an input to the TREET. We tested the estimation
optimization algorithm with shorter input length and longer input length, than
the demanded one. Figure 9(a) shows the analysis of the attention weights
related to the $\widehat{\mathsf{D}}_{Y_{l}|Y^{l-1}X^{l}\|\widetilde{Y}_{l}}$
model with $l=130$ input steps for the estimation, at the late stages of
training close to convergence, while Figure 9(b) considers the same for $l=90$
input steps. It is evident from Figure 9 that TREET captures relevant
information even from very long time series, while shorter length will lead
larger error of information extraction which cause a worse estimation result,
which is rational due to trimming of critical information.

$60$$70$$80$$90$TimeAmplitudeHeart
Rate$200$$220$$240$$260$$280$$300$$320$$340$$360$$380$$400$$-0.5$$0$$0.5$$1$$1.5$$\cdot
10^{4}$TimeAmplitudeBreath Rate

(a) Sampled sequence of Apnea dataset.

$3$$4$$5$$6$$7$$8$$9$$10$$11$$12$$13$$14$$15$$0$$0.5$$1$$1.5$$2$$2.5$$3$$3.5$$4$$4.5$$5$$5.5$$6$$\cdot
10^{-2}$Label Length $\displaystyle[l]$$\displaystyle TE_{X\to
Y}(k,l=2)$Breath $\rightarrow$ HeartHeart $\rightarrow$ Breath

(b) TE estimation for variable length history of $Y$.

Figure 10: Transfer Entropy estimation on physiological data. Apnea dataset
consists of heart rate and breath rate and we seek to find the information
flow for patient diagnosing. Patient with Apnea suffer from breathing
cessation which leads to alterations in heart rate during sleep.

Additionaly, we compare TREET with DINE for different input lengths on the
GMA(100) channel. Important to note that DINE can technically deal with long
sequences even if given a shorter backpropgation through time (bptt) input
length then the process memory is due to state propagation. However, as seen
in Table I, the DINE struggles to estimate the right capacity under long
memory, while TREET successfully estimates capacity with low error rate.
Nonetheless, when the TREET memory is shorter than the channel memory, its
performance significantly degrades.

| Memory | Estimated | Absolute
---|---|---|---
Model | Length | Capacity [nits] | Error (%)
DINE | 90 | 0.3435 | 15.3 (%)
TREET | 90 | 0.2962 | 26.9 (%)
DINE | 130 | 0.3333 | 17.8 (%)
TREET | 130 | 0.3851 |  5.0 (%)
TABLE I: GMA(100), with 0 SNR, capacity estimation with different lengths.
Note that memory length in DINE is the bptt length for LSTM, while in TREET it
is the sequence input length.

###  V-C Transfer Entropy Estimation in Physiological Data

Motivated by the results in [4, Chapter 7.1], we tested the TREET on the Apnea
dataset from Santa Fe Time Series
Competition333https://physionet.org/content/santa-fe/1.0.0/ [43, 44]. This is
a multivariate data that has been recorded from a diseased patient of Apnea in
a sleep laboratory, which is a condition characterized by brief, involuntary
pauses in breathing, particularly during sleep. Each sample consists of three
different variables from a specific time, with a sample rate of 2 Hz. The
three features are heart rate, chest volume (which is the respiration force)
and blood oxygen concentration. A sampled sequence of heart rate and breath
rate (chest volume) is presented in Figure 10(a).

Determining the interaction between different physiological features is
crucial for diagnosing diseases by revealing causal connections within the
human body, enabling targeted diagnostics and personalized treatments
according to identified risk factors. Therefore, we applied the TREET to
determine the magnitude and direction of information transfer in the given
setting. We used both heart rate and breath features and compared the TE
measurements, $\mathsf{TE}_{\mathsf{Breath}\to\mathsf{Heart}}(k,2)$ and
$\mathsf{TE}_{\mathsf{Heart}\to\mathsf{Breath}}(k,2)$ for variable length $k$
that is the $Y$ process’s history observations length. The results are
presented in Figure 10(b) and are aligned with the results of [4] (for
$k,l=2$) and extend it further with tests on variate history length of $Y$
process.

Notably, for every considered $k$, the TE from the breath process to the heart
process is consistently higher, aligning with the diagnosis of Apnea disorder
(abrupt cessation of breathing during sleep). In terms of information flow,
our results indicate that the breathing process transfers more information
regarding the behavior of the heart rate process, than the opposite direction,
since the value of the estimated TE is consistently higher in comparison,
while in information theory it essentially reflect that the $X$ process has
more to reveal about $Y$’s current value than $Y$’s past itself.

Furthermore, in the influence direction
$\mathsf{TE}_{\mathsf{Breath}\to\mathsf{Heart}}$, we can infer that increasing
the of visible heart rate history samples decreases the information transfer,
since a longer history of heart rate provides more insight into future
outcomes than the instantaneous breath value. However,
$\mathsf{TE}_{\mathsf{Heart}\to\mathsf{Breath}}$ shows that for variate $k$
values, the majority of TE results indicate that the heart process barely
affect the breathing process in Apnea patients.

The conclusion of this experiment is very valuable and validates what is known
to science - although it is known that when the heart rate increases the
breath rate increases also, for Apnea patients, the sudden cessation of muscle
movement during inhalation, i.e., the stopping of breathing, affects the heart
rate, in contrast to what is observed in healthy individuals. This experiment
supports the results of [4] and add a new insight about the information
transfer decreasing for increasing context of $Y$ process (which reveals more
information about $Y$, thus will lead to lower value of TE as mentioned
before).

##  VI Conclusions and Future Work

This work presented a modified attention-based architecture to estimate TE,
for a class of ergodic and stationary processes. We devise a DV-based neural
TE estimator, proved its consistency and described the proposed novel modified
attention mechanism for the task at hand. We then developed an optimizer of
estimated TE, which was then leveraged for the estimation of channel capacity.
We concluded by studying the TREET application in causal features analysis on
the Apnea dataset and compared its performance with RNN-based mechanisms.

With the increasing popularity of sequential data to most contemporary fields
of machine learning, we plan to leverage the TREET for information theoretic
analysis and design of architectures, through the lenses of causal information
transfer. Examples of such applications include enhancing predictive models,
refining feature selection processes, reconstructing complex networks,
improving anomaly detection capabilities, and optimizing decision-making in
dynamic environments. Furthermore, building upon the work in [45], we plan to
extend the TREET optimization scheme to sequential data compression tasks.

##  VII Appendix

###  VII-A Lemma TE Equals to DI

###### Lemma 3

Define markov property $Z_{n}-Z_{n-1}-Z^{n-2}$ as
$P_{Z_{n}|Z^{n-1}}=P_{Z_{n}|Z_{n-1}}$. Let $\mathbb{X}$ and $\mathbb{Y}$ be
two jointly stationary processes such that the following markov property
holds,

$\displaystyle Y_{l}-Y^{l-1}-Y^{0}_{-\infty},$ $\displaystyle
Y_{l}-(X^{l},Y^{l-1})-(X^{0}_{-\infty},Y^{0}_{-\infty}),$

for $l\in\mathbb{N}$. Then,

$\mathsf{TE}_{X\to Y}(m)=\mathsf{TE}_{X\to Y}(l)$ (31)

for $m\geq l$ and

$\mathsf{TE}_{X\to Y}(m)=\mathsf{I}\left(\mathbb{X}\to\mathbb{Y}\right).$ (32)

proof: Let $\mathbb{X}$, $\mathbb{Y}$ be two jointly stationary processes that
pose the markov property with $l\in\mathbb{N}$ and $m\geq l$, then

$\displaystyle\mathsf{TE}_{X\to Y}(l)$
$\displaystyle=\mathsf{I}\left(X^{l};Y_{l}|Y^{l-1}\right)$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{=}}\mathsf{I}\left(X^{l}_{l-m};Y_{l}|Y^{l-1}_{l-m}\right)$
$\displaystyle\stackrel{{\scriptstyle(b)}}{{=}}\mathsf{I}\left(X^{m};Y_{m}|Y^{m-1}\right)$
$\displaystyle=\mathsf{TE}_{X\to Y}(m),$ (33)

where (a) is true from the markovity, and (b) transition is index shift, and
is valid due to stationarity of the process. Observing the DI rate,

$\displaystyle\mathsf{I}(\left(\mathbb{X}\to\mathbb{Y}\right)$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{=}}\lim_{n\to\infty}\frac{1}{n}\sum_{i=1}^{n}\mathsf{I}\left(X^{i};Y_{i}|Y^{i-1}\right)$
$\displaystyle=\lim_{n\to\infty}\frac{1}{n}\sum_{i=1}^{n}\left[\mathsf{h}\left(Y_{i}|Y^{i-1}\right)-\mathsf{h}\left(Y_{i}|X^{i},Y^{i-1}\right)\right]$
$\displaystyle\stackrel{{\scriptstyle(b)}}{{=}}\lim_{n\to\infty}\mathsf{h}\left(Y_{n}|Y^{n-1}\right)-\mathsf{h}\left(Y_{n}|X^{n},Y^{n-1}\right)$
$\displaystyle=\lim_{n\to\infty}\mathsf{I}\left(X^{n};Y_{n}|Y^{n-1}\right)$
$\displaystyle\stackrel{{\scriptstyle(c)}}{{=}}\lim_{n\to\infty}\mathsf{TE}_{X\to
Y}(n)$ $\displaystyle\stackrel{{\scriptstyle(d)}}{{=}}\mathsf{TE}_{X\to
Y}(m),$ (34)

where the limit in (a) exists whenever the joint process is stationary,
transition (b) is valid because the limit exists for each series of
conditional entropies, and since conditioning reduces entropy, the limit for
each normalized sum of series is
$\lim_{n\to\infty}\mathsf{h}\left(Y_{n}|Y^{n-1}\right),\lim_{n\to\infty}\mathsf{h}\left(Y_{n}|X^{n},Y^{n-1}\right)$,
respectively [40, Theorem 4.2.1]. Since the limit exists for the conditional
MI, transition (c) is valid by definition of TE, and the TE with limit of
parameter $m$ exists, and (d) is from (VII-A) for $n\geq m$. Concluding the
proof. $\square$

###  VII-B Proof of theorem 3

The proof following the steps of representation step - represents TE as a
subtraction of two DV potentials, estimation step - proves that the DV
potentials is achievable by empirical mean of a given set of samples, and
approximation step - shows that the estimator converges to TE with the
corresponding memory parameter $l$. Thus concluding that our estimator is a
consistent estimator for the TE.

Let $\\{X_{i},Y_{i}\\}_{i\in\mathbb{Z}}$ be the two values of a process
$\mathbb{X},\mathbb{Y}$ respectively, and $\mathbb{P}$ be the stationary
ergodic measure over $\sigma(\mathbb{X},\mathbb{Y})$. Define
$P_{X^{n},Y^{n}}:=\mathbb{P}|_{\sigma(X^{n},Y^{n})}$ as the $n$-coordinate
projection of $\mathbb{P}$, where $\sigma(X^{n},Y^{n})$ is the
$\sigma$-algebra generated by $(X^{n},Y^{n})$. Let $D_{n}=(X^{n},Y^{n})\sim
P_{X^{n},Y^{n}}$. Let $\widetilde{Y}\sim{Q}_{Y}$ be an absolutely continuous
PDF, independent of $\\{(X_{i},Y_{i})\\}_{i\in\mathbb{Z}}$ and its
distribution noted as $\widetilde{P}_{Y}$. The proof is divided to three steps
- variational representation, estimation from samples and functional
approximation.

####  VII-B1 Representation of TE

In order to write TE as a difference between two KL divergences, recall Lemma
1 \- let,

$\displaystyle\mathsf{D}_{Y_{l}|Y^{l-1}\|\widetilde{Y}_{l}}:=\mathsf{D}_{\mathsf{KL}}\left(P_{Y_{l}|Y^{l-1}}\|\widetilde{P}_{Y_{l}}\Big{|}P_{Y^{l-1}}\right),$
(35a)
$\displaystyle\mathsf{D}_{Y_{l}|Y^{l-1}X^{l}\|\widetilde{Y}_{l}}:=\mathsf{D}_{\mathsf{KL}}\left(P_{Y_{l}|Y^{l-1}X^{l}}\|\widetilde{P}_{Y_{l}}\Big{|}P_{Y^{l-1}X^{l}}\right).$
(35b)

Then

$\mathsf{TE}_{X\to
Y}(l)=\mathsf{D}_{Y_{l}|Y^{l-1}X^{l}\|\widetilde{Y}_{l}}-\mathsf{D}_{Y_{l}|Y^{l-1}\|\widetilde{Y}_{l}}.$
(36)

This lemma is proved in Section VII-C. Estimating the KL divergence is
applicable with the DV representation 2

$\displaystyle\mathsf{D}_{Y_{l}|Y^{l-1}\|\widetilde{Y}_{l}}=$
$\displaystyle\sup_{f_{y}:\Omega_{\mathcal{Y}}\to\mathbb{R}}\mathbb{E}\left[f_{y}\left(Y^{l}\right)\right]$
$\displaystyle-\log\mathbb{E}\left[e^{f_{y}\left(Y^{l-1},\widetilde{Y}_{l}\right)}\right],$
(37a) where $\Omega_{\mathcal{Y}}=\mathcal{Y}^{l}$. For the other term, the DV
representation is,
$\displaystyle\mathsf{D}_{Y_{l}|Y^{l-1}X^{l}\|\widetilde{Y}_{l}}=$
$\displaystyle\sup_{f_{xy}:\Omega_{\mathcal{X}\times\mathcal{Y}}\to\mathbb{R}}\mathbb{E}\left[f_{xy}\left(X^{l},Y^{l}\right)\right]$
$\displaystyle-\log\mathbb{E}\left[e^{f_{xy}\left(X^{l},Y^{l-1},\widetilde{Y}_{l}\right)}\right],$
(37b)

where
$\Omega_{\mathcal{X}\times\mathcal{Y}}=\mathcal{X}^{l}\times\mathcal{Y}^{l}$.
The next section refers to (37a) but the claims are the same for (37b).

####  VII-B2 Estimation

By the DV representation, the supremum in (37a) achieved for

$\displaystyle f_{y,l}^{\star}$
$\displaystyle:=\log\left(\frac{dP_{Y^{l}}}{d(P_{Y^{l-1}}\otimes\widetilde{P}_{Y_{l}})}\right)$
$\displaystyle=\log p_{Y_{l}|Y^{l-1}}-\log\widetilde{p}_{Y_{l}},$ (38)

where the last equality holds due to $P_{Y^{l}}\ll
P_{Y^{l-1}}\otimes\widetilde{P}_{Y_{l}}$, both measures have Lebesgue
densities. While it is not mandatory to select the reference measurement as
uniform, choosing a uniform reference measurement can result in a constant
that can be neutralized to obtain the likelihood function of $Y_{l}|Y^{l-1}$.
This approach allows for simplification and facilitates the estimation
process. Empirical means can estimate the expectations in (37a), while
applying the generalized Birkhoff theorem [46], stated next:

######  Theorem 4 (The generalized Birkhoff theorem)

Let $T$ be a metrically transitive 1 - 1 measure preserving transformation of
the probability space $(\Omega,\mathcal{F},\mathbb{P})$ onto itself. Let
$g_{0}(\omega),g_{1}(\omega),\ldots$ be a sequence of measurable functions on
$\Omega$ converging a.s. to the function $g(\omega)$ such that
$\mathbb{E}[\sup_{i}|g_{i}|]\leq\infty$. Then,

$\frac{1}{n}\sum_{i=1}^{n}g_{i}(T^{i}\omega)\stackrel{{\scriptstyle
n\to\infty}}{{\longrightarrow{}}}\mathbb{E}[g],\quad\mathbb{P}-a.s.$ (39)

By applying Theorem 4 and for any $\epsilon>0$ and sufficiently large $n$, we
have

$\displaystyle\left|\mathbb{E}\left[f_{y,l}^{\star}\left(Y^{l}\right)\right]-\frac{1}{n}\sum_{i=0}^{n-1}f_{y,l}^{\star}\left(Y_{i}^{i+l}\right)\right|<\frac{\epsilon}{8},$
(40a)
$\displaystyle\left|\log\left(\mathbb{E}\left[e^{f_{y,l}^{\star}\left(Y^{l-1},\widetilde{Y}_{l}\right)}\right]\right)\right.$
$\displaystyle\quad\left.-\log\left(\frac{1}{n}\sum_{i=0}^{n-1}e^{f_{y,l}^{\star}\left(Y_{i}^{i+l-1},\widetilde{Y}_{i+l}\right)}\right)\right|<\frac{\epsilon}{8},$
(40b)

where $\\{f_{y,l}^{\star}\\}$ is the function of $l$ time steps, that achieves
the supremum of
$\mathsf{D}_{\mathsf{KL}}\left(P_{Y_{l}|Y^{l-1}}\|\widetilde{P}_{Y_{l}}\Big{|}P_{Y^{l-1}}\right)$.
Convergence achieved from the generalized Brikhoff theorem, where the series
of functions is the fixed function $f_{y,l}^{\star}$,

$\displaystyle\mathbb{E}_{n}\left[f_{y,l}^{\star}\left(Y^{l}\right)\right]:=\frac{1}{n}\sum_{i=0}^{n-1}f_{y,l}^{\star}\left(Y_{i}^{i+l}\right),$
(41a)
$\displaystyle\mathbb{E}_{n}\left[e^{f_{y,l}^{\star}\left(Y^{l-1},\widetilde{Y}_{l}\right)}\right]:=\frac{1}{n}\sum_{i=0}^{n-1}e^{f_{y,l}^{\star}\left(Y_{i}^{i+l-1},\widetilde{Y}_{i+l}\right)}.$
(41b)

####  VII-B3 Approximation

Last step is to approximate the functional space with the space of
transformers. Recall that set of causal transformer architectures
$\mathcal{G}_{\mathsf{ctf}}^{Y}:=\mathcal{G}_{\mathsf{ctf}}^{(d_{y},1,l,v_{y})},\mathcal{G}_{\mathsf{ctf}}^{XY}:=\mathcal{G}_{\mathsf{ctf}}^{(d_{x}+d_{y},1,l,v_{xy})}$
for given $l,d_{y},d_{x},v_{y},v_{xy}\in\mathbb{N}$. Define

$\displaystyle\widehat{\mathsf{D}}_{Y_{l}|Y^{l-1}\|\widetilde{Y}_{l}}(D_{n}):=\sup_{g_{y}\in\mathcal{G}_{\mathsf{tf}}^{Y}}\frac{1}{n}\sum_{i=0}^{n-1}g_{y}\left(Y_{i}^{i+l}\right)$
$\displaystyle\qquad\qquad\quad-\log\left(\frac{1}{n}\sum_{i=0}^{n-1}e^{g_{y}\left(Y_{i}^{i+l-1},\widetilde{Y}_{i+l}\right)}\right),$
(42)

where the DV functions are transformers
$g_{y}\in\mathcal{G}_{\mathsf{ctf}}^{Y},g_{xy}\in\mathcal{G}_{\mathsf{ctf}}^{XY}$.
We want to prove that for a given $\epsilon>0$

$\left|\widehat{\mathsf{D}}_{Y_{l}|Y^{l-1}\|\widetilde{Y}_{l}}(D_{n})-\mathsf{D}_{Y_{l}|Y^{l-1}\|\widetilde{Y}_{l}}\right|\leq\frac{\epsilon}{2}.$
(43)

From Theorem 2, we obtain

$\displaystyle\mathbb{E}\left[f_{y,l}^{\star}\left(Y^{l}\right)\right]=\mathsf{D}_{Y_{l}|Y^{l-1}\|\widetilde{Y}_{l}},$
(44a)
$\displaystyle\mathbb{E}\left[f_{y,l}^{\star}\left(Y^{l-1},\widetilde{Y}_{l}\right)\right]=1.$
(44b)

Thus, we bound the term following
$\left|\widehat{\mathsf{D}}_{Y_{l}|Y^{l-1}\|\widetilde{Y}_{l}}(D_{n})-\mathbb{E}\left[f_{y,l}^{\star}\left(Y_{i}^{i+l}\right)\right]\right|$.
By the identity $\log(x)\leq x-1,\forall x\in\mathbb{R}_{\geq 0}$ we gain

$\displaystyle\left|\widehat{\mathsf{D}}_{Y_{l}|Y^{l-1}\|\widetilde{Y}_{l}}(D_{n})-\mathbb{E}\left[f_{y,l}^{\star}\left(Y^{l}\right)\right]\right|$
$\displaystyle=\left|-\mathbb{E}\left[f_{y,l}^{\star}\left(Y^{l}\right)\right]+\sup_{g_{y}\in\mathcal{G}_{\mathsf{tf}}^{Y}}\left\\{\frac{1}{n}\sum_{i=0}^{n-1}g_{y}\left(Y_{i}^{i+l}\right)\right.\right.$
$\displaystyle\qquad\qquad\qquad\left.\left.-log\left(\frac{1}{n}\sum_{i=0}^{n-1}e^{g_{y}\left(Y_{i}^{i+l-1},\widetilde{Y}_{i+l}\right)}\right)\right\\}\right|$
$\displaystyle\leq\left|1-\mathbb{E}\left[f_{y,l}^{\star}\left(Y^{l}\right)\right]+\sup_{g_{y}\in\mathcal{G}_{\mathsf{tf}}^{Y}}\left\\{\frac{1}{n}\sum_{i=0}^{n-1}g_{y}\left(Y_{i}^{i+l}\right)\right.\right.$
$\displaystyle\qquad\qquad\qquad\quad\left.\left.-\left(\frac{1}{n}\sum_{i=0}^{n-1}e^{g_{y}\left(Y_{i}^{i+l-1},\widetilde{Y}_{i+l}\right)}\right)\right\\}\right|$
$\displaystyle\leq\left|+\mathbb{E}\left[f_{y,l}^{\star}\left(Y^{l-1},\widetilde{Y}_{l}\right)\right]-\mathbb{E}\left[f_{y,l}^{\star}\left(Y^{l}\right)\right]\right.$
$\displaystyle\quad\left.+\sup_{g_{y}\in\mathcal{G}_{\mathsf{tf}}^{Y}}\left\\{\frac{1}{n}\sum_{i=0}^{n-1}g_{y}\left(Y_{i}^{i+l}\right)\right.\right.$
$\displaystyle\qquad\qquad\qquad\quad\left.\left.-\left(\frac{1}{n}\sum_{i=0}^{n-1}e^{g_{y}\left(Y_{i}^{i+l-1},\widetilde{Y}_{i+l}\right)}\right)\right\\}\right|.$
(45)

Due to (40), there exists $N\in\mathbb{N}$ such that $\forall n>\mathbb{N}$

$\displaystyle\left|\mathbb{E}\left[f_{y,l}^{\star}\left(Y^{l}\right)\right]-\mathbb{E}_{n}\left[f_{y,l}^{\star}\left(Y^{l}\right)\right]\right|<\frac{\epsilon}{8},$
(46a)
$\displaystyle\left|\left(\mathbb{E}\left[e^{f_{y,l}^{\star}\left(Y^{l-1},\widetilde{Y}_{l}\right)}\right]\right)-\mathbb{E}_{n}\left[e^{f_{y,l}^{\star}\left(Y^{l-1},\widetilde{Y}_{l}\right)}\right]\right|<\frac{\epsilon}{8}.$
(46b)

Plugging (46) to (45) gives

$\displaystyle\left|\widehat{\mathsf{D}}_{Y_{l}|Y^{l-1}\|\widetilde{Y}_{l}}(D_{n})-\mathsf{D}_{Y_{l}|Y^{l-1}\|\widetilde{Y}_{l}}\right|$
$\displaystyle\leq\frac{\epsilon}{4}+\Bigg{|}-\mathbb{E}_{n}\left[e^{f_{y,l}^{\star}(Y^{l-1},\widetilde{Y}_{l})}\right]-\mathbb{E}_{n}\left[f_{y,l}^{\star}\left(Y^{l}\right)\right]$
$\displaystyle\quad+\sup_{g_{y}\in\mathcal{G}_{\mathsf{tf}}^{Y}}\left\\{\frac{1}{n}\sum_{i=0}^{n-1}g_{y}\left(Y_{i}^{i+l}\right)\right.$
$\displaystyle\qquad\qquad\qquad\left.-\left(\frac{1}{n}\sum_{i=0}^{n-1}e^{g_{y}\left(Y_{i}^{i+l-1},\widetilde{Y}_{i+l}\right)}\right)\right\\}\Bigg{|}.$
(47)

Since the empirical mean of $f_{y,l}^{\star}$ is converging to the expected
mean, it is uniformly bounded by some $M\in\mathbb{R}_{\geq 0}$. Since the
exponent function is Lipschitz continuous with Lipschitz constant $\exp{[M]}$
on the interval $(-\infty,M]$, we obtain

$\displaystyle\frac{1}{n}\sum_{i=1}^{n}e^{f_{y,l}^{\star}\left(\widetilde{Y}_{i+l},Y_{i}^{i+l-1}\right)}-e^{g_{y}\left(\widetilde{Y}_{l},Y_{i}^{i+l-1}\right)}$
$\displaystyle\leq
e^{M}\frac{1}{n}\sum_{i=1}^{n}\Big{|}f_{y,l}^{\star}\left(\widetilde{Y}_{i+l},Y_{i}^{i+l-1}\right)$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad-
g_{y}\left(\widetilde{Y}_{i+l},Y_{i}^{i+l-1}\right)\Big{|}.$ (48)

Definition 4 is a sub-functions class of the continuous sequence to sequence
functions class, and applies to the causal transformer. Thus, concludes that
the causal transformer $g\in\mathcal{G}_{\mathsf{ctf}}^{(d_{i},d_{o},l,v)}$
also applies for the universal approximation theorem, for continuous vector-
values functions. Moreover, for our case the output sequence is a scalar
value, $\mathcal{U}\subset\mathbb{R}^{l\times
d_{i}},\mathcal{Z}\subset\mathbb{R}^{1\times d_{o}}$. For given $\epsilon,M,l$
and $n$, denote $g_{y}^{\star}\in\mathcal{G}_{\mathsf{ctf}}^{(d_{y},1,l,v)}$
the Causal transformer, such that the approximation error is uniformly bounded
$\exp{[-M]}\times{\epsilon}/{4}$ for the final time prediction out from the
model. Finally, combining Theorem 1 we have

$\displaystyle\left|\widehat{\mathsf{D}}_{Y_{l}|Y^{l-1}\|\widetilde{Y}_{l}}(D_{n})-\mathsf{D}_{Y_{l}|Y^{l-1}\|\widetilde{Y}_{l}}\right|$
$\displaystyle\leq\left(1+e^{M}\right)\frac{1}{n}\sum_{i=1}^{n}\Big{|}f_{y,l}^{\star}\left(\widetilde{Y}_{i+l},Y_{i}^{i+l-1}\right)$
$\displaystyle\quad\qquad\qquad\qquad\qquad-
g_{y}^{\star}\left(\widetilde{Y}_{i+l},Y_{i}^{i+l-1}\right)\Big{|}+\frac{\epsilon}{4}$
$\displaystyle\leq\frac{\epsilon}{2}.$ (49)

This concludes the proof of (37a). For (37b), note that

$\displaystyle f_{y,l}^{\star}$
$\displaystyle:=\log\left(\frac{dP_{Y^{l}}}{d(P_{X^{l}Y^{l-1}}\otimes\widetilde{P}_{Y_{l}})}\right)$
$\displaystyle=\log p_{Y_{l}|X^{l}Y^{l-1}}-\log\widetilde{p}_{Y_{l}},$ (50)

achieves the supremum. Following the same claims for
$\widehat{\mathsf{D}}_{Y_{l}|Y^{l-1}X^{l}\|\widetilde{Y}_{l}}(D_{n})$

$\left|\widehat{\mathsf{D}}_{Y_{l}|Y^{l-1}X^{l}\|\widetilde{Y}_{l}}(D_{n})-\mathsf{D}_{Y_{l}|Y^{l-1}X^{l}\|\widetilde{Y}_{l}}\right|\leq\frac{\epsilon}{2},$
(51)

where

$\displaystyle\widehat{\mathsf{D}}_{Y_{l}|Y^{l-1}X^{l}\|\widetilde{Y}_{l}}(D_{n})$
$\displaystyle:=\sup_{g_{xy}\in\mathcal{G}_{\mathsf{tf}}^{XY}}\frac{1}{n}\sum_{i=0}^{n-1}g_{xy}\left(Y_{i}^{i+l},X_{i}^{i+l}\right)$
$\displaystyle\qquad\qquad-\log\left(\frac{1}{n}\sum_{i=0}^{n-1}e^{g_{xy}\left(Y_{i}^{i+l-1},X_{i}^{i+l},\widetilde{Y}_{i+l}\right)}\right).$
(52)

With (51) and (49) we end the proof. $\blacksquare$

###  VII-C Proof of lemma 1

Recall

$\displaystyle\mathsf{D}_{Y_{l}|Y^{l-1}\|\widetilde{Y}_{l}}:=\mathsf{D}_{\mathsf{KL}}\left(P_{Y_{l}|Y^{l-1}}\|\widetilde{P}_{Y_{l}}\Big{|}P_{Y^{l-1}}\right),$
(53a)
$\displaystyle\mathsf{D}_{Y_{l}|Y^{l-1}X^{l}\|\widetilde{Y}_{l}}:=\mathsf{D}_{\mathsf{KL}}\left(P_{Y_{l}|Y^{l-1}X^{l}}\|\widetilde{P}_{Y_{l}}\Big{|}P_{Y^{l-1}X^{l}}\right).$
(53b)

By expanding the first term we obtain,

$\displaystyle\mathsf{D}_{Y_{l}|Y^{l-1}\|\widetilde{Y}_{l}}$
$\displaystyle=\mathbb{E}_{P_{Y^{l-1}}}\left[\mathsf{D}_{\mathsf{KL}}(P_{Y_{l}|Y^{l-1}}\|P_{\widetilde{Y}_{l}})\right]$
$\displaystyle=\int_{\mathcal{Y}^{l-1}}\Bigg{[}\int_{\mathcal{Y}_{l}}\log\left(\frac{P(y_{l}|y^{l-1})}{\widetilde{P}(y^{l})}\right)$
$\displaystyle\qquad\qquad\qquad
P(y_{l}|y^{l-1})d(y_{l})\Bigg{]}P(y^{l-1})d(y^{l-1})$
$\displaystyle=\int_{\mathcal{Y}^{l}}\log\left(\frac{P(y_{l}|y^{l-1})}{\widetilde{P}(y_{l})}\right)P(y^{l})d(y^{l})$
$\displaystyle=\int_{\mathcal{X}^{l}\mathcal{Y}^{l}}\log\left(\frac{P(y_{l}|y^{l-1})}{\widetilde{P}(y_{l})}\right)P(x^{l}y^{l})d(x^{l}y^{l}),$
(54)

and the second term,

$\displaystyle\mathsf{D}_{Y_{l}|Y^{l-1}X^{l}\|\widetilde{Y}_{l}}$
$\displaystyle=\mathbb{E}_{P_{X^{l}Y^{l-1}}}\left[\mathsf{D}_{\mathsf{KL}}(P_{Y_{l}|X^{l}Y^{l-1}}\|P_{\widetilde{Y}_{l}})\right]$
$\displaystyle=\int_{\mathcal{X}^{l}\mathcal{Y}^{l-1}}\Bigg{[}\int_{\mathcal{Y}_{l}}\log\left(\frac{P(y_{l}|x^{l}y^{l-1})}{\widetilde{P}(y_{l})}\right)$
$\displaystyle\qquad\qquad\quad
P(y_{l}|x^{l}y^{l-1})d(y_{l})\Bigg{]}P(x^{l}y^{l-1})d(x^{l}y^{l-1})$
$\displaystyle=\int_{\mathcal{X}^{l}\mathcal{Y}^{l}}\log\left(\frac{P(y_{l}|x^{l}y^{l-1}}{\widetilde{P}(y_{l})}\right)P(x^{l}y^{l})d(x^{l}y^{l}),$
(55)

where $\widetilde{P}_{Y}$ is a reference density function and is absolutely
continuous on $\mathcal{Y}$. Subtructing the two terms resulting,

$\displaystyle\mathsf{D}_{Y_{l}|Y^{l-1}X^{l}\|\widetilde{Y}_{l}}-\mathsf{D}_{Y_{l}|Y^{l-1}\|\widetilde{Y}_{l}}$
$\displaystyle=\int_{\mathcal{X}^{l}\mathcal{Y}^{l}}\log\left(\frac{P(y_{l}|x^{l}y^{l-1})}{P(y_{l}|y^{l-1}}\right)P(x^{l}y^{l})d(x^{l}y^{l})$
$\displaystyle=h(Y_{l}|Y^{l-1})-h(Y_{l}|X^{l}Y^{l-1})$
$\displaystyle=\mathsf{TE}_{X\to Y}(l).\quad\square$ (56)

###  VII-D Proof of lemma 2

The optimal TE, $\mathsf{TE}_{X\to Y}^{\star}(l)$, by non-finite memory
process $\mathbb{X}$, is given by

$\mathsf{TE}_{X\to
Y}^{\star}(l):=\lim_{n\to\infty}\sup_{P_{X^{l}_{-n}}}\mathsf{TE}_{X\to Y}(l).$
(57)

For any $n>0$ we have,

$\displaystyle\sup_{P_{X^{l}_{-n}}}\mathsf{TE}_{X\to Y}(l)$
$\displaystyle=\sup_{P_{X^{0}_{-n}}P_{X^{l}|X^{0}_{-n}}}\mathsf{I}\left(X^{l};Y_{l}|Y^{l-1}\right)$
$\displaystyle=\sup_{P_{X^{0}_{-n}}P_{X^{l}|X^{0}_{-n}}}\int_{\mathcal{X}^{l}_{-n}\mathcal{Y}^{l}}P(x^{0}_{-n})P(x^{l},y^{l}|x^{0}_{-n})$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\underbrace{\mathsf{I}\left(X^{l};Y_{l}|Y^{l-1}\right)}_{\text{independent
of }x^{0}_{-n}}d(x^{l}x^{0}_{-n}y^{l})$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{=}}\sup_{P_{X^{l}}}\int_{\mathcal{X}^{l}\mathcal{Y}^{l}}P(x^{l},y^{l})\mathsf{I}\left(X^{l};Y_{l}|Y^{l-1}\right)d(x^{l}y^{l})$
$\displaystyle=\sup_{P_{X^{l}}}\mathsf{TE}_{X\to Y}(l),$ (58)

where (a) follows from the fact that the conditional MI does not depend on
$x^{0}_{-n}$, thus we can eliminate the integral of $\mathcal{X}_{-n}^{0}$
that does not affect the conditional MI, as for the supremum. Since (VII-D) is
true for any $n\in\mathbb{N}$, with (57) we get,

$\mathsf{TE}_{X\to Y}^{\star}(l)=\sup_{P_{X^{l}}}\mathsf{TE}_{X\to
Y}(l).\quad\square$ (59)

###  VII-E Implementation Parameters

For the channel capacity estimation model, we trained the optimization
procedure with a limit of 200 epochs. We utilized a batch size of 1024, a
learning rate of $8\times 10^{-3}$, and the Adam optimizer [47]. The
transformer architecture comprises one attention layer followed by a FF layer.
The attention mechanism is implemented with a single head, having a dimension
of 32 neurons. The dimension of the FF layer is 64, featuring ELU activation
[48]. The input sequence length is set to 59, and we back-propagate through
sequences of length 30, corresponding to the parameter $l$ in TE. Each epoch
generates 100K samples of channel input-output tuples, with random noise
uniformly distributed over the bounds of a given batch. The noise distribution
generator (NDG) also employs a transformer structure with one attention layer
and one FF layer, maintaining the same dimensional specifications. The
learning rate for the NDG is set to $8\times 10^{-4}$. It is trained at a rate
of once every four epochs.

For the Apnea dataset feature analysis, the learning rate was adjusted to
$1\times 10^{-4}$. Additionally, we configured the model to use 1 head with a
dimension of 16 for the attention mechanism, and the FF layer dimension was
set to 32, with ELU activation. All experiments were conducted on Nvidia RTX
3090 GPUs.

## References

  * [1] Thomas Schreiber.  Measuring information transfer.  Physical review letters, 85(2):461, 2000.
  * [2] JM Nichols.  Examining structural dynamics using information flow.  Probabilistic Engineering Mechanics, 21(4):420–433, 2006.
  * [3] Ping Duan, Fan Yang, Tongwen Chen, and Sirish L Shah.  Direct causality detection via the transfer entropy approach.  IEEE transactions on control systems technology, 21(6):2052–2066, 2013.
  * [4] T. R. J. Bossomaier, Lionel C. Barnett, Michael S. Harré, and Joseph T. Lizier.  An introduction to transfer entropy.  In Cambridge International Law Journal, 2016.
  * [5] Chulhee Yun, Srinadh Bhojanapalli, Ankit Singh Rawat, Sashank J Reddi, and Sanjiv Kumar.  Are transformers universal approximators of sequence-to-sequence functions?  arXiv preprint arXiv:1912.10077, 2019.
  * [6] Haixu Wu, Jiehui Xu, Jianmin Wang, and Mingsheng Long.  Autoformer: Decomposition transformers with auto-correlation for long-term series forecasting.  Advances in Neural Information Processing Systems, 34:22419–22430, 2021.
  * [7] Haoyi Zhou, Shanghang Zhang, Jieqi Peng, Shuai Zhang, Jianxin Li, Hui Xiong, and Wancai Zhang.  Informer: Beyond efficient transformer for long sequence time-series forecasting.  In Proceedings of the AAAI conference on artificial intelligence, volume 35, pages 11106–11115, 2021.
  * [8] Tian Zhou, Ziqing Ma, Qingsong Wen, Xue Wang, Liang Sun, and Rong Jin.  Fedformer: Frequency enhanced decomposed transformer for long-term series forecasting.  In International Conference on Machine Learning, pages 27268–27286. PMLR, 2022.
  * [9] Qingsong Wen, Tian Zhou, Chaoli Zhang, Weiqi Chen, Ziqing Ma, Junchi Yan, and Liang Sun.  Transformers in time series: A survey, 2023.
  * [10] Boris Gourévitch and Jos J. Eggermont.  Evaluating information transfer between auditory cortical neurons.  Journal of Neurophysiology, 97(3):2533–2543, 2007.  PMID: 17202243.
  * [11] Shinya Ito, Michael E Hansen, Randy Heiland, Andrew Lumsdaine, Alan M Litke, and John M Beggs.  Extending transfer entropy improves identification of effective connectivity in a spiking cortical network model.  PloS one, 6(11):e27431, 2011.
  * [12] M. D. Madulara, P. A. B. Francisco, S. Nawang, D. C. Arogancia, C. J. Cellucci, P. E. Rapp, and A. M. Albano.  Eeg transfer entropy tracks changes in information transfer on the onset of vision.  International Journal of Modern Physics: Conference Series, 17:9–18, 2012.
  * [13] Max Lungarella and Olaf Sporns.  Mapping information flow in sensorimotor networks.  PLoS computational biology, 2(10):e144, 2006.
  * [14] Greg Ver Steeg and Aram Galstyan.  Information transfer in social media.  In Proceedings of the 21st international conference on World Wide Web, pages 509–518, 2012.
  * [15] Sondo Kim, Seungmo Ku, Woojin Chang, and Jae Wook Song.  Predicting the direction of us stock prices using effective transfer entropy and machine learning techniques.  IEEE Access, 8:111660–111682, 2020.
  * [16] Paolo Bonetti, Alberto Maria Metelli, and Marcello Restelli.  Causal feature selection via transfer entropy.  arXiv preprint arXiv:2310.11059, 2023.
  * [17] Murray Rosenblatt.  Remarks on some nonparametric estimates of a density function.  The annals of mathematical statistics, pages 832–837, 1956.
  * [18] Lyudmyla F Kozachenko and Nikolai N Leonenko.  Sample estimate of the entropy of a random vector.  Problemy Peredachi Informatsii, 23(2):9–16, 1987.
  * [19] Alexander Kraskov, Harald Stögbauer, and Peter Grassberger.  Estimating mutual information.  Physical review E, 69(6):066138, 2004.
  * [20] Clive WJ Granger.  Investigating causal relations by econometric models and cross-spectral methods.  Econometrica: journal of the Econometric Society, pages 424–438, 1969.
  * [21] Lionel Barnett and Terry Bossomaier.  Transfer entropy as a log-likelihood ratio.  Physical review letters, 109(13):138105, 2012.
  * [22] Lionel Barnett, Adam B Barrett, and Anil K Seth.  Granger causality and transfer entropy are equivalent for gaussian variables.  Physical review letters, 103(23):238701, 2009.
  * [23] M. D. Donsker and S. R. S. Varadhan.  Asymptotic evaluation of certain Markov process expectations for large time. iv.  Communications on Pure and Applied Mathematics, 36(2):183–212, 1983\.
  * [24] M. I. Belghazi et. al.  Mutual information neural estimation.  In International Conference on Machine Learning, pages 531–540. PMLR, 2018.
  * [25] Dor Tsur, Ziv Aharoni, Ziv Goldfeld, and Haim Permuter.  Neural estimation and optimization of directed information over continuous spaces.  IEEE Transactions on Information Theory, 2023.
  * [26] J. Zhang, O. Simeone, Z. Cvetkovic, E. Abela, and M. Richardson.  Itene: Intrinsic transfer entropy neural estimator.  arXiv preprint arXiv:1912.07277, 2019.
  * [27] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Łukasz Kaiser, and Illia Polosukhin.  Attention is all you need.  Advances in neural information processing systems, 30, 2017.
  * [28] Jiaming Song and Stefano Ermon.  Understanding the limitations of variational mutual information estimators.  arXiv preprint arXiv:1910.06222, 2019.
  * [29] David McAllester and Karl Stratos.  Formal limitations on the measurement of mutual information.  In International Conference on Artificial Intelligence and Statistics, pages 875–884. PMLR, 2020.
  * [30] Kwanghee Choi and Siyeong Lee.  Combating the instability of mutual information-based losses via regularization.  In Uncertainty in Artificial Intelligence, pages 411–421. PMLR, 2022.
  * [31] Ying Liu and Selin Aviyente.  The relationship between transfer entropy and directed information.  In 2012 IEEE Statistical Signal Processing Workshop (SSP), pages 73–76. IEEE, 2012.
  * [32] J. Massey.  Causality, feedback and directed information.  In Proc. Int. Symp. Inf. Theory Applic.(ISITA-90), pages 303–305. Citeseer, 1990.
  * [33] G. Kramer.  Directed information for channels with feedback, volume 11.  Citeseer, 1998.
  * [34] K. Hornik, M. Stinchcombe, and H. White.  Multilayer feedforward networks are universal approximators.  Neural networks, 2(5):359–366, 1989.
  * [35] S. Molavipour, H. Ghourchian, G. Bassi, and M. Skoglund.  Neural estimator of information for time-series data with dependency.  Entropy, 23(6):641, 2021.
  * [36] Ke Bai, Pengyu Cheng, Weituo Hao, Ricardo Henao, and Larry Carin.  Estimating total correlation with mutual information estimators.  In International Conference on Artificial Intelligence and Statistics, pages 2147–2164. PMLR, 2023.
  * [37] Roland L. D.  General formulation of Shannon’s main theorem in information theory.  American mathematical society translations, 33:323–438, 1963.
  * [38] Sekhar Tatikonda and Sanjoy Mitter.  The capacity of channels with feedback.  IEEE Transactions on Information Theory, 55(1):323–349, 2008.
  * [39] A. El Gamal and Y. H. Kim.  Network information theory.  Cambridge university press, 2011.
  * [40] T. M. Cover and J. A. Thomas.  Elements of Information Theory.  Wiley, New-York, 2nd edition, 2006.
  * [41] S. Yang, A. Kavcic, and S. Tatikonda.  On the feedback capacity of power-constrained gaussian noise channels with memory.  IEEE Transactions on Information Theory, 53(3):929–954, 2007.
  * [42] O. Sabag, V. Kostina, and B. Hassibi.  Feedback capacity of mimo gaussian channels.  arXiv preprint arXiv:2106.01994, 2021.
  * [43] DR Rigney, AL Goldberger, WC Ocasio, Y Ichimaru, GB Moody, and RG Mark.  Multi-channel physiological data: description and analysis.  In AS Weigend and NA Gershenfeld, editors, Time Series Prediction: Forecasting the Future and Understanding the Past, pages 105–129. Addison-Wesley, Reading, MA, 1993.
  * [44] A Goldberger, L Amaral, L Glass, J Hausdorff, PC Ivanov, R Mark, and HE Stanley.  Physiobank, physiotoolkit, and physionet: Components of a new research resource for complex physiologic signals.  Circulation [Online], 101(23):e215–e220, 2000.
  * [45] Dor Tsur, Bashar Huleihel, and Haim Permuter.  Rate distortion via constrained estimated mutual information minimization.  In 2023 IEEE International Symposium on Information Theory (ISIT), pages 695–700. IEEE, 2023.
  * [46] L. Breiman.  The individual ergodic theorem of information theory.  The Annals of Mathematical Statistics, 28(3):809–811, 1957.
  * [47] Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980, 2014.
  * [48] Djork-Arné Clevert, Thomas Unterthiner, and Sepp Hochreiter.  Fast and accurate deep network learning by exponential linear units (elus).  arXiv preprint arXiv:1511.07289, 2015.