# A Low-Complexity ADMM-based Massive MIMO Detectors via Deep Neural Networks

###### Abstract

An alternate direction method of multipliers (ADMM)-based detectors can
achieve good performance in both small and large-scale multiple-input
multiple-output (MIMO) systems. However, due to the difficulty of choosing the
optimal penalty parameters, their performance is limited. This paper presents
a deep neural network (DNN)-based massive MIMO detection method which can
overcome the above limitation. It exploits the unfolding technique and learns
to estimate the penalty parameters. Additionally, a computationally cheaper
detector is also proposed. The proposed methods can handle the higher-order
modulation signals. Numerical results are presented to demonstrate the
performances of the proposed methods compared with the existing works.

Index Terms—  Deep learning, sharing-ADMM, MIMO detection

##  1 Introduction

Several deep neural networks (DNN)-based detectors are proposed in the
literature to improve the efficiency of existing traditional multi-input
multi-output (MIMO) detectors [1] [2]. These detectors are designed by using
the “deep unfolding” technique [3]. This technique considers the iterations of
existing algorithms to be a network layer; and then learns from the data to
estimate optimal parameters for the underlying algorithm. For example, the
work in [1] proposed a network called “DetNet” that unfolds the iterations of
the projected gradient descent. Similar idea is used in the works of [4]–[6].

Further, detectors based on the alternating direction multiplier (ADMM)
algorithm are also proposed, both in the traditional and deep learning
approaches. In such cases, since ADMM is used to solve the optimal maximum
likelihood (ML) problem with some constraints, the resulting detectors can
achieve a sub-optimal performance. For instance, the authors of [7] proposed
an infinity norm constrained traditional detector known as “ADMIN” based on
the ADMM. Further, the work in [8] proposed a network called “ADMM-Net” by
unfolding the iterations of the ADMM algorithm. The ADMM-Net has a similar
problem structure as the ADMIN, but it is limited to the BPSK and QPSK
modulations. Furthermore, recently, another ADMM-based detector known as “PS-
ADMM” is proposed in [9]. In that work, the authors have proposed a penalized
ML detection formulation be solved via the sharing-ADMM method that can handle
high-order QAM signals efficiently. They have also shown that the theoretical
performance of the PS-ADMM is guaranteed. However, since there is no
analytical method for obtaining the optimal penalty parameters [11] [12], it
is challenging to improve the performance in the approaches mentioned above.

Inspired by the work of [8] and [9], and using the unfolding method, we
propose an efficient DNN-based ADMM detector. Our work extends the existing
algorithm of [9] in two different forms. First, we unfold the iterations of
the PS-ADMM and design a DNN architecture that can efficiently approximate the
optimal penalty parameters. Second, we customize the standard multilayer
perceptron (MLP) [13] into the ADMM and utilize different feature selection
strategies to obtain improved performance with lower computational complexity.

##  2 System Model and Problem Formulation

We consider a MIMO system in which a base station (BS) is equipped with
$M_{c}$ antenna arrays to serve a user terminal equipped with $K_{c}$ number
of antennas, where $M_{c}>K_{c}$. We assume that a transmit symbol vector
$\tilde{{\bf s}}$ is drawn from a set $\mathcal{S}_{c}$ in a rectangular
$4^{q}$-QAM complex alphabets, where $q$ is a positive integer. Then the
complex received signal vector $\tilde{{\bf y}}$ can be written as

$\displaystyle\tilde{\mathbf{y}}=\tilde{\mathbf{H}}\tilde{\mathbf{s}}+\tilde{\bm{\upsilon}},$
(1)

where ${\tilde{\bm{\upsilon}}}\in\mathbf{\mathbb{C}}^{M_{c}}$ is a zero mean,
unit variance complex additive white Gaussian noise (AWGN), and
$\tilde{\mathbf{H}}$ is a flat fading i.i.d Rayleigh channel. We transform (1)
to an equivalent real model

$\displaystyle\mathbf{y}=\mathbf{H}\mathbf{s}+\bm{\upsilon},$ (2)

where

$\displaystyle\mathbf{s}\\!=\\!\begin{bmatrix}{\rm Re}(\tilde{\mathbf{s}})\\\
{\rm Im}(\tilde{\mathbf{s}})\end{bmatrix}\
\mathbf{y}\\!=\\!\begin{bmatrix}{\rm Re}(\tilde{\mathbf{y}})\\\ {\rm
Im}(\tilde{\mathbf{y}})\end{bmatrix}\ \mathbf{H}\\!=\\!\begin{bmatrix}{\rm
Re}(\tilde{\mathbf{H}})\ &{\rm Im}(\tilde{\mathbf{H}})\\\ -{\rm
Im}(\tilde{\mathbf{H}})\ &{\rm Re}(\tilde{\mathbf{H}})\end{bmatrix},$

$\bm{\upsilon}\in\mathbb{R}^{M}$ is the AWGN. The real domain MIMO size and
symbol sets can be written as $(M,K)=(2M_{c},2K_{c})$, and
$\mathcal{S}=\\{{\rm Re}({\mathbf{\mathcal{S}}_{c}})\\}\bigcup\\{{\rm
Im}({\mathbf{\mathcal{S}}_{c}})\\}$ respectively.

Generally for any given $4^{q}$-QAM signal in the set $\mathcal{S}$ the
transmit symbol can be written as

$\displaystyle\mathbf{s}=\sum_{i=1}^{q}2^{i-1}\mathbf{z}_{i},$ (3)

where ${\bf z}_{i}\in\\{-1,1\\}^{K}$ for all $i=1,2,\cdots,q$. Now let us
introduce an auxiliary variable ${\bf x}\in\mathbb{R}^{K}$, and consider the
binary vectors $\\{{\bf z}_{i}\\}$ as penalty terms. Then by relaxing the
constraint of the classical ML problem, the MIMO detection problem can be
written as

$\displaystyle\begin{split}&\underset{{\bf x},{\bf
z}_{i}}{\min}\frac{1}{2}\|\mathbf{y}-\mathbf{H}\mathbf{x}\|^{2}_{2}-\frac{1}{2}\sum_{i=1}^{q}\alpha_{i}\|\mathbf{z}_{i}\|^{2}_{2}\\\
&\text{s.t.}\quad\mathbf{x}-\sum_{i=1}^{q}2^{i-1}\mathbf{z}_{i}=0,\quad{\bf
z}_{i}\in\mathcal{Q}^{K}\end{split}$ (4)

where $\mathcal{Q}=[-1,1]$ is the box constraint that relaxes the binary set
$\\{-1,1\\}^{K}$, and $\alpha_{i}\geq 0$ are the penalty parameters. The
augmented Lagrangian corresponding to (4) can be written as

$\displaystyle\begin{split}\mathcal{L}_{\rho}(\\{{\bf
z}_{i}|i=1,2,\cdots,q\\},\mathbf{x},\bm{\lambda})=\frac{1}{2}\|\mathbf{y}-\mathbf{H}\mathbf{x}\|^{2}_{2}-\frac{1}{2}\sum_{i=1}^{q}\alpha_{i}\|\mathbf{z}_{i}\|^{2}_{2}\\\
+\bm{\lambda}^{T}\left(\mathbf{x}-\sum_{i=1}^{q}2^{i-1}\mathbf{z}_{i}\right)+\frac{\rho}{2}\|\mathbf{x}-\sum_{i=1}^{q}2^{i-1}\mathbf{z}_{i}\|^{2}_{2},\end{split}$
(5)

where $\bm{\lambda}\in\mathbb{R}^{K},\,\rho,\,T$ are the dual variable,
penalty parameter, and transpose operation respectively. By defining the
scaled dual variable as ${\bf u}=\frac{\bm{\lambda}}{\rho}$, and solving (5)
through ADMM in its equivalent scaled-form, we can obtain the following
sequence of solutions [9]

$\displaystyle{\bf z}_{i}^{t+1}:=$
$\displaystyle{\Pi}_{[-1,1]}\bigg{(}\frac{2^{i-1}\rho}{4^{i-1}\rho-\alpha_{i}}\bigg{(}{\bf
x}^{t}-\sum_{j^{\prime}=1,i\neq 1}^{i-1}2^{j^{\prime}-1}{\bf
z}_{j^{\prime}}^{t+1}$ $\displaystyle-\sum_{j=i+1,q>1}^{q}2^{j-1}{\bf
z}_{j}^{t}+{\bf u}^{t}\bigg{)}\bigg{)},i=1,2,\cdots,q$ (6a) $\displaystyle{\bf
x}^{t+1}:=$ $\displaystyle\bigg{(}{\bf H}^{T}{\bf H}+\rho{\bf
I}\bigg{)}^{-1}\bigg{(}{\bf H}^{T}{\bf y}$
$\displaystyle+\rho\bigg{(}\sum_{i=1}^{q}2^{i-1}\mathbf{z}^{t+1}_{i}-{\bf
u}^{t}\bigg{)}\bigg{)},$ (6b) $\displaystyle{\bf u}^{t+1}:=$
$\displaystyle{\bf u}^{t}+{\bf
x}^{t+1}-\sum_{i=1}^{q}2^{i-1}\mathbf{z}^{t+1}_{i},$ (6c)

where $t=1,2,\cdots,\tau$ is the number of iteration, ${\Pi}(\cdot)$ is a
projection.

From the above procedure, we can see that the detection framework in (6) can
achieve a better solution and theoretically-guaranteed convergence [12].
However, one of the issues with this strategy is the choice of penalty
parameters, as there is no theoretical procedure for choosing the optimal
values [11]. Besides that, since (6b) involves the matrix inversion operation,
it incurs higher computation cost, to the overall complexity of the detector.

To overcome these challenges, we propose DNN-based solutions by transforming
(6) into the deep network architecture. Unlike the work in [8], which is
limited to BPSK and QPSK signals, our method can detect QAM signals in higher
orders efficiently. In addition, we will present a low-complexity learning
strategy for ADMM-based detection.

##  3 The proposed DNN-based ADMM Detectors

In this section, we present two distinct detectors based on the ADMM solutions
(6).

We begin by describing the deep unfolding procedure that can transform (6)
into sequentially connected layers of DNN as follows. Suppose that ${\bf
w}^{\ell}_{1}\in\mathbb{R}^{K}$ denotes a linear transformation vector, where
$\ell=1,2,\cdots,L$ is the number of layers. Then we define

$\displaystyle{\bf
w}_{i}^{\ell+1}=\frac{2^{i-1}\rho}{4^{i-1}\rho-\alpha_{i}}\bigg{(}{\bf
x}^{\ell}-\sum_{{j=1,i\neq 1}}^{i-1}2^{j-1}{\bf z}_{j}^{\ell+1}$
$\displaystyle-\sum_{j=i+1}^{q}2^{j-1}{\bf z}_{j}^{\ell}+{\bf
u}^{\ell}\bigg{)},i=1,2,\cdots,q$ (7)

as linear transformations at the $\ell$-th layer, with a set of trainable
parameters $\bm{\theta}=\\{\\{\alpha_{i}\\}_{i=1}^{q},\rho\\}$, where ${\bf
x}^{\ell},{\bf u}^{\ell},\\{{\bf z}_{j}^{\ell}\\}^{q}_{j=i+1}$ are the inputs.
Next, we untie the projection $\Pi$ and define a pointwise activation function
as

$\displaystyle{\rm sgnlin}({\bf w})=\min(\max(-1,{\bf w}),1),$ (8)

where ${\bf w}\in\mathbb{R}^{K}$ is any vector.

On the basis of the above transformations, we subsequently present the network
architectures of the proposed detectors.

###  3.1 ADMM-PSNet

The goal of this network is to learn the set of optimal penalty parameters
$\bm{\theta}$, from the data generated by the model (6). Its network
architecture is constructed as follows.

$\displaystyle\text{Initialize }{\bf x}^{0}=$ $\displaystyle{\bf 0},{\bf
u}^{0}={\bf 0},\\{{\bf z}_{j}^{0}\\}^{q}_{j=i+1}={\bf 0},i=1,2,\cdots,q;$
$\displaystyle\text{For }\ell=$ $\displaystyle 1,2,\cdots,L:$
$\displaystyle{\bf z}_{i}^{\ell+1}=$ $\displaystyle{\rm sgnlin}({\bf
w}_{i}^{\ell+1}),$ (9a) $\displaystyle{\bf
x}^{\ell+1}\big{(}\bm{\theta}\big{)}=$ $\displaystyle\bigg{(}{\bf H}^{T}{\bf
H}+\rho{\bf I}\bigg{)}^{-1}\bigg{(}{\bf H}^{T}{\bf y}$
$\displaystyle+\rho\big{(}\sum_{i=1}^{q}2^{i-1}\mathbf{z}^{\ell+1}_{i}-{\bf
u}^{\ell}\big{)}\bigg{)},$ (9b) $\displaystyle{\bf u}^{\ell+1}=$
$\displaystyle{\bf u}^{\ell}+{\bf
x}^{\ell+1}\big{(}\bm{\theta}\big{)}-\sum_{i=1}^{q}2^{i-1}\mathbf{z}^{\ell+1}_{i}.$
(9c)

This network is trained in an end-to-end fashion. In the forward pass, the
parameter $\bm{\theta}$ is randomly initialized and kept constant in each
layer. Then information propagates through each layer according to the update
rules in (9). At the output ( $L$-th layer), it computes ${\bf x}^{L}$ and
back propagates to update the parameter $\bm{\theta}$ by minimizing the loss
function

$\displaystyle\mathcal{L}\big{(}\bm{\theta}\big{)}=\frac{1}{m}\sum_{i^{\prime}=1}^{m}\big{\|}{\hat{{\bf
x}}}^{L(i)}\big{(}\bm{\theta}\big{)}-{\bf s}^{(i)}\big{\|}_{2}^{2},$ (10)

through the stochastic gradient descent (SGD), where $m$ denotes the number of
training examples. This process will proceed until the acceptable error limit
has been met. After the training is completed, ADMM-PSNet uses the trained
parameters for online detection.

###  3.2 ADMM-HNet

This network introduces a hidden layer MLP in each iteration of (6). We begin
the discussion by noting the computationally costly term $\left({\bf
H}^{T}{\bf H}+\rho{\bf I}\right)^{-1}$ in the ${\bf x}$-update (6b). Our goal
is to obtain an accurate estimate of the ${\bf x}$-update without analytically
computing this term by utilizing the approximation power of the MLP networks
[13]. To achieve this, first, we assume that the ${\bf x}$-update (6b) is
unknown and describe it as a mapping function $\psi$:

$\displaystyle{\bf x}^{\ell+1}=\psi^{\ell}\big{(}{\bf y},{\bf H},\\{{\bf
z}_{i}^{\ell+1}\\}_{i=1}^{q},\rho,{\bf
u}^{\ell}\big{)},\quad\ell=1,2,\cdots,L.$ (11)

Then, to efficiently approximate the function $\psi^{\ell}$ we introduce an
MLP network $g_{n}^{\ell}$ that has one hidden layer and $n$ number of units
into the $\ell$-th layer. The key task here is to obtain a proper input
(features) of the $g_{n}^{\ell}$ that leads to an effective learning
environment without requiring costly computations. One way to do this would be
to combine all the inputs at each iteration; i.e. stacking the set ( ${\bf
y},{\bf H},{\rho},\\{{\bf z}_{i}^{\ell+1}\\}_{i=1}^{q},{\bf u}^{\ell}$ ) as a
feature vector at each ${\bf x}$-update. However, since the received signal
${\bf y}$ and the channel ${\bf H}$ are constants through ADMM iterations, no
new information will be provided by these parameters in each layer.
Consequently, the learning would be stuck, at some local minima, since ${\bf
y}$ and ${\bf H}$ are dominating the input features. As such, in this work, we
introduce the following procedure to obtain efficient training features.

Instead of directly using ${\bf y}$ and ${\bf H}$ in the input, let us
introduce a residual variable ${\bf r}$, defined as

$\displaystyle{\bf r}^{\ell+1}={\bf y}-{\bf H}{{\bf x}}^{\ell},$ (12)

in each iteration. Then we define a new feature vector ${\bf
a}_{1}^{\ell}\in\mathbb{R}^{K}$ as

$\displaystyle{\bf a}_{1}^{\ell}=\frac{1}{M}\big{(}{\bf H}^{T}{\bf
r}^{\ell+1}\big{)}+{{\bf x}}^{\ell}.$ (13)

Notice that the first term in (13) works as successive noise cancellation, and
the second term is the signal of interest from the previous layer. As the
number of layers increases, since ${\bf a}_{1}^{\ell}$ becomes more associated
with the estimated signal, it can be regarded as a subset of the potential
variables to the input. In addition, we choose

$\displaystyle{\bf
a}^{\ell}_{2}=\rho\big{(}\sum_{i=1}^{q}2^{i-1}\mathbf{z}^{\ell+1}_{i}-{\bf
u}^{\ell}\big{)},$ (14)

as another potential input variable since it naturally updates through the
ADMM iterations. Hence, the potential input feature set for the $\ell$-th
layer can be obtained by vertically stacking the above features as

$\displaystyle{\bf a}^{\ell}=[{\bf a}_{1}^{\ell T},{\bf a}_{2}^{\ell T}]^{T}.$
(15)

Notice that the dimension of each subset feature vector in (15) is $K\times
1$. Hence, the size of the input feature vector at each layer can be defined
as $d_{0}=2K$.

Given the above feature vector, the architecture of the $g_{n}^{\ell}$ can be
described as follows. Suppose that ${\bf W}_{1}^{\ell+1}\in\mathbb{R}^{n\times
d_{0}},{\bf W}_{2}^{\ell+1}\in\mathbb{R}^{K\times n}$ and ${\bf
b}^{\ell}_{1}\in\mathbb{R}^{n},{\bf b}^{\ell}_{2}\in\mathbb{R}^{K}$ denote the
weight matrices and the corresponding bias terms respectively. Then the
forward pass of the $g_{n}^{\ell}$ at the $\ell$-th layer can be explicitly
defined as

$\displaystyle\begin{split}{\bf t}^{\ell+1}=\Gamma\big{(}{\bf
W}_{1}^{\ell+1}{\bf a}^{\ell}+{\bf b}_{1}^{\ell+1}\big{)},\,\hat{{\bf
x}}^{\ell+1}={\bf W}_{2}^{\ell+1}{\bf t}^{\ell+1}+{\bf
b}_{2}^{\ell+1},\end{split}$ (16)

where $\Gamma(.)$ is a pointwise ReLU activation: $\max(.,0)$, and ${\bf
t}^{\ell+1}\in\mathbb{R}^{n}$ is an intermediate variable. Each network
$g_{n}^{\ell}$ will be trained by minimizing the loss

$\displaystyle\mathcal{L}\big{(}\bm{\Theta}^{\ell+1}\big{)}=\frac{1}{m}\sum_{i^{\prime}=1}^{m}\big{\|}g_{n}^{\ell}({\bf
a}^{\ell(i)};\bm{\Theta}^{\ell+1})-{\bf s}^{(i)}\big{\|}_{2}^{2},$ (17)

via SGD, where $\bm{\Theta}^{\ell+1}=\\{{\bf W}_{1}^{\ell+1},{\bf
W}_{2}^{\ell+1},{\bf b}_{1}^{\ell+1},{\bf b}_{2}^{\ell+1}\\}$. Algorithm 1
summarizes the training procedure of the ADMM-HNet.

Input: ${\bf y},{\bf H},L,q,m$, $\\{\\{\alpha_{i}\\}_{i=1}^{q},\rho\\}$;

1 Initialization: $\hat{{\bf x}}^{0}={\bf H}^{T}{\bf y}/M,{\bf u}^{0}={\bf
0},\\{{\bf z}_{j}^{0}\\}^{q}_{j=i+1}={\bf 0}$, $\bm{\Theta}^{\ell}=\\{{\bf
W}_{1}^{\ell},{\bf W}_{2}^{\ell},{\bf b}_{1}^{\ell},{\bf
b}_{2}^{\ell}\\}^{L}_{\ell=1}$;

2 for  _$\ell=1,2,\cdots,L$ _ do

3 Update $\\{{\bf z}_{i}^{\ell+1}\\}$ via (9a);

4 Generate the input features through (13)-(15);

5 Construct the network $g_{n}^{\ell}$ according to (16);

6 for  _$i^{\prime}=1,2,\cdots,m$_ do

7 Construct the dataset ${\bf D}=\\{{\bf a}^{\ell(i^{\prime})},{\bf
s}^{(i^{\prime})}\\}$;

8 Train the $g_{n}^{\ell}$ for some batchsize;

9 end for

10 Return $\hat{{\bf x}}^{\ell+1}=g_{n}^{\ell}\big{(}{\bf
a}^{\ell},{\bm{\Theta}^{\ell+1}}\big{)}$;

11 Update the dual variable: ${\bf u}^{\ell+1}={\bf u}^{\ell}+{\hat{{\bf
x}}}^{\ell+1}-{\bf z}^{\ell+1}$;

12

13 end for

Output: Return $\hat{{\bf x}}^{L}$.

Algorithm 1 Pseudo code fo the ADMM-HNet.

We conclude this discussion with the following technical remarks on penalty
parameters and computational complexity. (1) In the ADMM-HNet above, since the
penalty parameters do not appear in the network layers (see equation (16)), it
is not convenient to learn them explicitly in the hidden layers. As such, we
have pre-trained them in our implementation using the ADMM-PSNet model. (2)
The computational complexity of the ADMM-PSNet shall be the same as that of
the PS-ADMM and ADMIN detectors. However, it is significantly reduced for the
ADMM-HNet. From the procedures in section 3.2 and Algorithm 1, we can see that
the complexity of this detector grows approximately in the order of
$O(MK+L(MK+3Kn))$ where $n$ is the number of units.

##  4 Numerical Results

In this section, different numerical results are presented to demonstrate the
effectiveness of the proposed methods. All detectors are implemented in python
3.6 and Keras API on top of TensorFlow. The ADMM-PSNet and ADMM-HNet are
trained on $10,000$ and $90,000$ for samples, respectively. The ADMM-PSNet is
trained in an end-to-end manner with $L\\!=\\!30,{\rm epochs\\!=\\!10,000}$,
while the ADMM-HNet, is trained in a layer-wise with $\rm
batchsize\\!=\\!1024,L\\!=\\!30,\rm epochs\\!=\\!2000$.

Fig. 1 shows the generalization ability of our newly proposed ADMM-HNet in
terms of its symbol-error-rate (SER) performance.

Fig. 1: The SER performances of the ADMM-HNet for a 16-QAM, $M=64,K=16$ for
different number of layers.

We run a Monte Carlo simulation on it with a newly generated data to examine
its performance at different layers. It is clear from the figure that as the
number of layers increases the performance also increased. Particularly when
$L$ is in the range $1$ to $20$, the performance increased significantly. On
the other hand, when $L$ increases beyond $20$, we can observe that it
achieved nearly equal performances. This implies that, unlike some existing
detectors, the ADMM-HNet does not require higher number of layers to achieve a
competitive performance; i.e. with a properly selected input feature, and
learned penalty parameters it can converge at early stage (such as $L=20$).
Next, to examine the effectiveness of the proposed methods, we compare their
performances with the existing detectors.

Fig. 2 illustrates the SER versus SNR performance comparison between the
proposed and existing detectors for the 16-QAM, $M\\!=\\!64,K\\!=\\!16$ MIMO
system. As we can see the ADMM-PSNet performance is nearly similar to that of
the PS-ADMM detector, whereas the ADMM-HNet outperformed the remaining
detectors. Notice that the number of layers for this detector is 20, while the
others run up to 30.

Fig. 2: The SER performance comparison for a 16-QAM, $M=64,K=16$.

Fig. 3 depicts the SER versus SNR performances of the detectors for the
64-QAM, $M\\!=\\!64,K\\!=\\!16$ MIMO case. In this figure, we can see that the
ADMM-PSNet achieved slightly better performance than the PS-ADMM under the
higher SNR regime. On the other hand, ADMM-HNet performed almost similar to
ADMM-PSNet and PS-ADMM under the lower SNR regime and improved performance
under the higher SNR regime.

Fig. 3: The SER performance comparison for a 64-QAM, $M\\!=\\!64,K\\!=\\!16$.

These imply that with properly selected input features, the proposed ADMM-HNet
can achieve a very competitive performance with a moderate number of layers
such as 20. Moreover, as discussed in section 3.2, since its computational
complexity is cheaper compared to the others, this detector can efficiently
serve the massive MIMO system.

We conclude our discussion by comparing the average runtime between ADMM-based
detectors in the table 1. These results are measured on the CPU: Intel(R)
Core(TM) i3-7100U CPU $3.90GHz\times 2$.

Table 1: The average run time comparison measured in seconds; “*” stands for “ADMM”. Mod, $M\times K$ | Detectors
---|---
| *-HNet | PS-* | DetNet | ADMIN
16-QAM, $128\times 64$ | 0.0029 | 0.0093 | 0.015 | 0.091
64-QAM, $64\times 16$ | 0.0018 | 0.0071 | 0.0083 | 0.076

In this table, we can see that the ADMM-HNet reported faster time in both MIMO
configurations. This is due to the fact that it is constructed from a cascade
of a one-hidden-layer network that only computes a matrix-vector
multiplication. The others, however, compute the multiplication and inversion
of the matrix in each layer (iteration), and hence, they are relatively
slower.

##  5 Conclusion

In this paper, DNN massive MIMO detection methods that can work with the
higher-order QAM signals are proposed based on the ADMM algorithm. The ADMM-
PSNet is trained to estimate the penalty parameters that result in improved
performance. On the other hand, the ADMM-HNet is designed to reduce
computational complexity by eliminating the matrix inversion and
multiplication operations. The numerical results reveal that the proposed
detectors can significantly improve the performance of the existing methods.
In Particular, the ADMM-HNet has shown improved performance with less
computational complexity.

## References

  * [1] N. Samuel, T. Diskin, and A. Wiesel, “Learning to detect,” _IEEE Trans. Signal Process._ , vol. 67, no. 10, pp. 2554–2564, Feb. 2019.
  * [2] H. He, C. Wen, S. Jin, and G. Y. Li, “Model-driven deep learning for MIMO detection,” _IEEE Trans. Signal Process._ , vol. 68, pp. 1702–1715, Feb. 2020.
  * [3] J. Hershey, J. Le Roux, and F. Weninger, “Deep unfolding: Model-based inspiration of novel deep architectures,” _arXiv preprint arXiv:1409.2574,_ 2014\.
  * [4] M. Khani, M. Alizadeh, J. Hoydis, and P. Fleming, “Adaptive neural signal detection for massive MIMO,” _IEEE Trans. Wireless Commun._ , vol. 19, no. 8, pp. 5635–5648, May 2020.
  * [5] X. Tan, W. Xu, K. Sun, Y. Xu, Y. Be’ery, X. You, and C. Zhang, “Improving massive MIMO message passing detectors With deep neural network,” _IEEE Trans. Veh. Technol._ , vol. 69, no. 2, pp. 1267–1280, Dec. 2020.
  * [6] Y. Wei, M.-M. Zhao, M. Hong, M.-J. Zhao, and M. Lei, “Learned conjugate gradient descent network for massive MIMO detection,” _IEEE Trans. Signal Process._ , vol. 68, pp. 6336–6349, Nov. 2020.
  * [7] S. Shahabuddin, M. Juntti, and C. Studer, “ADMM-based infinity norm detection for large MU-MIMO: Algorithm and VLSI architecture,” in _Proc. IEEE Int. Symp. Circuits Syst._ , 2017, pp. 1–4.
  * [8] M. Un, M. Shao, W. Ma, and P. C. Ching, “Deep MIMO detection using ADMM unfolding,” in _Proc. IEEE Data Sci. Workshop_ , Jun. 2019, pp. 333–337.
  * [9] Q. Zhang and Y. Wang, “Designing massive MIMO detector via PS-ADMM approach,” _arXiv preprint arXiv:2010.06798_ , Oct. 2020.
  * [10] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein, “Distributed optimization and statistical learning via the alternating direction method of multipliers,” _Foundations Trends Mach. Learn._ , vol. 3, no. 1, pp. 1–122, 2011.
  * [11] B. Wohlberg, “ADMM penalty parameter selection by residual balancing,” _arXiv preprint arXiv:1704.06209_ , Apr. 2017.
  * [12] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein, “Distributed optimization and statistical learning via the alternating direction method of multipliers,” _Foundations Trends Mach. Learn._ , vol. 3, no. 1, pp. 1–122, 2011.
  * [13] K. Hornik, M. Stinchcombe, H. White, “Multilayer feedforward networks are universal approximators,” _Neural networks_ , vol. 2, no. 5, pp. 359–366, Jul. 1989.


# Detecting Causality in the Frequency Domain with Cross-Mapping Coherence

Zsigmond Benkő  HUN-REN Wigner Research Centre for Physics, Department of
Computational Sciences, Budapest, Hungary  Bálint Varga  HUN-REN Wigner
Research Centre for Physics, Department of Computational Sciences, Budapest,
Hungary  Semmelweis University, Semmelweis Doctoral School, János Szentágothai
Neurosciences Division, Budapest, Hungary  Marcell Stippinger  HUN-REN Wigner
Research Centre for Physics, Department of Computational Sciences, Budapest,
Hungary  Zoltán Somogyvári  HUN-REN Wigner Research Centre for Physics,
Department of Computational Sciences, Budapest, Hungary

###### Abstract

Understanding causal relationships within a system is crucial for uncovering
its underlying mechanisms. Causal discovery methods, which facilitate the
construction of such models from time-series data, hold the potential to
significantly advance scientific and engineering fields.

This study introduces the Cross-Mapping Coherence (CMC) method, designed to
reveal causal connections in the frequency domain between time series. CMC
builds upon nonlinear state-space reconstruction and extends the Convergent
Cross-Mapping algorithm to the frequency domain by utilizing coherence metrics
for evaluation. We tested the Cross-Mapping Coherence method using simulations
of logistic maps, Lorenz systems, Kuramoto oscillators, and the Wilson-Cowan
model of the visual cortex. CMC accurately identified the direction of causal
connections in all simulated scenarios. When applied to the Wilson-Cowan
model, CMC yielded consistent results similar to spectral Granger causality.

Furthermore, CMC exhibits high sensitivity in detecting weak connections,
demonstrates sample efficiency, and maintains robustness in the presence of
noise.

In conclusion, the capability to determine directed causal influences across
different frequency bands allows CMC to provide valuable insights into the
dynamics of complex, nonlinear systems.

##  1 Introduction

###  1.1 Causal discovery from time series

A causal understanding of the world is crucial to building robust and
explainable models for prediction and control. Where interventions are not
possible, observational causal discovery methods help to determine the
underlying causal structures from static or time-series data. Despite its
apparent simplicity, this is a non-trivial task because there can be multiple
causal structures between the variables, including unidirectional causation,
bidirectional causality, hidden confounders, or even independence [1, 2, 3, 4,
5, 6]. Further complicating the matter is the fact that the underlying systems
can exhibit either stochastic or deterministic characteristics. The task of a
causal discovery algorithm is to identify the actual structure from
observational data, if possible (Fig. 1).

For time series, there are many approaches to discovering causal structures,
such as the methods based on the predictive causality principle or the theory
of nonlinear dynamical systems. Statistical and information-theoretic methods
work by following the principle of predictive causality [7, 8, 9, 10, 11, 12,
13, 3, 14, 15, 16]. According to this, given two time series, one is a cause
of the other if its values help in predicting the future values of the other.
Granger Causality (GC) [8] and Transfer Entropy (TE) [9, 10] are two members
of this family with frequency-domain extensions such as Spectral Granger
Causality [11, 12] and Transfer Entropy Spectrum [13]. These methods work well
for linear systems and for unidirectional couplings, but despite recent
extensions [3, 14, 15, 16], nonlinearity and hidden confounders may pose a
significant obstacle to the detection of causal relationships.

The nonlinear dynamics approach leverages the properties of dynamical systems
to reveal causal connections [17, 18, 4, 19, 20, 21]. Mathematical theorems
[22, 23, 24] enable the reconstruction of system states from time series,
allowing for the comparison of trajectories to discern causal relationships
among variables. Prominent methodologies within this domain include Convergent
Cross-Mapping (CCM) [17, 18, 4], the technique of recurrence plots [19], and
causal inference based on manifold dimensions [20, 2] or complexity measures
[6], all of which are predominantly applied within the time domain. These
methodologies are efficient in identifying weak directed and circular
couplings, as well as revealing latent confounding factors within nonlinear
systems, although their applicability is restricted to time-domain
information. The extension to the spectral domain introduces a novel dimension
for causal discovery. In this study, we pursue this objective by extending CCM
to the frequency domain.

Figure 1: Causal discovery from time series. Causal discovery algorithms try
to identify causal structures between time series, whether it is independence,
unidirectional direct or circular causation, or hidden confounding.

###  1.2 Dynamical systems and Convergent Cross Mapping

CCM leverages the mathematical framework of nonlinear dynamical systems to
identify causal relationships from time series data. In this section, we
review the fundamental concepts of non-linear dynamical systems and the
methodology of causality detection using the CCM technique.

A dynamical system constitutes a self-mapping of the state space. The
subsequent state is determined by the current state in conjunction with the
mapping. This mapping may exhibit deterministic or stochastic characteristics,
the temporal parameter may be continuous or discrete, and the system’s
dynamics can be either linear or nonlinear. For discrete time dynamical
systems, the next state is generated by the update rule:

$s_{t+1}=f(s_{t})$ (1)

where $f:\mathcal{S}\rightarrow\mathcal{S}$ is the mapping. Illustrative
instances of nonlinear discrete-time dynamical systems encompass the logistic
map, the tent map [25], and the Labos neuron model [26], among others.

In the continuous-time case we can describe the dynamics by differential
equations:

$\dot{s}=\varphi(s)$ (2)

where the time derivative of the state $s$ defines a flow in the state space.
The trajectories of the system can be obtained from an initial condition
$s_{0}$ by solving the differential equations
$s_{T}=\int_{0}^{T}\varphi(s_{t})\mathrm{d}t$. Prominent examples of nonlinear
continuous-time dynamical models include the Lorenz system [27], the Rössler
system [28], the Lotka–Volterra equations [29, 30], and the Wilson–Cowan model
of cortical dynamics [31].

In practical scenarios, direct observation of a system’s state is often
unfeasible, with only time-series data being accessible. This one-dimensional
data shows only a, possibly non-linear, projection of the state, and the
reconstruction of the underlying high-dimensional state from these
observations is thus a pivotal task. Fortunately, Takens’ theorem offers a
robust mathematical framework for state reconstruction from time series data
[22]. Time Delay Embedding constitutes a methodology for non-linear state
space reconstruction, wherein one augments a scalar time series
$x(t)\in\mathbb{R}$ to reconstruct the state $s(t)\in\mathbb{R}^{d}$ of the
underlying dynamical system:

$X(t)=[x(t),x(t-\tau),...,x(t-(E-1)\tau)]$ (3)

where $X(t)\in\mathbb{R}^{E}$ denotes the reconstructed state, with
$E\in\mathbb{N}$ and $\tau\in\mathbb{R}$ representing the embedding dimension
and the embedding delay respectively. Provided that $E>2d$, the reconstruction
occurs with a high probability, where $d$ signifies the intrinsic dimension of
the attractor dynamics.

In the context of two interacting systems, where one system exerts influence
over the other, the driver encodes its state into the state of the driven
system in accordance with the skew product theorem [23]. The state of the
driven system encapsulates comprehensive information about the state of the
driver, whereas the converse does not hold, thereby giving the chance to
identify directed causal links.

Convergent Cross-Mapping (CCM) is a causal discovery algorithm grounded in
Takens’ and skew product theorems [17]. CCM deduces directed causal
relationships between systems based on two time series (Box 1). The algorithm
comprises three fundamental stages: state-space reconstruction, prediction,
and evaluation (Fig. 2). The algorithm’s input consists of two time series,
with the objective of investigating the interconnection between them (Fig. 2a,
b, f). In the reconstruction phase, CCM employs time-delay embedding on the
time series to reconstruct the state of the underlying dynamics (Fig. 2a, c,
g). Consequently, the manifold of each system is reconstructed. The embedding
dimension ($E$) and embedding delay ($\tau$) are critical parameters of this
methodology. In the prediction phase, CCM employs weighted k-nearest neighbors
(kNN) regression models to predict time series from the reconstructed
manifolds in both directions (Fig. 2d, h). The inputs to the regression model
consist of data points within one of the reconstructed state-spaces, while the
outputs are the estimates for the corresponding time series. The number of
neighbors ($k$) is an integer parameter, determined as $E+1$. Subsequently,
the evaluation phase scrutinizes the performance of the prediction models to
infer the causal relationships between the underlying systems (Fig. 2e, i).
The coefficient of determination $R^{2}$ (the squared Pearson correlation)
serves as the standard metric for evaluation. The prediction and evaluation
steps are iterated across varying sample sizes to ascertain whether the
coefficient of determination converges with increasing sample size (Fig. 2j).
Should $R^{2}$ converge to a significantly positive value, the predicted
variable may be inferred as a driver of the model input variable.
Additionally, temporal shifting of the time series and subsequent reanalysis
can elucidate the propagation delay of the effect (Fig. 2k).

### Box 1 : Computing the CCM Score

0. Time series $x_{t},y_{t}$ 1. Time Delay Embedding ($X(t),Y(t)$) 2. Neighborhood search $k$NN around each point in X, store neighbor indices ($\tau_{j}^{t}$) and distances ($d_{tj}$) 3. Compute weights: $w_{tj}=\gamma_{tj}/\sum_{j=1}^{k}\gamma_{tj}$ , $\gamma_{tj}=\exp\left(-\frac{d_{tj}}{d_{t1}}\right)$ 4. Mapping with computed weights: $\hat{y}_{t}=\sum_{j=1}^{k}w_{tj}y\left(\tau_{j}^{t}\right)$, where $\tau_{j}^{t}$ is the time index of the $j$th nearest neighbor of $X(t)$ in the reconstructed state-space. 5. Evaluation of predicted values ($R^{2}=\mathrm{Corr}(\hat{y},y)^{2}$) 6. repeat $2\rightarrow 5$ with higher library lengths (sample size), and check whether $R^{2}$ values are converging or not. If $R^{2}$ converges, then $Y\rightarrow X$. Swap the variables and repeat the computation to check in the other direction for a causal effect. These steps fit into the basic structure of state-space reconstruction ($1.$), prediction ($2-4.$), and evaluation ($5-6.$) scheme shown in Fig. 2 and Fig. 3.

Figure 2: Overview of Convergent Cross Mapping. (a) Overview of the CCM
algorithm from a macroscopic perspective and (b - k) demonstrative stages for
the coupled logistic map system $x\rightarrow y$. (b, f) The inputs consist of
individual time series observations for the $x$ (blue) and $y$ (orange)
variables. (c, g) State-space reconstruction employing time-delay embedding.
The three-dimensional plots depict the reconstructed states of the subsystems
$X$ (blue) and $Y$ (orange). (d, h) Predicted time series $\hat{x}$ (blue) and
$\hat{y}$ (orange). (e, i) Evaluation of the the predictions using the
coefficient of determination metrics ($R^{2}$). The models exhibited high
prediction-quality for the $x$ (blue) variable, whereas low prediction-quality
for $y$ (orange) variable. (j) Convergence as the function of library length.
The $R^{2}$ values are converging in the $x\rightarrow y$ direction (blue),
but do not converge in the $y\rightarrow x$ direction (orange), thereby
indicating a unidirectional causal link. (k) Cross-mapping quality as the
function of time shift between the time series (Cross-Mapping Function). The
negative semi-axis denotes the causal side, where the cause precedes the
effect. The positive semi-axis represents the anti-causal side, where the
predictability of the effect from the cause can be identified. The peak on the
blue curve at zero signifies unidirectional coupling from $x$ to $y$ with
minimal propagation delay. The absence of peaks in the reversed direction
indicates no causal link in the $y$ to $x$ direction (orange curve).

Over the past decade, numerous advancements have been made to the foundational
CCM algorithm to augment its capabilities and enhance its robustness. The
fundamental concepts of mutual cross-mapping [32] are effectively implemented
by the CCM algorithm [17] and have been subsequently refined for asymmetric
inference [33]. Additional extensions have been developed to address temporal
variations in interactions [34, 35], short parallel time series [36], and
potential delays in causal effects [18, 37]. The delay in causal effect is
ascertained by identifying the peak performance of the CCM function with
respect to temporal shifts, wherein the two time series are temporally
displaced relative to one another (Fig. 2k, Fig. 3b). If the optimal shift
resides on the causal side of the time-shift axis, indicating that the cause
precedes the effect, this signifies a delayed causal effect. Conversely, if
the highest peak of the CCM function is located on the anticausal side, it
denotes a Granger causality-like predictability relationship between the
variables [4, 38]. In the presence of predictability, the automatic
determination of directed causal links becomes intricate with the current CCM
methodology, necessitating measures to disentangle causal peaks from the
predictability effects.

We delineate additional advancements within the reconstruction-prediction-
evaluation framework. To enhance state space reconstruction in the CCM
paradigm, methodologies such as Gaussian processes [39, 40], reservoir
computing [41], recurrent neural networks [42], and frequency domain embedding
[43] have been proposed.

To enhance the evaluation and prediction phases, efforts were made to automate
the evaluation procedure and refine the evaluation metric. To eliminate the
subjective element in determining the convergence, automatic fitting of
convergence rates was implemented [44]. Multiple attempts were made to develop
a significance test for the framework, although with limited success [45, 46,
37, 47, 48]. Alternative prediction methodologies and evaluation metrics have
been proposed as modifications to the foundational algorithm to enhance
robustness, including symbolic CCM [49], Convergent Cross Sorting [50],
Gaussian Process priors [51, 52], Neural Shadow Mapping [53], Robust CCM [54],
and a metalevel approach via the Continuity scaling method [55].

These extensions predominantly addressed the bivariate context, and
alternative methodologies were developed to generalize CCM to multivariate
signals, including partial cross-mapping [56] and other techniques [57, 58,
46].

Despite substantial advancements in the Convergent Cross-Mapping (CCM)
algorithm, progress in the domain of spectral causal discovery remains
limited. To date, only two methodologies, Frequency-Domain CCM (FDCCM)[43] and
Cross-Frequency Symbolic CCM (CFSCCM)[59], have made strides in this area.
Frequency-Domain Convergent Cross-Mapping utilizes a frequency-domain state
representation for the reconstruction of state space. However, this approach
is not designed to discern frequency-specific causal interactions. Cross-
Frequency Symbolic CCM, on the other hand, employs the Hilbert transform to
decompose temporal signals and evaluates CCM between the resultant signal
modes, thereby examining phase-amplitude couplings.

In this study, we introduce Cross-Mapping Coherence (CMC), an extension of
Convergent Cross Mapping (CCM) into the frequency domain for the
identification of frequency-specific causal relationships between time series.
Additionally, we propose a novel methodology to disentangle predictability
from dynamic causality. In the following sections, we delineate Cross-Mapping
Coherence and demonstrate its efficacy on coupled systems such as the Logistic
map, the Lorenz system, Kuramoto oscillators, and the Wilson-Cowan model.
Furthermore, we evaluate its robustness with respect to the choice of
embedding parameters and the presence of observational noise. We contextualize
CMC within the framework of spectral causal discovery methods, highlighting
its merits and limitations. Lastly, we provide a comprehensive description of
simulations and analyses in the Methods section.

##  2 Results

###  2.1 Cross-Mapping Coherence

Figure 3: Cross-Mapping Coherence detects causal relationship in the frequency
domain and filters out spurious causality by peak prominence calculation. (a)
A comparative overview of Cross-Mapping Coherence (CMC) and Convergent Cross
Mapping (CCM). Both methodologies analyze two time series and are comprised of
three principal stages: state reconstruction, prediction, and evaluation.
Although the state reconstruction and prediction phases are identical in both
methods, CMC diverges in the evaluation phase. Specifically, CMC utilizes the
coherence metric to compare predicted and actual values, thereby offering
spectral insights into the causal connection. (b) Significant peak
determination with peak prominence calculation on the Cross-Mapping Coherence
Function. The coherence value is computed for each temporal shift, followed by
an examination of the resultant coherence function across individual frequency
bands. An illustrative coherence function is presented at a specific frequency
(blue curve). Peaks may manifest on the anti-causal side of the function
(Granger peak; blue dot) if the state of the effect can be inferred from the
state of the cause. This Granger peak may hinder the causal inference
procedure, so if predictability is present, one has to handle maximum
determination with special care. To reduce the impact of Granger peaks, we
identify peaks on the causal side of the graph, and calculate the prominence
value of each peak to quantify the causal effect.

We expand the Convergent Cross-Mapping (CCM) method to the frequency domain,
allowing us to identify causal links in specific spectral bands. Cross-Mapping
Coherence (CMC) follows the same basic steps as CCM – reconstruction,
prediction, and evaluation – but the evaluation step is altered (Fig. 3a).

To apply CCM, one takes two time series, applying time delay embedding to
reconstruct the state from both (Fig. 2). Then we employ the cross-mapping
step with the k-NN method to predict each other’s states. Finally, we evaluate
the quality of the prediction by computing the coefficient of determination
between the actual and predicted values. To determine the causal effect delay,
we shift the time series relative to each other, compute the CCM values, and
find the maximum of the resulting CCM function.

With CMC, the reconstruction and prediction steps remain the same, but the
evaluation process is different (Fig. 3). Rather than using the coefficient of
determination, we calculate the coherence between the predicted and actual
time series (Eq. 4), providing spectral information about the quality of the
prediction:

$\mathrm{coh}_{x\hat{x}}(f)=\frac{|S_{x\hat{x}}(f)|}{\sqrt{S_{xx}(f)S_{\hat{x}\hat{x}}(f)}}$
(4)

where the coherence $\mathrm{coh}_{x\hat{x}}$ at a specific frequency $f$ is
defined as the magnitude of the averaged cross-spectral density between the
time series $x$ and the predicted values $\hat{x}$, normalized by the square
root of the power spectral densities $S_{xx}$ and $S_{\hat{x}\hat{x}}$ of the
respective signals [60].

To determine the causal effect delay, we shift the time series relative to
each other and recompute the cross-mapping coherence metric resulting in a 2D
shift- and frequency-dependent CMC function.

###  2.2 Effect Delay Determination with Peak Prominence

The time delay dependence of the CCM, as well as CMC, raises the problem of
significant maxima detection: while we can easily assume that higher
correlation and coherence values correspond to stronger connections, simply
taking the maximum of the CCM or the CMC functions does not produce
satisfactory results (Fig. 3b). The problem is, that high values can appear
not only because the cause can be reconstructed from the consequence, but
because the consequence is predictable based on the cause. Consequently, high
peaks can occur at positive time delays, where the consequence precedes the
cause, resulting in apparent contradictions; we call these peaks Granger
peaks. While it is clear that these peaks should be excluded from the
evaluation, they make it more difficult to select the proper peaks, as side
peaks on the causal semi-axis (where cause precedes consequence) can be formed
due to random fluctuations on the shoulder of Granger peaks on the positive
semi-axis. To make the task even more complex, both causal and Granger peaks
can appear simultaneously, especially in the case of bidirectional coupling.

To find the amplitude of the most significant peak on the negative semi-axis,
even in the presence of a possible Granger peak on the positive semi-axis, the
prominence of each local peak on the delay-dependent CMC functions is
calculated for each frequency band. The prominence measure of a peak is used
especially in geography and it expresses how much higher a peak is over its
area of dominance. Specifically, the prominence of a peak corresponds to its
height over the minimum, which should be crossed to get to a higher peak,
while the prominence of the absolute maximum is equal to its elevation.

We measure the strength and delay of the causal relationship, for each
frequency band, by the prominence and delay of the most prominent peak of the
coherence values, where it precedes or co-occurs with the consequence. This
case holds if the time index of the peak is less than the embedding time
window (E in Eq. 3), wherein positive indices are permitted in recognition of
the temporal uncertainty induced by time delay embedding. Thus, if the time
index of the absolute maximum is less than $E$, then the momentary causal
strength is expressed by its value, but if the time index of the absolute
maximum is above $E$, the causal strength is less, corresponding only to the
prominence of a peak below $E$. In this way, we only get high causality if a
well-identifiable independent peak occurs below $E$.

The prominence-based computation facilitates the identification of a genuine
peak on the negative semi-axis, whilst preventing the erroneous recognition of
minor local maxima situated on the shoulder of a more substantial peak on the
positive semi-axis.

In the following subsections, we demonstrate the effectiveness of CMC on time
series of simulated systems.

###  2.3 Logistic maps

Figure 4: Causal discovery on the logistic map systems with diverse coupling
scenarios using Cross-Mapping Coherence. The columns of the figure illustrate
various coupling configurations: unidirectional coupling $\rightarrow$ (a, e,
i, m), circular coupling $x\leftrightarrow y$ (b, f, j, n), hidden common
driver $x\Lsh\Rsh y$ (c, g, k, o), and independence $x\perp y$ (d, h, l, p),
respectively. The direction of the inferred relationship is color-coded: the
$x\rightarrow y$ direction is indicated by blue, whereas the $y\rightarrow x$
direction is marked by orange colors. (a-d) The time-delayed CCM function
reveals causal relations from bivariate time series. A pronounced peak
indicates the presence of causal links in both the unidirectional (a) and
circular (b) scenarios. However, the absence of distinct peaks in the hidden
common driver (c) and independent (d) cases mark the lack of direct causal
connections. (e-l) CMC functions for the coupled logistic map systems show
spectral information for the inferred $x\rightarrow y$ (e-h) and $y\rightarrow
x$ (i-l) directions. The vertical bands on the two-dimensional CMC functions
in the case of unidirectional and circular couplings indicate a roughly
frequency-independent causal link between the time series (e, f, j). The
absence of such a band indicates that there were no direct causal link to be
inferred (i). Here the colormaps are normalized to the $[0,1]$ interval. (m-p)
The peak prominence of CMC values show the proper causal couplings uniformly
along the frequency axis.

To test the CMC method, we performed CCM and CMC analyses on bivariate
logistic map systems with unidirectional, circular, hidden driver, and
independent couplings.

The time-delayed CCM analysis recovered the correct direct causal connections
from the time series, and the CMC investigation revealed near-uniform spectral
contributions (Fig. 4). In the unidirectional and circular coupling cases, we
found high coherence values at every frequency band, while we observed near-
zero coherence values for the hidden driver and independent case.

Figure 5: Dependence of CMC on simulation parameters for logistic maps with
unidirectional coupling. (a) Dependence on time-series length ($D=2$,
$C_{X\rightarrow Y}=0.15$, $\sigma_{\varepsilon}=0$). The CMC function is
convergent, the causal link is hardly detectable for $L=400$, but clearly
visible at $L=700$ simulation steps. (b) Dependence on coupling strength
($L=2000$, $D=2$, $\sigma_{\varepsilon}=0$). The causal link is detected by
CMC even in case of such weak coupling as $C_{X\rightarrow Y}=0.05$. (c)
Dependence on signal-to-noise ratio ($L=10000$, $D=2$, $C_{X\rightarrow
Y}=0.15$ ). CMC can detect the causal link when the signal-to-noise ratio is
at least between $5$ and $10$. Figure 6: Dependence of the CMC function on the
embedding dimension. Time delay embedding introduces temporal uncertainty in
the determination of the effect-propagation delay, thus higher embedding
dimensions result in a plateau in the CMC function proportional to the
embedding window size ($w=(D-1)\tau$).

We investigated the dependence of CMC on simulation parameters for the
unidirectional coupling case (Fig. 5). First, we looked at CMC as a function
of time series length by computing the performance on various sample sizes
from $L=400$ to $L=5000$ with $E=2$ and $\tau=1$ embedding parameters. CMC was
convergent: the directional coupling is not detectable on low sample size
($L=400$), but the pattern of connection becomes clear from a sample size of
$L=1000$ (Fig. 5a). Second, to see how cross-mapping performance depends on
the coupling strength between the systems, we simulated unidirectionally
coupled logistic map systems ($L=2000$, $E=2$, $\tau=1$) by varying the value
of the coupling parameter $C_{X\rightarrow Y}$ in the range $[0,0.2]$. CMC
proved to be very sensitive, as it could detect even weak coupling strengths
($C_{X\rightarrow Y}=0.05$) between the subsystems (Fig. 5b). Third, the noise
robustness of the CMC method was tested by adding observational noise to the
time series ($L=2000$, $E=2$, $\tau=1$), varying the signal to noise ratio
from $1000$ to $2$. CMC robustly detected the coupling down to the signal to
noise value of $10$ (Fig. 5c).

We also examined the dependence of CMC performance on the embedding parameters
for the unidirectionally coupled logistic map system (Fig. 6). Varying the
embedding dimension between $D=2-10$ ($L=10000$, $\tau=1$), we found that as
the dimension increased, a plateau formed on the CMC function along the time-
shift axis proportional to the size of the embedding window, while the
amplitude of the coherence decreased.

###  2.4 Lorenz system

Figure 7: Causal discovery of the Lorenz system in various coupling scenarios.
The columns contain different couplings like unidirectional link $\rightarrow$
(a, d, g, j), circular coupling $\leftrightarrow$ (b, e, h, k), and
independence $\perp$ (c, f, j, l), respectively. (a-c) The time-delayed CCM
score reveals causal relations from bivariate time series. A peak shows the
presence of causal links in the unidirectional (a) and circular (b) cases,
whereas the absence of the peak marks the lack of connection in the
independent case (c). (d-i) CMC function hints the absence of causal links.
CMC scores for the $x_{1}\rightarrow x_{2}$ direction (d-f) are drawn as blue
(2nd row) and the $x_{2}\rightarrow x_{1}$ direction (g-i) as orange (3rd
row). (j-l) The peak prominence CMC values show the true causal couplings
along the frequency axis.

We carried out the CCM and CMC analysis on coupled Lorenz systems to test CMC
on continuous-time dynamics (Fig. 7).

The time-delayed CCM correctly detected the couplings, and CMC returned the
frequency-dependent interactions. In the unidirectional case, a peak can be
identified at $0$ timeshift for the $x\rightarrow y$ direction, and the power
decreases from lower to higher frequencies. In the case of circular coupling,
back-and-forth connections are detected, featuring a decreasing coherence
spectrum with a small peak between $1$ and $2$ Hz. When the subsystems were
independent, the CCM and CMC values were close to zero.

Figure 8: CMC causal discovery on a network of Kuramoto oscillators. (a)
Network structure of coupled Kuramoto oscillators; the $z$ variable drives the
mutually coupled $x$ and $y$ variables. (b) Illustrative time series from the
three Kuramoto oscillators showing sinusoid waveforms of the $x$ (blue), $y$
(orange) and $z$ (green) oscillators with differing central frequencies of
$59$ Hz, $40$ Hz and around $10.5$ Hz, respectively. (c) Coherence between the
signals returns the coupling between the subsystems for the signal pairs. The
coherence has a single peak in for the $z-x$ (blue) and $z-y$ (orange) pairs
around $10$ Hz, whereas for the $x-z$ pair the coherence function has three
peaks around the central frequencies. (d - o) The three columns present the
causality analysis results for the $x-z$ (d, g, j, m), the $y-z$ (e, h, k, n)
and the $x-y$ (f, i, l, o) pairs. (d - f) CCM identifies directional couplings
between the subsystems. Specifically, it shows unidirectional coupling from
$z$ to (d) $x$ and (e) $y$, and bidirectional coupling for the $x-y$ pair (f).
(g - l) CMC scores reveal the frequency-specificity of causal links. Causal
effect of $z$ on $x$ and $y$ shows up as increased CMC around $10$ Hz (g, h),
whereas the lack of feedback is indicated by near-zero CMC values in the
backward direction (j, k). (i, l) CMC detects circular causality for the $x-y$
pair on distinct frequency channels. In the $x\rightarrow y$ direction the CMC
function indicates causal coupling at $59$ Hz and around $10$ Hz (i).
Meanwhile in the reverse $y\rightarrow x$ direction a backward causal link
appears at $40$ Hz and around $10$ Hz (l). (m-o) Reducing the CMC function to
the maximal peak, CMC values per frequency bands reveal the actual causal
couplings along the frequency axis.

###  2.5 Kuramoto models

We demonstrate how CMC revealed frequency-specific interactions on a small
network of $3$ interconnected Kuramoto oscillators. The network was designed
such that there is a common driver ($z$) with a central frequency of
$10.50422624$ Hz and there are two other oscillators, $x$ and $y$ with
respective central frequencies of $59$ and $40$ Hz, that also drive each
other.

For the $z$-$x$ and $z$-$y$ pairs, CMC showed maximum coherence at around $10$
Hz with $z$ as the driver (Fig. 8). For the $x$-$y$ pair the $x\rightarrow y$
link had a peak at $60$ Hz and $10$ Hz, whereas the opposite direction had a
peak at $40$ Hz and $10$ Hz.

We found that CMC could recover each causal connection properly in a
frequency-specific way, even the circular coupling in different frequency
bands, but it could not recover the proper effect delay with certainty. Due to
the strong periodicity of these oscillators, both CCM and CMC had multiple
peaks of similar amplitude that were periodic in the time shifts between the
two time series. Thus, the proper time delay of the causal effects could not
be inferred. However, this did not affect the inference of the directionality
and the frequency dependency of the causal effects.

###  2.6 Wilson-Cowan model

Figure 9:  Causal discovery with CCM and CMC on a Wilson-Cowan model of
cortical areas V1 and V4. (a) Schematic illustration of the network. Each area
is modeled by $2\times 2$ stochastic differential equations. The coupling
strengths are chosen to be biologically appropriate. (b) Time-delayed CCM
analysis of the Wilson-Cowen dataset ($E=9,\tau=1$). The solid curves show the
average performance over reconstructed state coordinates and realizations
($N=10$) in the function of temporal shift, whereas dashed lines denotes the
score for the first and the last coordinates of the embeddings, with the
transition between these latter two marked as shading giving a sense about the
temporal uncertainty of effect propagation delay. The $V4\rightarrow V1$
connection is stronger at the $-8$ timesteps shift (orange). (c, e) Shift-
frequency CMC analysis shows a strong peak at the $5-20$ Hz range for the
$V4\rightarrow V1$ connection on (e) and the corresponding Granger-peak on the
anti-causal side on (c). (d, f) The same peak is present on the band-
normalized CMC functions on (f), but an additional peak becomes visible in the
opposite direction ($V1\rightarrow V4$) at the $20-50$ Hz range on (d). (g, h)
CMC prominences reveal frequency-specific causal interactions between V1 and
V4 both for the raw (g) and the normalized (h) case. From V1 to V4, the causal
effect appears in the $20-50$ Hz range, whereas the V4 to V1 coupling takes
shape in the $1-20$ Hz frequency band.

We applied CMC on the Wilson-Cowan model of a pair of cortical areas, V1 and
V4. The model is a system of stochastic differential equations, that simulates
mean-field cortical dynamics as coupled supragranular and infragranular
excitatory and inhibitory cell populations, $2\times 2$ in each area. The
areas are linked together with a biologically realistic connection structure
(Fig. 9a, [61]).

CMC revealed frequency-dependent interactions in the V1-V4 model (Fig. 9). The
basic time-delayed CCM analysis returned a pronounced effect from V4 to V1,
without information on the frequency of the connections (Fig. 9b). The time
delay of the effect is around $-6$ time-steps, that is at $32$ ms. CMC found
frequency-dependent connections, a feedforward V1$\rightarrow$V4 link in the
$25-50$ Hz band, and a feedback V4$\rightarrow$V1 link in the $1-20$ Hz band
(Fig. 9 c-h). This directionality in the coupling showed up if we took the
maximal prominence value of coherence over the time-shift axis(Fig. 9c, e, g),
and became even more visible if we normalized the coherence values for each
frequency band beforehand (Fig. 9d, f, h). The time shifts of the maximal
prominence values were in the $-2$ and $-8$ shifts for the V1$\rightarrow$V4
and V4$\rightarrow$V1 links respectively for the unnormalized CMC function,
and at the $-1$ and $-6$ shifts for the normalized CMC analysis. Both results
show that the V1$\rightarrow$V4 feedforward connection has smaller effect
propagation delay.

##  3 Discussion

Cross-Mapping Coherence (CMC) reveals frequency-specific interactions, but
there are several ways to form a frequency-dependent method from Convergent
Cross-Mapping (CCM) by making modifications at various points in the
algorithm, such as the embedding, the prediction, or the evaluation steps.
Here we chose to change the evaluation metrics from correlation to coherence,
because it is a minimal working modification to the original algorithm which
achieves a frequency-dependent measure of causality.

We introduced the peak prominence metrics to determine causal effect delays in
the frequency- and shift-dependent CMC functions. The frequency-dependent
extension of the cross-mapping method rendered the evaluation of the delay-
dependent cross-mapping functions a complex task. Thus, visual examination of
the peaks for all frequency bands was not a viable solution in determining the
causal effect delays in each frequency band. To handle this challenge, we
introduced a maximum detection method, based on the prominence measure of the
peaks of the cross-mapping functions, applied in all frequency bands
separately. Our simulations with the Wilson-Cowan model highlight the
importance of this type of peak evaluation and maximum detection. If one would
simply take the maxima of the CMC function on the negative semi-axis (where
cause precedes the consequences), one would get similarly high values in both
directions on all frequency bands. If done this way, the evaluation would
almost certainly return false circular results in the case of a unidirectional
connection. But closer examination reveals that, on many frequencies, the high
values on the negative semi-axis are actually the upward slope of higher peaks
on the positive semi-axis, thus evaluating the predictability of the
consequence, not the reconstructability of the cause (Fig. 9c-f). Our
prominence-based solution retrieved the proper directionality for
unidirectional connections, while not losing the ability to detect
bidirectional connections.

CMC identifies causal relationships in the frequency domain; however, it is
not unique, as other spectral methodologies, such as spectral Granger
Causality (SGC) and Transfer Entropy Spectrum (TES), exhibit analogous
characteristics and functionalities. In particular, both CMC and SGC showed
the same frequency-specific directional coupling in the Wilson-Cowan model
[61, 62]. However, given that SGC and TES are based on the predictive
causality principle, whereas CMC is rooted in the theory of dynamical systems,
CMC introduces a novel theoretical and practical aspect for discovering causal
relationships from time-series data.

Frequency Domain CCM (FDCCM) and Cross-Frequency Symbolic CCM (CFSCCM) are two
alternative extensions of the CCM method that utilize spectral information
[43, 59]. FDCCM uses the Fourier transform at the embedding step, but does not
provide explicit frequency-specific information about the causal link. CFSCCM
investigates the cross-frequency coupling between fast and slow signal modes.
It decomposes the signal into modes, then computes amplitudes and phases by
Hilbert transform, and finally, it computes symbolic CCM between amplitudes
and phases. In contrast, CMC applies a different approach, it runs on the time
domain signals and uses spectral information only at the evaluation step by
the coherence metric.

CMC shows convergence with increasing data length; in the unidirectionally
coupled example, this convergence seemed independent of frequency bands. CMC
may need longer measurements for systems with smoothly varying dynamics to
achieve convergence, especially in lower frequency bands. To ensure the
required number of recurrences, at least 15-20 cycles should occur during the
investigated section to make the coherence computation meaningful.

The CMC method proved very sensitive, detecting even weak couplings. This
feature is inherited from CCM since it was originally developed to detect
causal relations between weakly coupled nonlinear and non-separable dynamical
systems [17].

The CMC method demonstrates robust performance in the presence of
observational noise, maintaining efficacy up to a signal-to-noise ratio of
$10$. Although this degree of noise tolerance may be adequate for specific
applications, information-theoretic methods such as SGC and TES may exhibit
superior performance on noisy datasets, as evidenced by recent time domain
information-theoretical methods [63].

The choice of embedding parameters affects the shape of the CMC function. The
embedding dimension must be large enough to embed the dynamics of the
underlying system. Provided that this condition is satisfied, only the shape
of the CMC function is influenced by the choice of the specific parameter
value within a reasonable range, the conclusion about the direction of causal
influence does not change. The chosen embedding dimension determines the size
of the embedding window, which introduces uncertainty in the temporal labeling
of the reconstructed states. Consequently, this causes uncertainty in the
propagation delay of the causal effect, manifesting itself as a plateau in the
CMC time-shift function. This plateau complicates the precise determination of
the effect propagation delay; hence, we recommend using the minimal embedding
dimension that satisfies the necessary criteria. Additionally, the embedding
delay must be calibrated to the fastest frequency component in the signal, to
mitigate the risk of aliasing.

A current weakness of CMC is that there is no statistical framework for
significant link-detection backed by the theory of nonlinear dynamics. This is
a caveat compared to SGC and TES, where the statistical formulation and
detection are already well-established [12]. There were several previous
attempts to create this statistical framework for topological causal discovery
methods [45, 46, 37, 47, 48], but these tools neglect significant information
about the bundle structure of coupled dynamical systems [23, 20]. Therefore, a
proper statistical test is yet to be created.

Another interesting extension of CMC would be to detect time-varying
connections by swapping simple coherence metrics for wavelet coherence. This
would allow us to perform time-frequency analysis to investigate temporally
varying frequency-specific causal relations between time series.

Finally, within the scope of this study, we introduced the CMC method to
detect causality in the frequency domain using simulated datasets. Further
research should include the application of the proposed method on diverse
data, with the objective of elucidating novel insights into the underlying
mechanisms of nonlinear dynamical systems.

##  4 Conclusions

In this study, we proposed the Cross-Mapping Coherence (CMC) method for
detecting frequency-dependent directional and circular causal interactions
from time series data. We demonstrated the efficacy of CMC on various
nonlinear model-systems, including logistic maps, Lorenz systems, Kuramoto
oscillators, and the Wilson-Cowan model of V1-V4 cortical areas.

The CMC method extends the Convergent Cross-Mapping method by changing
performance evaluation from the coefficient of determination to coherence, and
automates the detection of causal effect delay with peak prominence metrics.

We demonstrated that CMC identifies complex causal interactions between
nonlinear systems. CMC revealed uniform causality along each frequency band
for coupled logistic maps, frequency-dependent interactions with higher power
at lower frequencies for the Lorenz systems, and frequency-specific
interactions in coupled Kuramoto oscillators. Moreover, it identified high
frequency feedforward and low frequency feedback connections in the Wilson-
Cowan model.

We investigated how detection performance depends on sample size, coupling
strength, and observation noise. In addition, we inspected the effects of the
choice of embedding parameters and gave recommendations on how to choose them
optimally.

In summary, our findings reveal that the novel CMC method is efficient at
revealing interactions between nonlinear systems in the frequency domain. Our
results may have important implications for various fields, including
neuroscience, Earth systems science, and engineering, where understanding
complex nonlinear interactions is crucial.

##  5 Methods

The simulations and analyses were carried out using the python numpy
package[64], visualization were generated by matplotlib[65].

###  5.1 Data Simulations

####  5.1.1 Logistic maps

Coupled logistic maps can be generated according to the equations:

$x_{i}\left[t+1\right]\;=\;\mathcal{F}\left(r_{i}\;\
x_{i}\left[t\right]\;\left(1-\sum_{j}A_{ij}\;x_{j}\left[t\right]\right)\right)$
(5)

where $r_{i}$ is the rate-parameter of the $i$th logistic map, $A_{ij}$ is the
coupling matrix with $A_{ii}=1$ and $\mathcal{F}$ is a mirroring boundary
condition to constrain the data-point to the $[0,1]$ interval.

We simulated $n=10000$ datapoints with unidirectional, circular,common driver,
and independent connection topologies. The simulation parameters are shown in
Table 1.

scenario | $r_{1}$ | $r_{2}$ | $r_{3}$ | $A_{21}$ | $A_{12}$ | $A_{13}$ | $A_{23}$
---|---|---|---|---|---|---|---
$\rightarrow$ | $3.9902032398544094$ | $3.9900842430866197$ | $0$ | $0.05$ | $0.0$ | $0.0$ | $0.0$
$\leftrightarrow$ | $3.9903787118484475$ | $3.9900528401775186$ | $0$ | $0.05$ | $0.05$ | $0.0$ | $0.0$
$\Lsh\Rsh$ | $3.9903839016316964$ | $3.9904120926448896$ | $3.9904110403001893$ | $0.0$ | $0.0$ | $0.05$ | $0.05$
$\perp$ | $3.9903770705735107$ | $3.9907504255884914$ | $0$ | $0.0$ | $0.0$ | $0.0$ | $0.0$
Table 1: The values of the rate and coupling parameters of the coupled
logistic maps.

We tested the length dependence of CMC on unidirectionally coupled logistic
map systems on $L=400$, $700$, $1000$, $2000$, $5000$ dataset lengths. The
simulations were carried out with the rate parameters $r_{1}=3.99097965$ and
$r_{2}=3.99024767$ and the coupling parameter was set to $A_{21}=0.05$.

The coupling strength sensitivity was checked on the same kind of coupled
systems by varying the coupling strength parameter $A_{21}$ to take the
values: $0$, $0.05$, $0.1$, $0.15$, and $0.2$ with rate parameters
$r_{1}=r_{2}=3.99$ and data set length $L=2000$.

The robustness for observational noise was investigated in the signal-to-noise
ratios of $1000$, $100$, $10$, $5$, and $2$ by adding Gaussian white noise to
the signal with corresponding standard deviation. The rate parameters were
$r_{1}=3.99097965$ and $r_{2}=3.99024767$, and the coupling parameter was set
to $A_{21}=0.15$ with the data length as $L=2000$.

The effect of the embedding parameters on the inference results was analyzed
by observing the shape of the CMC function depending on the embedding
dimension in the range of [2, 10] with two-wise steps and with a data length
of $L=20000$. The simulation parameters were the same as in the noise
robustness analysis with the noise level set to zero.

####  5.1.2 Lorenz system

We simulated coupled Lorenz systems ($T=5000,\Delta t=10^{-3}$) according to
the equations:

$\begin{split}\dot{x}_{1}&=\sigma_{1}\left[(y_{1}-x_{1})+\kappa_{1}(y_{2}-x_{1})\right]\\\
\dot{y}_{1}&=x_{1}(\rho_{1}-z_{1})-y_{1}\\\
\dot{z}_{1}&=x_{1}y_{1}-\beta_{1}z_{1}\\\ \\\
\dot{x}_{2}&=\sigma_{2}\left[(y_{2}-x_{2})+\kappa_{2}(y_{1}-x_{2})\right]\\\
\dot{y}_{2}&=x_{2}(\rho_{2}-z_{2})-y_{2}\\\
\dot{z}_{2}&=x_{2}y_{2}-\beta_{2}z_{2}\\\ \end{split}$ (6)

where parameters are signed as Greek lowercase letters, and the coupling is
taking effect through the $\kappa_{i}$ terms. The parameter values are shown
in Table 2 and Table 3. The integration was performed using the 4th order
Runge-Kutta method with the scipy odeint function[66].

parameter | value
---|---
$\sigma_{1}$ | $10.000$
$\rho_{1}$ | $27.000$
$\beta_{1}$ | $2.667$
$\sigma_{2}$ | $10.209$
$\rho_{2}$ | $25.900$
$\beta_{2}$ | $2.652$
Table 2: fix parameter values for the Lorenz simulations. | $\rightarrow$ | $\leftrightarrow$ | $\perp$
---|---|---|---
$\kappa_{1}$ | $0.0$ | $0.1$ | $0.0$
$\kappa_{2}$ | $0.1$ | $0.1$ | $0.0$
Table 3: Coupling parameters for the Lorenz simulations

####  5.1.3 Kuramoto

We simulated a small network of Kuramoto oscillators to demonstrate
directional couplings on different frequency bands with simulation step
$\Delta t=5\times 10^{-3}$, $N=20000$:

$\begin{split}\dot{\theta}_{i}&=2\pi\nu_{i}+\frac{K_{i}}{m}\sum_{j=1}^{m}\sin(\theta_{i}-\theta_{j})+\eta_{i}\\\
x_{i}&=\sin\theta_{i}\end{split}$ (7)

where $m=3$, the basic frequencies are $\nu_{i}=[10.50422624,59,40]$ Hz, the
coupling constants are $K_{i}=[0,3,4.3]$ , and
$\eta_{i}\sim\mathcal{N}(0,0.1)$ is a white noise term. The actual observed
time series are the values of $x_{i}$. The numerical integration of the
stochastic differential equation was carried out using the sdeint package[67].

####  5.1.4 Wilson-Cowan model

A nonlinear Wilson-Cowan-type neural mass model was employed to simulated the
coupled cortical areas. The activity of a single area is described by four
uniquely coupled stochastic differential equations, representing a pair of
excitatory and inhibitory populations in two laminar compartments (supra- and
infragranular) of the forms:

$\tau_{E}\frac{dr_{E}}{dt}=-r_{E}+\Theta(I_{E}^{net}+I_{E}^{ext})+\sqrt{\tau_{E}}\xi_{E}(t),$
(8)

$\tau_{I}\frac{dr_{I}}{dt}=-r_{I}+\Theta(I_{I}^{net}+I_{I}^{ext})+\sqrt{\tau_{I}}\xi_{I}(t),$
(9)

where $r_{E,I}$ are dimensionless activities (i.e. mean firing rates) of an
excitatory and an inhibitory population, respectively, $\tau_{E,I}$ are the
corresponding time constants, $\xi_{E,I}$ are Gaussian noise terms with zero
mean, and $\sigma_{E,I}$ strength and $\Theta(x)=x/(1-\exp^{-x})$ is the
transduction function (used instead of a sigmoid function), while
$I_{E,I}^{net}$ and $I_{E,I}^{ext}$ represent the input from other parts of
the network and from exterior sources (e.g from the thalamus, which is not
included as an individual node), respectively. See the supplementary materials
of Mejias et al. (2016) [61] for a detailed description of the original model,
or Varga et al. (2021) [62] for our specific implementation.

###  5.2 Data analysis steps

####  5.2.1 Logistic maps

The logistic map data does not need any preprocessing; we embedded the time
series and ran the CMC analysis in a $[-20,20]$ time delay range with $E=2$,
$\tau=1$ embedding parameters. We plotted the results on the shift-frequency
axis (Fig. 4). We used the same embedding parameter values for the data
length, coupling sensitivity, and noise-robustness analyses.

####  5.2.2 Lorenz systems

For the Lorenz dataset, we subsampled the time series by a factor of $100$ and
then applied the time delay embedding procedure on the $x_{i}$ variables from
the subsystems ($E=7$, $tau=1$). We conducted the CMC analysis in the $[-2,2]$
s range to discover time-delayed interactions (Fig. 7).

####  5.2.3 Kuramoto models

For the Kuramoto dataset, we chose the first 10000 points and ran the
embedding procedure and subsequent CMC analysis ($E=5,\tau=1$) on the dataset
with a $[-0.1,0.1]$ s time-shift range (Fig. 8).

####  5.2.4 Wilson-Cowan model

For the Wilson-Cowan dataset, we embedded the time series from V1 and V4
($E=9,\tau=1$) and conducted the CMC analysis in a $[-0.2,0.2]$ s time-shift
range to detect time-delayed connections over the frequency bands. We carried
out this analysis on $N=10$ instances and then averaged the results. We
applied prominence calculations per frequency band on the causal side of the
timeshift-frequency plot. Also, we averaged over the instances the per-
frequency normalized results to transform out the effect of amplitude
differences on different frequency bands and performed the prominence
reduction technique.

##  6 Acknowledgements

This research was founded by the Hungarian National Research, Development and
Innovation Office K135837 (Z.S.) and the Hungarian Research Network, HUN-REN
supported the project under the grant SA-114-2021 (Z.S.).

## References

  * [1] Judea Pearl  “Causality: Models, reasoning, and inference”  Cambridge University Press, 2000
  * [2] Zsigmond Benkő et al.  “Bayesian inference of causal relations between dynamical systems”  In _Chaos, Solitons & Fractals_ 185, 2024, pp. 115142  DOI: https://doi.org/10.1016/j.chaos.2024.115142
  * [3] Jakob Runge  “Causal network reconstruction from time series: From theoretical assumptions to practical estimation”  In _Chaos_ , 2018  DOI: 10.1063/1.5025050
  * [4] Zsigmond Benkő et al.  “Causal relationship between local field potential and intrinsic optical signal in epileptiform activity in vitro”  In _Scientific Reports_ 9.1, 2019, pp. 5171  DOI: 10.1038/s41598-019-41554-x
  * [5] Marcell Stippinger et al.  “Causal Discovery of Stochastic Dynamical Systems: A Markov Chain Approach”  In _Mathematics_ 11.4, 2023  DOI: 10.3390/math11040852
  * [6] Marcell Stippinger et al.  “CCDH: Complexity based Causal Discovery of Hidden common cause in time series”  In _Chaos, Solitons & Fractals_ 176, 2023, pp. 114054  DOI: https://doi.org/10.1016/j.chaos.2023.114054
  * [7] N. Wiener  “The Theory of Prediction”  In _Modern Mathematics for Engineers_ 58 McGraw-Hill Book Company, Inc., New York, 1956, pp. 323–329
  * [8] CWJ Granger  “Investigating causal relations by econometric models and cross-spectral methods”  In _Econometrica: Journal of the Econometric Society_ 37.3, 1969, pp. 424–438
  * [9] Thomas Schreiber  “Measuring information transfer.” MAG ID: 2041782669  In _Physical Review Letters_ , 2000  DOI: 10.1103/physrevlett.85.461
  * [10] Michael Wibral et al.  “Measuring Information-Transfer Delays” MAG ID: 2016918012  In _PLOS ONE_ , 2013  DOI: 10.1371/journal.pone.0055809
  * [11] John Geweke  “Measurement of Linear Dependence and Feedback between Multiple Time Series”  In _Journal of the American Statistical Association_ 77.378, 1982, pp. 304–313  DOI: 10.1080/01621459.1982.10477803
  * [12] Steven L Bressler and Anil K Seth  “Wiener-Granger causality: a well established methodology.”  In _NeuroImage_ 58.2 Elsevier Inc., 2011, pp. 323–9  DOI: 10.1016/j.neuroimage.2010.02.059
  * [13] Yang Tian, Yaoyuan Wang, Ziyang Zhang and Pei Sun  “Fourier-domain transfer entropy spectrum”  In _Physical Review Research_ 3.4, 2021, pp. L042040  DOI: 10.1103/PhysRevResearch.3.L042040
  * [14] Jakob Runge et al.  “Detecting causal associations in large nonlinear time series datasets”  In _Science Advances_ , 2019  DOI: 10.1126/sciadv.aau4996
  * [15] Alex Tank et al.  “Neural Granger Causality.”  In _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 2021  DOI: 10.1109/tpami.2021.3065601
  * [16] Ricards Marcinkevics and Julia E. Vogt  “Interpretable Models for Granger Causality Using Self-explaining Neural Networks”  In _International Conference on Learning Representations_ , 2021  DOI: null
  * [17] George Sugihara et al.  “Detecting causality in complex ecosystems.”  In _Science (New York, N.Y.)_ 338.6106, 2012, pp. 496–500  DOI: 10.1126/science.1227079
  * [18] Hao Ye, Ethan R. Deyle, Luis J. Gilarranz and George Sugihara  “Distinguishing time-delayed causal interactions using convergent cross mapping”  In _Scientific Reports_ 5, 2015, pp. 14750  DOI: 10.1038/srep14750
  * [19] Yoshito Hirata and Kazuyuki Aihara  “Identifying hidden common causes from bivariate time series: A method using recurrence plots” MAG ID: 2040971434  In _Physical Review E_ , 2010  DOI: 10.1103/physreve.81.016203
  * [20] Zsigmond Benkő et al.  “Complete Inference of Causal Relations between Dynamical Systems” MAG ID: 3006999686  In _arXiv: Quantitative Methods_ , 2020
  * [21] Zsigmond Benkő et al.  “Manifold-adaptive dimension estimation revisited.” MAG ID: 3047731977  In _ArXiv_ , 2020
  * [22] Floris Takens  “Detecting strange attractors in turbulence”  In _Dynamical Systems and Turbulence, Warwick 1980_ 898, 1981, pp. 366–381  DOI: 10.1007/BFb0091924
  * [23] J Stark  “Delay embeddings for foreced systems. I. Deterministic forcing”  In _Journal of Nonlinear Science_ 9, 1999, pp. 255–332  DOI: 10.1007/s003329900072
  * [24] Ethan R Deyle and George Sugihara  “Generalized theorems for nonlinear state space reconstruction.”  In _PloS one_ 6.3, 2011, pp. e18295  DOI: 10.1371/journal.pone.0018295
  * [25] Robert M. May and George F. Oster  “Bifurcations and Dynamic Complexity in Simple Ecological Models”  In _The American Naturalist_ 110.974 [University of Chicago Press, American Society of Naturalists], 1976, pp. 573–599  URL: http://www.jstor.org/stable/2459579
  * [26] Elemer Labos  “Conditions for periodic and aperiodic behavior of formal neural nets”  In _Biophysical Chemistry_ 50, 1994, pp. 217–223  DOI: 10.1016/0301-4622(94)85033-X
  * [27] Edward N Lorenz  “Deterministic nonperiodic flow”  In _Journal of atmospheric sciences_ 20.2, 1963, pp. 130–141
  * [28] O.E. Rössler  “An equation for continuous chaos”  In _Physics Letters A_ 57.5, 1976, pp. 397–398  DOI: https://doi.org/10.1016/0375-9601(76)90101-8
  * [29] Alfred J. Lotka  “ELEMENTS OF PHYSICAL BIOLOGY”  In _Science Progress in the Twentieth Century (1919-1933)_ 21.82 Temporary Publisher, 1926, pp. 341–343  URL: http://www.jstor.org/stable/43430362
  * [30] Vito Volterra  “Variations and Fluctuations of the Number of Individuals in Animal Species living together”  In _ICES Journal of Marine Science_ 3.1, 1928, pp. 3–51  DOI: 10.1093/icesjms/3.1.3
  * [31] Jack D. Cowan, Jeremy Neuman and Wim Drongelen  “Wilson–Cowan Equations for Neocortical Dynamics”  In _The Journal of Mathematical Neuroscience_ 6.1, 2016, pp. 1  DOI: 10.1186/s13408-015-0034-5
  * [32] Steven J. Schiff et al.  “Detecting dynamical interdependence and generalized synchrony through mutual prediction in a neural ensemble”  In _Physical Review E_ 54.6, 1996, pp. 6708–6724  DOI: 10.1103/PhysRevE.54.6708
  * [33] James M. McCracken and Robert S. Weigel  “Convergent cross-mapping and pairwise asymmetric inference”  In _Physical Review E_ 90.6, 2014, pp. 062903  DOI: 10.1103/PhysRevE.90.062903
  * [34] Huanfei Ma, Kazuyuki Aihara and Luonan Chen  “Detecting Causality from Nonlinear Dynamics with Short-term Time Series” MAG ID: 2055273782  In _Scientific Reports_ , 2015  DOI: 10.1038/srep07464
  * [35] Xinlei Ge and Aijing Lin  “Dynamic causality analysis using overlapped sliding windows based on the extended convergent cross-mapping”  In _Nonlinear Dynamics_ 104.2, 2021, pp. 1753–1765  DOI: 10.1007/s11071-021-06362-x
  * [36] Adam Thomas Clark et al.  “Spatial convergent cross mapping to detect causal relationships from short time series” MAG ID: 2025836487  In _Ecology_ , 2015  DOI: 10.1890/14-1479.1
  * [37] Yu Huang, Christian L.. Franzke, Naiming Yuan and Zuntao Fu  “Systematic identification of causal relations in high-dimensional chaotic systems: application to stratosphere-troposphere coupling”  In _Climate Dynamics_ 55.9-10, 2020, pp. 2469–2481  DOI: 10.1007/s00382-020-05394-0
  * [38] Jinxian Deng et al.  “Causalized convergent cross-mapping and its approximate equivalence with directed information in causality analysis”  In _PNAS Nexus_ 3.1, 2023, pp. pgad422  DOI: 10.1093/pnasnexus/pgad422
  * [39] Guanchao Feng et al.  “Improving Convergent Cross Mapping for Causal Discovery with Gaussian Processes”  In _ICASSP 2020 - 2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ Barcelona, Spain: IEEE, 2020, pp. 3692–3696  DOI: 10.1109/ICASSP40776.2020.9053166
  * [40] Guanchao Feng, J. Quirk and Petar M. Djuric  “Detecting Causality using Deep Gaussian Processes”  In _2019 53rd Asilomar Conference on Signals, Systems, and Computers_ Pacific Grove, CA, USA: IEEE, 2019, pp. 472–476  DOI: 10.1109/IEEECONF44664.2019.9048963
  * [41] Yu Huang, Zuntao Fu and Christian L.. Franzke  “Detecting causality from time series in a machine learning framework”  In _Chaos: An Interdisciplinary Journal of Nonlinear Science_ 30.6, 2020, pp. 063116  DOI: 10.1063/5.0007670
  * [42] Edward De Brouwer and Adam Arany  “LATENT CONVERGENT CROSS MAPPING”, 2021
  * [43] Sandeep Avvaru and Keshab K. Parhi  “Effective Brain Connectivity Extraction by Frequency-Domain Convergent Cross-Mapping (FDCCM) and its Application in Parkinson’s Disease Classification”  In _IEEE Transactions on Biomedical Engineering_ , 2023, pp. 1–11  DOI: 10.1109/TBME.2023.3250355
  * [44] Dan Mønster et al.  “Causal inference from noisy time-series data — Testing the Convergent Cross-Mapping algorithm in the presence of noise and external influence”  In _Future Generation Computer Systems_ 73, 2017, pp. 52–62  DOI: 10.1016/j.future.2016.12.009
  * [45] Yalin Wang et al.  “Improved CCM for variable causality detection in complex systems”  In _Control Engineering Practice_ 83, 2019, pp. 67–82  DOI: 10.1016/j.conengprac.2018.10.005
  * [46] P V Krishna and Arun K. Tangirala  “Inferring Direct Causality from Noisy Data using Convergent Cross Mapping”  In _2019 58th Annual Conference of the Society of Instrument and Control Engineers of Japan (SICE)_ Hiroshima, Japan: IEEE, 2019, pp. 1523–1528  DOI: 10.23919/SICE.2019.8859922
  * [47] Sarah Cobey and Edward B. Baskerville  “Limits to causal inference with state-space reconstruction for infectious disease.” MAG ID: 2230868847  In _PLOS ONE_ , 2016  DOI: 10.1371/journal.pone.0169050
  * [48] Anna Krakovská and Filip Hanzely  “Testing for causality in reconstructed state spaces by an optimized mixed prediction method”  In _Physical Review E_ 94.5, 2016, pp. 052203  DOI: 10.1103/PhysRevE.94.052203
  * [49] Xinlei Ge and Aijing Lin  “Symbolic convergent cross mapping based on permutation mutual information”  In _Chaos, Solitons & Fractals_ 167, 2023, pp. 112992  DOI: 10.1016/j.chaos.2022.112992
  * [50] Leo Breston et al.  “Convergent cross sorting for estimating dynamic coupling”  In _Scientific Reports_ 11.1, 2021, pp. 20374  DOI: 10.1038/s41598-021-98864-2
  * [51] Ameer Ghouse, Luca Faes and Gaetano Valenza  “Inferring directionality of coupled dynamical systems using Gaussian process priors: Application on neurovascular systems”  In _Physical Review E_ 104.6, 2021, pp. 064208  DOI: 10.1103/PhysRevE.104.064208
  * [52] Ameer Ghouse and Gaetano Valenza  “Inferring Parsimonious Coupling Statistics in Nonlinear Dynamics with Variational Gaussian Processes” arXiv:2203.03868 [stat]  arXiv, 2022  URL: http://arxiv.org/abs/2203.03868
  * [53] Matthew J. Vowels, Necati Cihan Camgoz and Richard Bowden  “Shadow-Mapping for Unsupervised Neural Causal Discovery” MAG ID: 3188528620, 2021, pp. 1740–1743
  * [54] Emiliano Díaz et al.  “Inferring causal relations from observational long-term carbon and water fluxes records”  In _Scientific Reports_ 12.1, 2022, pp. 1610  DOI: 10.1038/s41598-022-05377-7
  * [55] Ying, Xiong et al.  “Continuity scaling: A rigorous framework for detecting and quantifying\n causality accurately” MAG ID: 4221144396, 2022
  * [56] Siyang Leng et al.  “Partial cross mapping eliminates indirect causal influences”  In _Nature Communications_ 11.1, 2020, pp. 2632  DOI: 10.1038/s41467-020-16238-0
  * [57] Niels Berkel et al.  “Modeling interaction as a complex system”  In _Human–Computer Interaction_ 36.4, 2021, pp. 279–305  DOI: 10.1080/07370024.2020.1715221
  * [58] S Nithya and Arun K Tangirala  “Multivariable Causal Analysis of Nonlinear Dynamical Systems using Convergent Cross Mapping”  In _2021 Seventh Indian Control Conference (ICC)_ IEEE, 2021, pp. 436–441
  * [59] Danlei Gu, Aijing Lin and Guancen Lin  “Detection of Attention Deficit Hyperactivity Disorder in children using CEEMDAN-based cross frequency symbolic convergent cross mapping”  In _Expert Systems with Applications_ 226, 2023, pp. 120105  DOI: 10.1016/j.eswa.2023.120105
  * [60] André M. Bastos and Jan-Mathijs Schoffelen  “A Tutorial Review of Functional Connectivity Analysis Methods and Their Interpretational Pitfalls”  In _Frontiers in Systems Neuroscience_ 9, 2016  DOI: 10.3389/fnsys.2015.00175
  * [61] Jorge F. Mejias, John D. Murray, Henry Kennedy and Xiao-Jing Wang  “Feedforward and feedback frequency-dependent interactions in a large-scale laminar network of the primate cortex”  In _Science Advances_ 2.11 American Association for the Advancement of Science, 2016  DOI: 10.1126/sciadv.1601335
  * [62] Bálint Varga et al.  “Network Path Convergence Shapes Low-Level Processing in the Visual Cortex”  In _Frontiers in Systems Neuroscience_ 15, 2021, pp. 645709  DOI: 10.3389/fnsys.2021.645709
  * [63] Diego Bueso, Maria Piles and Gustau Camps-Valls  “Explicit Granger causality in kernel Hilbert spaces”  In _Phys. Rev. E_ 102 American Physical Society, 2020, pp. 062201  DOI: 10.1103/PhysRevE.102.062201
  * [64] Charles R. Harris et al.  “Array programming with NumPy”  In _Nature_ 585.7825 Springer ScienceBusiness Media LLC, 2020, pp. 357–362  DOI: 10.1038/s41586-020-2649-2
  * [65] J.. Hunter  “Matplotlib: A 2D graphics environment”  In _Computing in Science & Engineering_ 9.3 IEEE COMPUTER SOC, 2007, pp. 90–95  DOI: 10.1109/MCSE.2007.55
  * [66] Pauli Virtanen et al.  “SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python”  In _Nature Methods_ 17, 2020, pp. 261–272  DOI: 10.1038/s41592-019-0686-2
  * [67] Matthew Aburn  “Critical fluctuations and coupling of stochastic neural mass models”, 2017  DOI: 10.14264/uql.2017.148