# Data-driven operator-theoretic analysis of weak interactions in synchronized
network dynamics

Yuka Hashimoto1, Masahiro Ikeda2,3, Hiroya Nakao4, and Yoshinobu Kawahara5,2
1\. NTT Network Service Systems Laboratories, NTT Corporation, Tokyo 180-8585,
Japan
2\. Center for Advanced Intelligence Project, RIKEN, Tokyo 103-0027, Japan
3\. Faculty of Science and Technology, Keio University, Yokohama 223-8522,
Japan
4\. School of Engineering, Tokyo Institute of Technology, Tokyo 152-8552,
Japan
5\. Institute of Mathematics for Industry, Kyushu University, Fukuoka
819-0395, Japan

## Abstract

We propose a data-driven approach to extracting interactions among oscillators
in synchronized networks. The phase model describing the network is estimated
directly from time-series data by solving a multiparameter eigenvalue problem
associated with the Koopman operator on vector-valued function spaces. The
asymptotic phase function of the oscillator and phase coupling functions are
simultaneously obtained within the same framework without prior knowledge of
the network. We validate the proposed method by numerical simulations and
analyze real-world data of synchronized networks.

#### Introduction.—

Network dynamics is important in various phenomena in nature and human society
in which multiple components interact, including economic network, traffic
network, power grids, neural circuits, and physiological systems. Various
model-based approaches have been developed to analyze network dynamics. One
relevant method is the phase model for a network of coupled oscillators, which
gives a low-dimensional description of interacting oscillators by reducing the
dimensionality of individual oscillators dynamics [26, 7, 1, 36, 37, 47, 39].
The phase model is represented by intrinsic frequencies of the oscillators and
phase coupling functions (PCF), which respectively characterize the dynamics
of individual oscillators and interactions between the oscillators.

Data-driven approaches to dynamical systems have been widely investigated
recently [46, 13, 49, 34, 6, 5, 45, 25, 43]. They enable us to understand the
dynamics using given time-series data without knowing the underlying dynamical
system that generates the data. In particular, the application of Koopman
operator theory for this purpose has attracted much attention [40, 27, 8, 48,
21, 44, 9, 19, 12, 23, 16, 18]. The Koopman operator gives a lifted
representation of a nonlinear dynamical system. Since the Koopman operator is
linear even though the original system is nonlinear, we can use linear algebra
to analyze the behaviors of the dynamics with given data. One successful
application of Koopman operator theory to data analysis is dynamic mode
decomposition (DMD) and extended DMD (EDMD) [40, 48, 21]. We can decompose
observables of the system by using the eigenvalues and eigenfunctions of the
Koopman operator and extract the long-term behavior of the system. Using DMD
and EDMD, we can extract the Koopman eigenvalues and eigenfunctions from time-
series data and predict long-term behaviors. Koopman operator theory has also
been applied to network dynamics. However, most existing works do not focus on
the data-driven setting or require partial information about the dynamical
system [31, 29, 30, 17].

The relationship between Koopman operator theory and phase models has been
investigated [42, 29]. For example, the intrinsic frequency and asymptotic
phase function of a limit-cycle oscillator are obtained from the eigenvalue
corresponding to the fundamental frequency and associated eigenfunction of the
Koopman operator, respectively. Therefore, by estimating the Koopman operator
and computing those eigenvalues and eigenfunctions, we can reconstruct the
phase model only with given data. However, most existing studies focus on
extracting the phase function, or equivalently the eigenfunction of the
Koopman operator from the data. There is no existing work on the simultaneous
reconstruction of the PCF of coupled oscillators in the same framework via the
Koopman operator for network dynamics.

In this paper, we propose a data-driven approach to extract the phase model
based on Koopman operators. We generalize the relationship between the Koopman
operator and the phase model for a single isolated oscillator to that for a
network of oscillators with interactions. To describe the dynamics of each
dynamical element in a network, we focus on the Koopman operator defined on a
vector-valued function space (instead of that on a scalar-valued function
space for describing a single isolated oscillator). We show that the phase
model with interactions is obtained by solving a generalized multiparameter
eigenvalue problem (GMEP) associated with the Koopman operator, instead of the
standard eigenvalue problem discussed in the existing studies. By estimating
the Koopman operator and the solution of the GMEP from given data, we can
obtain the common intrinsic frequency, phase function, and PCF between the
oscillators and reconstruct a phase model describing the network dynamics that
generate the data. We validate the proposed method using numerical examples
and apply it to real-world healthcare data.

#### Reduction of weakly interacting network dynamics.—

Figure 1: Reduction of coupled nonlinear oscillators to the phase model.

Let $\mathcal{X}=\mathbb{R}^{d}$ and $X_{1}(t),\ldots,X_{N}(t)\in\mathcal{X}$
be mutually interacting $N$ oscillatory elements:

$\displaystyle\frac{dX_{i}(t)}{dt}=F_{i}(X_{i}(t))+$
$\displaystyle\sum_{k=1}^{N}G_{i,k}(X_{i}(t),X_{k}(t))$ (1)

for $i=1,\ldots,N$, where $F_{i}\colon\mathcal{X}\to\mathcal{X}$ and
$G_{i,k}\colon\mathcal{X}\times\mathcal{X}\to\mathcal{X}$ represent the
individual dynamics of the element $i$ and effects of the element $k$ on the
element $i$, respectively. It is often the case that the intrinsic frequencies
of the elements are close to each other [37, 39]. Thus, we here assume that
the elements $1,\ldots,N$ are dominated by a common intrinsic frequency
$\omega$ and weakly interact. That is, we consider the situation where the
dynamics can be reduced to the following phase model:

$\displaystyle\frac{d\theta_{i}(t)}{dt}=\omega+$
$\displaystyle\Gamma_{i}(\theta_{i}(t)-\theta_{1}(t),\ldots,\theta_{i}(t)-\theta_{N}(t))$
(2)

for $i=1,\ldots,N$. Here, $\theta_{i}(t)\in[0,2\pi)$ is the phase
corresponding to $X_{i}(t)$ (i.e., the variable obtained from $X_{i}(t)$ via
some function $\mathcal{X}\to[0,2\pi)$) and $\Gamma_{i}$ is the PCF
representing effects of the phase differences on the element $i$ [35]. Note
that since the interactions are assumed to be weak, we can approximate the
effect only with the phase difference.

The phase model enables us to perform a detailed analysis of synchronization
due to interactions among elements in a network. To analyze synchronization
dynamics only with given times-series data, estimation of the phase model (2)
from data has been discussed [25, 43, 5, 45]. In these existing methods, to
obtain the phase model (2), the phase variable $\theta_{i}$ is first estimated
by, for example, the linear interpolation or the Hilbert transform [25, 43].
Then, the PCF $\Gamma_{i}$ is estimated, for example, by fitting the data to
the Fourier series [5, 45], by using kernel density estimation [25], or by
using dynamic Bayesian inference [43].

#### Data-driven estimation of the phase model via Koopman operator.—

While the traditional geometric perspective of dynamical systems describes the
topological organization of trajectories of the systems, the Koopman operator
gives another perspective of dynamical systems on the evolution of
observables. This perspective was originally introduced by Koopman [24] and
later extended to dissipative systems by Mezic [32]. We here formally discuss
a generic way of deriving the transformation of $X_{i}(t)$ to $\theta_{i}(t)$,
frequency $\omega$, and PCF $\Gamma_{i}$ in Eq. (2) using the Koopman operator
in a vector-valued function space.

We discretize the continuous time as $t_{0},t_{1},\ldots$ with time interval
$\Delta t$, that is, $t_{i}=t_{i-1}+\Delta t$ for $i=1,2,\ldots$. Let
$\mathcal{H}$ be a Hilbert space consisting of $N$-tuples of functions
$\mathbf{v}=[v_{1},\ldots,v_{N}]$, where $v_{i}$ is a complex-valued function
on $\mathcal{X}$. We can regard $\mathbf{v}$ as a vector-valued function that
maps $\mathbf{X}(t)=[X_{1}(t),\ldots,X_{N}(t)]$ to
$[v_{1}(X_{1}(t)),\ldots,v_{N}(X_{N}(t))]$. Then, the Koopman operator $K$ on
$\mathcal{H}$ with respect to the dynamical system (1) is a linear operator on
$\mathcal{H}$ satisfying

$\displaystyle K[v_{1}(X_{1}(t)),\ldots,v_{N}(X_{N}(t))]$
$\displaystyle\qquad=[{v}_{1}(X_{1}(t+\Delta t)),\ldots,v_{N}(X_{N}(t+\Delta
t))]$ (3)

for $[v_{1},\ldots,v_{N}]\in\mathcal{H}$.

First, we estimate the Koopman operator $K$. Let
$[x_{1,0},\ldots,x_{N,0}]=\mathbf{x}_{0},\ldots,[x_{1,T},\ldots,x_{N,T}]=\mathbf{x}_{T}\in\mathcal{X}^{N}$
be given data generated by the dynamical system (1) at time
$t_{0},\ldots,t_{T}$. We assume that $S$ sequences
$\\{\mathbf{x}^{1}_{0},\ldots,\mathbf{x}^{1}_{T}\\},\ldots,\\{\mathbf{x}^{S}_{0},\ldots,\mathbf{x}^{S}_{T}\\}$
of data are given. For the Hilbert space $\mathcal{H}$, we use a vector-valued
reproducing kernel Hilbert space [20]. We generate a subspace of $\mathcal{H}$
by the given data and estimate the Koopman operator $K$ on the subspace. The
estimation of the Koopman operator is proposed in, e.g. [15] and detailed in
Section E.

Next, we drive the phase model via the estimated Koopman operator. The
relationship between the eigenvalues and eigenfunctions of the Koopman
operator and the phase model has already been studied for $N=1$ [42, 29, 30].
Here, we generalize the discussion of the case of $N=1$ to that of $N>1$. We
consider the following generalized multiparameter eigenvalue problem (GMEP):
Find a sequence
$\\{(\lambda_{j},\\{a_{i,k}^{j}\\}_{i,k=1}^{N},\mathbf{u}^{j}=[u_{1}^{j},\ldots,u_{N}^{j}])\\}_{j=1}^{M}$
satisfying

$\displaystyle(\lambda_{1})^{-1}K\mathbf{u}^{1}\odot\cdots\odot(\lambda_{M})^{-1}K\mathbf{u}^{M}$
$\displaystyle=\bigg{(}\sum_{i,k=1}^{N}a^{1}_{i,k}B_{i,k}\mathbf{u}^{1}\bigg{)}\odot\cdots\odot\bigg{(}\sum_{i,k=1}^{N}a^{M}_{i,k}B_{i,k}\mathbf{u}^{M}\bigg{)}.$
(4)

Here, $B_{i,k}$ is a linear operator on $\mathcal{H}$ that maps
$[u_{1},\ldots,u_{k}]\in\mathcal{H}$ to the vector whose $i$th element is
$u_{k}$ and other elements are $0$. Moreover, for vector-valued functions
$\mathbf{u}=[u_{1},\ldots,u_{N}]$ and $\mathbf{v}=[v_{1},\ldots,v_{N}]$,
$\mathbf{u}\odot\mathbf{v}$ is defined as the element-wise product
$\mathbf{u}\odot\mathbf{v}=[u_{1}v_{1},\ldots,u_{N}v_{N}]$. See Section A for
further information about GMEP. Assume that there exists $\omega\in[0,2\pi)$
such that $\lambda_{j}=e^{\sqrt{-1}j\omega\Delta t}$. Then, by the definitions
of the Koopman operator and the operator $B_{i,k}$, we obtain

$\displaystyle\prod_{j=1}^{M}u^{j}_{i}(X_{i}(t+\Delta t))$
$\displaystyle\qquad=\prod_{j=1}^{M}\bigg{(}e^{\sqrt{-1}j\omega\Delta
t}\sum_{k=1}^{N}a_{i,k}^{j}u^{j}_{k}(X_{k}(t))\bigg{)}.$ (5)

Let $\theta^{j}_{i}(t)=\operatorname{arg}({u}^{j}_{i}(X_{i}(t)))$. The map
$\operatorname{arg}(u_{i}^{j}(\cdot))$ transforms $X_{i}(t)$ into the
corresponding phase $\theta_{i}(t)$. Assume in addition for $j=1,\ldots,M$,
$a_{i,k}^{j}\approx 1$ for $i=k$, $a_{i,k}^{j}\approx 0$ for $i\neq k$, that
is, the interactions are weak. Moreover, we normalize the eigenfunctions
$\mathbf{u}_{i}^{j}$ to equalize the effects of all the dynamical elements on
the phase model (see Section B). Then, we regard
$|{u}^{j}_{i}(X_{i}(t))|\approx 1$ for any $i=1,\ldots,N,\ j=1,\ldots,M$ and
any $t$. In addition, since $\lambda_{j}$ is represented as
$\lambda_{j}=e^{\sqrt{-1}j\omega\Delta t}$, we regard $\theta^{j}_{i}\approx
j\theta_{i}$, where $\theta_{i}=\theta^{1}_{i}$ (see Section B). Thus, we have

$\displaystyle\theta_{i}(t+\Delta t)\approx\omega\Delta t+\theta_{i}(t)$
$\displaystyle\qquad+\sum_{j=1}^{M}\operatorname{arg}\bigg{(}\sum_{k=1}^{N}{a}_{i,k}^{j}e^{\sqrt{-1}j(\theta_{k}(t)-\theta_{i}(t))}\bigg{)}.$
(6)

Therefore, we obtain a discretized version of the phase model (2) by setting

$\displaystyle{\Gamma}_{i}(\psi_{1},\ldots,\psi_{N})=\frac{1}{\Delta
t}\sum_{j=1}^{M}\operatorname{arg}\bigg{(}\sum_{k=1}^{N}{a}_{i,k}^{j}e^{-\sqrt{-1}j\psi_{k}}\bigg{)}.$
(7)

Note that $a_{i,k}^{j}$ for $j=1,\ldots,M$ describes the effect of the element
$k$ on the element $i$. See Section B for the details of solving the problem
(4) and Section C for the derivation of Eq. (6).

#### Estimation of phase coupling functions.—

A |  B |  C
---|---|---

D $M=1$ |  E $M=2$ |  F $M=3$ |  G $M=4$
---|---|---|---

Figure 2: The synthetic data generated from the Stuart–Landau model. (A) The
eigenvalues of the estimated Koopman operator in the complex plane (horizontal
axis: real, vertical axis: imaginary). (B–C) The phase function obtained from
the original model (B) and the estimated Koopman operator (C). The black line
represents the limit cycle. (D–G) The difference of the PCF $\Gamma_{d}(\psi)$
obtained from the original model (blue) and estimated Koopman operator
(orange) for $M=1,2,3,4$.

A |  B |  C
---|---|---

D $M=1$ |  E $M=2$ |  F $M=3$ |  G $M=4$
---|---|---|---

Figure 3: The synthetic data generated from the FitzHugh–Nagumo model. (A) The eigenvalues of the estimated Koopman operator in the complex plane. (B–C) The phase function obtained by the model-based approach [28] (B) and the estimated Koopman operator (C). (D–G) The difference of the PCF $\Gamma_{d}(\psi)$ obtained by the model-based approach (blue) and estimated Koopman operator (orange) for $M=1,2,3,4$. A |  B |  C
---|---|---
D |  E |  F
Figure 4: Real-world healthcare data of breathing. (A) Time series regarding
the RIP (blue) and airflow (orange) of a person while sleeping. (B) The
standard eigenvalues of the estimated Koopman operator in the complex plane.
(C) The difference of the PCF $\Gamma_{d}(\psi)$ obtained by the estimated
Koopman operator for $M=2$. (D–F) The phase $\theta_{1}(t)$ w.r.t. RIP data
(D), the phase $\theta_{2}(t)$ w.r.t. airflow (E), and the phase difference
$\theta_{1}(t)-\theta_{2}(t)$ obtained by the estimated Koopman operator.

By solving the GMEP with respect to the estimated Koopman operator, we can
estimate the PCFs. Here, we consider a network with two elements ($N=2$) in
$\mathcal{X}=\mathbb{R}^{2}$ whose dynamics are respectively described for
$i,j\in\\{1,2\\}$ by

$\displaystyle\frac{dX_{i}}{dt}$
$\displaystyle=F_{i}(X_{i})+G_{i,j}(X_{i},X_{j}).$ (8)

For the first example, we consider the Stuart–Landau model, where
$X_{i}(t)=[y_{i}(t),z_{i}(t)]\in\mathbb{R}^{2}$,
$F_{i}(X_{i})=F(X_{i})=[y_{i}-az_{i}-(y_{i}^{2}+z_{i}^{2})(y_{i}-bz_{i}),ay_{i}+z_{i}-(y_{i}^{2}+z_{i}^{2})(by_{i}+z_{i})]$,
$G_{i,j}(X_{i},X_{j})=G(X_{i},X_{j})=[\epsilon\tilde{G}(z_{i},z_{j}),0]$,
$a=2$, $b=1$, and $\epsilon=0.01$. The Stuart–Landau model is a simple model
of oscillators and a normal form of the supercritical Hopf bifurcation. We set
$\tilde{G}(z_{i},z_{j})=z_{j}-z_{i}$. We generated data according to Eq. (8)
and estimated the Koopman operator. We first computed standard eigenvalues of
the Koopman operator, which are illustrated in Fig. 2 A. There are pairs of
eigenvalues $\lambda^{j}_{1},\lambda^{j}_{2}\in\mathbb{C}$ that satisfy
$|\lambda_{i}^{j}|\approx 1$ and $\lambda_{i}^{j}\approx
e^{\sqrt{-1}j\omega_{0}\Delta t}$ for $\omega_{0}\approx 0.25$ and for $i=1,2$
and $j=1,2,\ldots$, which are colored in red in Fig. 2 A. This corresponds to
the fact that two dynamical elements are synchronized and the common frequency
is around $0.25/\Delta t$ (see Subsection B.2 for further details). In fact,
if $\epsilon=0$, then the frequency of each element is shown to be $a-b=1$. By
using the above eigenvalues and the corresponding eigenfunctions, we solve the
GMEP and obtain Eq. (6). See Sections A and G for details of the computations.

To analyze the time evolution of the phase difference between two dynamical
elements, the difference of the PCF is considered [35]. For Eq. (2), it is
defined as $\Gamma_{d}(\psi)={\Gamma}_{1}(-\psi)-{\Gamma}_{2}(\psi)$. Then,
the phase difference $\psi=\theta_{1}-\theta_{2}$ obeys

$\psi(t+\Delta t)\approx\psi(t)+\Delta t\Gamma_{d}(\psi(t)).$ (9)

For the Stuart–Landau model, we can obtain a phase function $\theta_{i}$ and a
PCF $\Gamma_{i}$ for $i=1,\ldots,N$ analytically from the original model [26,
47]. In Figs. 2 B and C, we show the phase functions globally, which are
obtained from the original model and the estimated Koopman operator,
respectively. Moreover, Figs. 2 D–G show the difference of the PCF
$\Gamma_{d}$ obtained by the model-based approach and the estimated Koopman
operator. Regarding the estimation of the function $\Gamma_{d}$ with the
Koopman operator, we set the value $M$ in Eq. (4) as $1\sim 4$. Since the true
PCF of the Stuart–Landau model obtained analytically is represented only with
a sine function, setting $M$ as $1$ is sufficient for obtaining an accurate
estimation. Therefore, the estimated function does not change even if $M$ is
set as a larger value than $1$ in this example.

For the second example, we consider the FitzHugh–Nagumo model, whose orbit and
interactions are more complex than the first example, where
$X_{i}(t)=[y_{i}(t),z_{i}(t)]$, $y_{i},z_{i}\in\mathbb{R}$,
$F_{i}(X_{i})=F(X_{i})=[y_{i}(y_{i}-c)(1-y_{i})-z_{i},\mu^{-1}(y_{i}-dz_{i})]$,
$G_{i,j}(X_{i},X_{j})=G(X_{i},X_{j})=[\epsilon\tilde{G}(z_{i},z_{j}),0]$,
$c=-0.1$, $d=0.5$, $\mu=10$, and $\epsilon=0.0025$. We set
$\tilde{G}(z_{i},z_{j})=z_{i}-z_{j}$, for which the PCF has higher-harmonic
components and the estimation is challenging. We estimated the Koopman
operator and obtained the phase function and difference of the PCF in the same
manner as the first example. The results are illustrated in Figs. 3 D–G. In
Figs. 3 B and C, we show the phase functions globally, which are obtained from
the original model and the estimated Koopman operator, respectively. As in the
first example, regarding the estimation of the function $\Gamma_{d}$ with the
Koopman operator, we set the value $M$ in Eq. (4) as $1\sim 4$. The true PCF
of the FitzHugh–Nagumo model is more complex than that of the Stuart–Landau
model. Therefore, unlike the first example, the estimated function changes and
becomes more accurate as $M$ becomes larger in this example.

For the third example, we consider real-world healthcare data of sleep
patterns [51, 50]. Fig. 4 A illustrates time series regarding the abdominal
respiratory inductance plethysmography (RIP) and airflow of a person while
sleeping. We focus on the synchronization and interaction of RIP and airflow
data. We set $y_{1}$ as a variable describing RIP data and $y_{2}$ as that
describing airflow data. In addition, we set $z_{i}(t)=y_{i}(t+1)$ and
$X_{i}(t)=[y_{i}(t),z_{i}(t)]$ for $i=1,2$. We estimated the Koopman operator
with the data and obtained the phase function and difference of the PCF
$\Gamma_{d}$ in the same manner as the first and second examples. The result
is illustrated in Fig. 4 C. There are two values, indicated by red points in
Fig. 4 C, of $\psi$ which satisfies $\Gamma_{d}(\psi)=0$ and
$\frac{d\Gamma_{d}(\psi)}{d\psi}<0$. These values are stable fixed points of
the phase difference $\psi$. Thus, $\psi$ stays around one of the fixed
points, and because of the effect of noise, $\psi$ sometimes slips to the
other fixed point. Indeed, Fig. 4 F shows the value of phase difference
$\psi=\theta_{1}-\theta_{2}$. The orange lines indicates the values of two
stable fixed points. We can see $\psi$ stays around two different values most
of the time.

#### Conclusion and discussion.—

In this paper, we proposed a data-driven approach to analyze network dynamics
by using Koopman operators on vector-valued function spaces. It was shown that
the phase model for network dynamics can be obtained by solving a GMEP
associated with the Koopman operator. To obtain solutions of the GMEP, we
estimate the Koopman operator with observed time-series data. This enables us
to understand synchronization and interactions in network dynamics. We applied
our proposed method to real-world healthcare data and cohabiting animals data
to analyze the synchronization and interactions.

This paper is the first step towards estimating the interactions described by
the phase model with operator-theoretic approaches. Further theoretical and
practical investigations are required. For example, the convergence of
solutions of the GMEP, noise robustness, dependence on data length, efficient
computations, and designs of the function space should be investigated.

We would like to thank Dr. Isao Ishikawa for constructive discussions. This
work was partially supported by JST CREST Grant Number JPMJCR1913.

## References

  * [1] J. A. Acebrón, L. L. Bonilla, C. J. Pérez Vicente, F. Ritort, and R. Spigler.  The kuramoto model: A simple paradigm for synchronization phenomena.  Rev. Mod. Phys., 77:137–185, 2005.
  * [2] F. V. Atkinson.  Multiparameter Eigenvalue Problems.  Mathematics in Science and Engineering v. 82. Academic Press, 1972.
  * [3] F. V. Atkinson and A. B. Mingarelli.  Multiparameter Eigenvalue Problems: Sturm-Liouville Theory.  CRC Press, 2011.
  * [4] P. Binding and P. J. Browne.  A variational approach to multiparameter eigenvalue problems for matrices.  SIAM J Math anal, 8(5):763–777, 1977.
  * [5] K. A. Blaha, A. Pikovsky, M. Rosenblum, M. T. Clark, C. G. Rusin, and J. L. Hudson.  Reconstruction of two-dimensional phase dynamics from experiments on coupled oscillators.  Phys Rev E, 84:046201–046201, 2011.
  * [6] W. Bomela, S. Wang, C.-A. Chou, and J.-S. Li.  Real-time inference and detection of disruptive eeg networks for epileptic seizures.  Sci Rep, 10:8653, 2020.
  * [7] E. Brown, J. Moehlis, and P. Holmes.  On the phase reduction and response dynamics of neural oscillator populations.  Neural Comput, 16(4):673–715, 2004.
  * [8] M. Budišić, R. Mohr, and I. Mezić.  Applied Koopmanism.  Chaos, 22:047510, 2012.
  * [9] N. Črnjarić-Žic, S. Maćešić, and I. Mezić.  Koopman operator spectrum for random dynamical systems.  J Nonlinear Sci, 30:2007–2056, 2020.
  * [10] S. Das and D. Giannakis.  On harmonic hilbert spaces on compact abelian groups.  arXiv [Preprint] https://arxiv.org/abs/1912.11664, accessed 22 April 2022.
  * [11] B. Dong, B. Yu, and Y. Yu.  A homotopy method for finding all solutions of a multiparameter eigenvalue problem.  SIAM J Matrix Anal Appl, 37(2):550–571, 2016.
  * [12] K. Fujii and Y. Kawahara.  Dynamic mode decomposition in vector-valued reproducing kernel Hilbert spaces for extracting dynamical structure among observables.  Neural Netw, 117:94–103, 2019.
  * [13] K. Fujii, N. Takeishi, B. Kibushi, M. Kouzaki, and Y. Kawahara.  Data-driven spectral analysis for coordinative structures in periodic human locomotion.  Sci Rep, 9:16755, 2019.
  * [14] C. W. J. Granger.  Investigating causal relations by econometric models and cross-spectral methods.  Econometrica, 37(3):424–438, 1969.
  * [15] Y. Hashimoto, I. Ishikawa, M. Ikeda, F. Komura, T. Katsura, and Y. Kawahara.  Analysis via orthonormal systems in reproducing kernel Hilbert $C^{*}$-modules and applications.  arXiv:2003.00738, 2020.
  * [16] Y. Hashimoto, I. Ishikawa, M. Ikeda, Y. Matsuo, and Y. Kawahara.  Krylov subspace method for nonlinear dynamical systems with random noise.  JMLR, 21(172):1–29, 2020.
  * [17] B. Heersink, M. A. Warren, and H. Hoffmann.  Dynamic mode decomposition for interconnected control systems.  arXiv [Preprint] https://arxiv.org/abs/1709.02883, accessed 19 June 2022.
  * [18] J. Hu and Y. Lan.  Koopman analysis in oscillator synchronization.  arXiv:2009.09227, 2020.
  * [19] I. Ishikawa, K. Fujii, M. Ikeda, Y. Hashimoto, and Y. Kawahara.  Metric on nonlinear dynamical systems with Perron-Frobenius operators.  In NeurIPS, 2018.
  * [20] H. Kadri, E. Duflos, P. Preux, S. Canu, A. Rakotomamonjy, and J. Audiffren.  Operator-valued kernels for learning from functional response data.  Journal of Machine Learning Research, 17(20):1–54, 2016.
  * [21] Y. Kawahara.  Dynamic mode decomposition with reproducing kernels for Koopman spectral analysis.  In NeurIPS, 2016.
  * [22] V. Khlobystov, B. Podlevskii, and O. Yaroshko.  One method of solving linear multiparameter eigenvalue problem.  IJSBAR, 7(1):53–65, 2015.
  * [23] S. Klus, F. Nüske, S. Peitz, J.-H. Niemann, C. Clementi, and C. Schütte.  Data-driven approximation of the Koopman generator: Model reduction, system identification, and control.  Phys D, 406:132416, 2020.
  * [24] B. O. Koopman.  Hamiltonian systems and transformation in Hilbert space.  Proc Natl Acad Sci USA, 17(5):315–318, 1931.
  * [25] B. Kralemann, M. Frühwirth, A. Pikovsky, M. Rosenblum, T. Kenner, J. Schaefer, and M. Moser.  In vivo cardiac phase response curve elucidates human respiratory heart rate variability.  Nat Commun, 4:2418, 2013.
  * [26] Y. Kuramoto.  Chemical Oscillations, Waves, and Turbulence.  Springer, 1984.
  * [27] J. N. Kutz.  Data-Driven Modeling & Scientific Computation: Methods for Complex Systems & Big Data.  OUP Oxford, 2013.
  * [28] A. Mauroy and I. Mezić.  On the use of fourier averages to compute the global isochrons of (quasi)periodic dynamics.  Chaos, 22(3):033112, 2012.
  * [29] A. Mauroy and I. Mezić.  Global computation of phase-amplitude reduction for limit-cycle dynamics.  Chaos, 28(7):073108, 2018.
  * [30] A. Mauroy, I. MeziÄ, and J. Moehlis.  Isostables, isochrons, and Koopman spectrum for the actionâangle representation of stable fixed point dynamics.  Phys D, 261:19–30, 2013.
  * [31] Z. Mei and T. Oguchi.  Network structure identification via koopman analysis and sparse identification.  NOLTA, IEICE, 13(2):477–492, 2022.
  * [32] I. Mezić.  Spectral properties of dynamical systems, model reduction and decompositions.  Nonlinear Dyn., 41:309–325, 2005.
  * [33] H. Q. Minh, L. Bazzani, and V. Murino.  A unifying framework in vector-valued reproducing kernel Hilbert spaces for manifold regularization and co-regularized multi-view learning.  Journal of Machine Learning Research, 17(25):1–72, 2016.
  * [34] F. Mori and H. Kori.  Noninvasive inference methods for interaction and noise intensities of coupled oscillators using only spike time data.  Proc Natl Acad Sci USA, 119(6):e2113620119, 2022.
  * [35] H. Nakao.  Phase reduction approach to synchronisation of nonlinear oscillators.  Contemporary Physics, 57(2):188–214, 2016.
  * [36] H. Nakao, T. Yanagita, and Y. Kawamura.  Phase-reduction approach to synchronization of spatiotemporal rhythms in reaction-diffusion systems.  Phys Rev X, 4:021032, 2014.
  * [37] H. Nakao, S. Yasui, M. Ota, K. Arai, and Y. Kawamura.  Phase reduction and synchronization of a network of coupled dynamical elements exhibiting collective oscillations.  Chaos, 28(4):045103, 2018.
  * [38] M. J. Paul, P. Indic, and W. J. Schwartz.  Social synchronization of circadian rhythmicity in female mice depends on the number of cohabiting animals.  Biol Lett, 11(6):20150204, 2015.
  * [39] B. Pietras and A. Daffertshofer.  Network dynamics of coupled oscillators and phase reduction techniques.  Phys Rep, 819:1–105, 2019.
  * [40] C. W. Rowley, I. Mezić, S. Bagheri, P. Schlatter, and D. S. Henningson.  Spectral analysis of nonlinear flows.  J Fluid Mech, 641:115â–127, 2009.
  * [41] B. Schölkopf and A. J. Smola.  Learning with Kernels: Support Vector Machines, Regularization, Optimization, and Beyond.  MIT Press, 2001.
  * [42] S. Shirasaka, W. Kurebayashi, and H. Nakao.  Phase-amplitude reduction of transient dynamics far from attractors for limit-cycling systems.  Chaos, 27(2):023119, 2017.
  * [43] T. Stankovski, V. Ticcinelli, P. V. E. McClintock, and A. Stefanovska.  Coupling functions in networks of oscillators.  New J Phys, 17:035002, 2015.
  * [44] N. Takeishi, Y. Kawahara, and T. Yairi.  Subspace dynamic mode decomposition for stochastic Koopman analysis.  Phy Rev E, 96:033310, 2017.
  * [45] I. Tokuda, Z. Levnajic, and K. Ishimura.  A practical method for estimating coupling functions in complex dynamical systems.  Philos Trans A Math Phys Eng Sci, 377(2160):20190015, 2019.
  * [46] S. Wang, E. D. Herzog, I. Z. Kiss, W. J. Schwartz, G. Bloch, M. Sebek, D. Granados-Fuentes, L. Wang, and J.-S. Li.  Inferring dynamic topology for decoding spatiotemporal structures in complex heterogeneous networks.  Proc Natl Acad Sci USA, 115(37):9300–9305, 2018.
  * [47] N. Watanabe, Y. Kato, S. Shirasaka, and H. Nakao.  Optimization of linear and nonlinear interaction schemes for stable synchronization of weakly coupled limit-cycle oscillators.  Phys Rev E, 100:042205, 2019.
  * [48] M. Williams, I. Kevrekidis, and C. Rowley.  A dataâdriven approximation of the Koopman operator: extending dynamic mode decomposition.  J Nonlinear Sci, 25:1307–â1346, 2015.
  * [49] D. Wilson.  A data-driven phase and isostable reduced modeling framework for oscillatory dynamical systems.  Chaos, 30(1):013121, 2020.
  * [50] T. Young, M. Palta, J. Dempsey, P. Peppard, F. Nieto, and K. Hla.  Burden of sleep apnea: rationale, design, and major findings of the Wisconsin sleep cohort study.  WMJ, 108(5):246–9, 2009.
  * [51] G. Q. Zhang, L. Cui, R. Mueller, S. Tao, M. Kim, M. Rueschman, S. Mariani, D. Mobley, and S. Redline.  The national sleep research resource: towards a sleep data commons.  J Am Med Inf Assoc, 25(10):1351–1358, 2018.

## A. Multiparameter eigenvalue problem

The multiparameter eigenvalue problem (MEP) is a generalization of the
standard eigenvalue problem and was developed originally for solving ordinary
differential equations with multiple parameters [3, 2]. Several algorithms to
solve this problem have also been investigated in numerical analysis [22, 11,
4]. In this paper, we further generalize the multiparameter eigenvalue problem
and get GMEP to derive the phase model.

## B. Computation of the solution of the problem (4)

### B.1. Solving the problem (4)

We explain the practical computation of finding $\lambda_{j}$,
$\\{a^{j}_{i,k}\\}_{i,k=1}^{N}$, and $\mathbf{u}^{j}$ for $j=1,\ldots,M$ in
Eq. (4). We first focus on the case where the PCF $\Gamma_{i}$ in the phase
model (2) equals the zero function for any $i=1,\ldots,N$, that is, there are
no interactions among dynamical elements. In this case, the relationship
between eigenvalues and eigenfunctions of the Koopman operator and the phase
model is directly deduced from the case of $N=1$ as follows. Let
$\lambda\in\mathbb{C}$ be an eigenvalue of the Koopman operator $K$ and let
$\mathbf{u}=[u_{1},\ldots,u_{N}]\in\mathcal{H}$ be a corresponding
eigenvector, that is,

$K\mathbf{u}=\lambda\mathbf{u}.$ (10)

By the definition of $K$, for $i=1,\ldots,N$, we have

$u_{i}(X_{i}(t+\Delta t))=\lambda u_{i}(X_{i}(t)).$ (11)

Let $\theta_{i}(t)=\operatorname{arg}(u_{i}(X_{i}(t)))$. Note that the map
$\operatorname{arg}(u_{i}(\cdot))$ transforms $X_{i}(t)$ into the
corresponding phase $\theta_{i}(t)$. Assume $|\lambda|=1$, that is,
$\lambda=e^{\sqrt{-1}\omega\Delta t}$ for some $\omega\in[0,2\pi)$. Then, by
Eq. (11), the identity $|u_{i}(X_{i}(t+\Delta t))|=|u_{i}(X_{i}(t))|$ holds.
Thus, we have

$\theta_{i}(t+\Delta t)=\omega\Delta t+\theta_{i}(t),$ (12)

which means that each element $i$ oscillates with frequency $\omega$, and
there are no interactions among the elements.

Now, to generalize the above discussion to the case where the PCF $\Gamma_{i}$
is not necessarily zero for some $i=1,\ldots,N$, we consider the following
multiparameter eigenvalue problem and its generalization (instead of the
standard eigenvalue problem (11)). The multiparameter eigenvalue problem (MEP)
[2, 3] is to find a value $\lambda\in\mathbb{C}$, a sequence
$\\{a_{i,k}\\}_{i,k=1}^{N}$ in $\mathbb{C}$, and a vector
$\mathbf{u}=[u_{1},\ldots,u_{N}]\in\mathcal{H}$ such that

$\lambda^{-1}K\mathbf{u}=\sum_{i,k=1}^{N}a_{i,k}B_{i,k}\mathbf{u}.$ (13)

Here, $B_{i,k}$ is a linear operator on $\mathcal{H}$ that maps
$[u_{1},\ldots,u_{k}]\in\mathcal{H}$ to the vector whose $i$th element is
$u_{k}$ and other elements are $0$. As we will explain, Eq. (13) corresponds
to the case where the PCF $\Gamma_{i}(\psi_{1},\ldots,\psi_{N})$ is
represented only with the functions
$e^{-\sqrt{-1}\psi_{1}},\ldots,e^{-\sqrt{-1}\psi_{N}}$ (see Eq. (7)). Since
the PCF is more complicated, the MEP (13) does not always have a solution.
Therefore, we consider the generalized problem (4), which we call GMEP. We
first solve the following minimization problem:

$\min_{|\lambda|=1,\mathbf{u}\in\mathcal{H},a_{i,k}\in\mathbb{C}}\bigg{\|}\lambda^{-1}K\mathbf{u}-\sum_{i,k=1}^{N}a_{i,k}B_{i,k}\mathbf{u}\bigg{\|},$
(14)

where $\mathbf{u}$ is normalized to be
$\mathbf{u}(\mathbf{X}(s))\approx\mathbf{1}$ for some time $s$. Here,
$\mathbf{1}\in\mathbb{C}^{N}$ is the vector whose elements are all $1$ and
$\|\cdot\|$ is the norm induced by the inner product in the Hilbert space
$\mathcal{H}$. See Subsection B.2. Normalization of eigenfunctions for the
details of the normalization.

We consider compensating for the error of the minimization problem (14). In
fact, by the definition of the operator $B_{i,k}$ for $i,k=1,\ldots,N$, we
have

$\sum_{k=1}^{N}a_{i,k}B_{i,k}\mathbf{u}(\mathbf{X}(t))=\sum_{k=1}^{N}a_{i,k}u_{k}(X_{k}(t)),$

so that Eq. (13) means the effect of element $k$ on element $i$ is linear with
respect to $u_{k}(X_{k}(t))$. We consider compensating for the error of the
problem (14) by considering the nonlinear effect of $u_{k}(X_{k}(t))$. For
this purpose, we use $M$ eigenfunctions
$\\{\mathbf{u}^{j}=[u_{1}^{j},\ldots,u_{N}^{j}]\\}_{j=1}^{M}$ of the standard
eigenvalue problem (11) whose $k$th elements $(u^{j}_{k}(X_{k}(\cdot))$ are
approximately equal to $u_{k}(X_{k}(\cdot))^{j}$ for $j=1,\ldots,M$. Let
$\omega\in[0,2\pi)$ satisfy $\lambda=e^{\sqrt{-1}\omega\Delta t}$, where
$\lambda$ is a solution of the problem (14). We assume that there exist $M$
sets of $N$ standard eigenvalues $\\{\lambda_{i}^{j}\\}_{i=1}^{N}$
($j=1,\ldots,M$) of $K$ which are represented as $\lambda_{i}^{j}\approx
e^{\sqrt{-1}j\omega\Delta t}$. See Subsection B.3. Assumption about the
eigenvalues of $K$ for the details and the validity of this assumption. Let
$\mathbf{u}^{1}=\mathbf{u}$ and $a^{1}_{i,k}=a_{i,k}$, which is the solution
of the problem (14). In addition, for $j=2,\ldots,M$, we set $\mathbf{u}^{j}$
so that it satisfies $K\mathbf{u}^{j}\approx e^{-\sqrt{-1}j\omega\Delta
t}\mathbf{u}^{j}$ and is normalized to be
$\mathbf{u}^{j}(\mathbf{X}(s))=\mathbf{1}$. See Subsection B.2. Normalization
of eigenfunctions for the details of the normalization. For $j=2,\ldots,M$, we
compute recursively $\\{a^{j}_{i,k}\\}_{i,k=1}^{N}$, which is the solution of
the minimization problem

$\displaystyle\min_{a_{i,k}^{j}\in\mathbb{C}}\bigg{\|}e^{-\sqrt{-1}\omega\Delta
t}K\mathbf{u}^{1}\odot\cdots\odot e^{-\sqrt{-1}j\omega\Delta
t}K\mathbf{u}^{j}$
$\displaystyle-\bigg{(}\sum_{i,k=1}^{N}a^{1}_{i,k}B_{i,k}\mathbf{u}^{1}\bigg{)}\odot\cdots\odot\bigg{(}\sum_{i,k=1}^{N}a^{j}_{i,k}B_{i,k}\mathbf{u}^{j}\bigg{)}\bigg{\|}.$
(15)

Note that we need additional structure of product in $\mathcal{H}$ for
satisfying $\mathbf{u}\odot\mathbf{v}\in\mathcal{H}$ for
$\mathbf{u},\mathbf{v}\in\mathcal{H}$. RKHSs that have the product structure
are discussed in [10], and we can choose one of these RKHSs as $\mathcal{H}$.
However, the norm in Eq. (15) does not need to be that in $\mathcal{H}$.
Therefore, another option is to embed $\mathbf{u}$ and $\mathbf{v}$ to a
Hilbert space that contains $\mathbf{u}\odot\mathbf{v}$ such as the tensor
product $\mathcal{H}\otimes\mathcal{H}$ of an RKHS $\mathcal{H}$. If the
objective function of the minimization problem (15) becomes $0$ with solutions
$\\{a_{i,k}^{j}\\}_{i,k=1}^{N}$ $(j=2,\ldots,M)$, then we obtain a solution of
Eq. (4).

### B.2. Normalization of eigenfunctions

Since we focus on multiple dynamical elements, we have to apply a
normalization for collecting the phases and equalizing the magnitudes of
effects of elements on the model. We first consider the case where the PCF
$\Gamma_{i}$ is the zero function for any $i=1,\ldots,N$. Assume that there
exists a time $s$ that satisfies $u_{i}(X_{i}(s))=1$ for any $i=1,\ldots,N$,
that is, all the phases $\theta_{i}$ coincide at $s$ as
$\theta_{i}(s)=\operatorname{arg}(u_{i}(X_{i}(s)))=0$. In fact, if there are
$N$ distinct eigenfunctions
$\mathbf{v}_{1},\ldots,\mathbf{v}_{N}\in\mathcal{H}$ of an eigenvalue
$\lambda$ and if
$\mathbf{v}_{1}(\mathbf{X}(s)),\ldots,\mathbf{v}_{N}(\mathbf{X}(s))\in\mathbb{C}^{N}$
are linearly independent, then there exists a unique
$(c_{1},\ldots,c_{N})\in\mathbb{C}^{N}$ such that
$\sum_{i=1}^{N}c_{i}\mathbf{v}_{i}(\mathbf{X}(s))=\mathbf{1}$ and
$\sum_{i=1}^{N}c_{i}\mathbf{v}_{i}$ is an eigenfunction of $K$ with respect to
the eigenvalue $\lambda$. In this case, if we set
$\mathbf{u}=\sum_{i=1}^{N}c_{i}\mathbf{v}_{i}$ in Eq. (11), by Eq. (12),
$\theta_{1}(t)=\cdots=\theta_{N}(t)$ hold for any $t$, that is, each element
$i$ oscillates with frequency $\omega$, and there are no interactions among
elements.

Next, we consider the case where the PCF $\Gamma_{i}$ is not necessarily zero
for some $i=1,\ldots,N$. Assume that there exists a time $s$ that satisfies
$u_{i}^{j}(X_{i}(s))\approx 1$ for any $i=1,\ldots,N$ and $j=1,\ldots,M$, that
is, all the phases $\theta_{i}^{j}$ approximately coincide at $s$ as
$\theta_{i}^{j}(s)=\operatorname{arg}(u_{i}^{j}(X_{i}(s)))\approx 0$ and all
the magnitudes $|u_{i}^{j}(X_{i}(s))|$ are nearly equal to $1$. In fact, for
$j=1,\ldots,M$, if there are $N$ distinct eigenfunctions
$\mathbf{v}_{1}^{j},\ldots,\mathbf{v}_{N}^{j}\in\mathcal{H}$ of the standard
eigenvalues each of which is nearly equal to $e^{-\sqrt{-1}j\omega\Delta t}$
and if
$\mathbf{v}_{1}^{j}(\mathbf{X}(s)),\ldots,\mathbf{v}_{N}^{j}(\mathbf{X}(s))\in\mathbb{C}^{N}$
are linearly independent, then there exists unique
$(c_{1}^{j},\ldots,c_{N}^{j})\in\mathbb{C}^{N}$ such that
$\sum_{i=1}^{N}c_{i}^{j}\mathbf{v}_{i}^{j}(\mathbf{X}(s))=\mathbf{1}$ and
$K\sum_{i=1}^{N}c_{i}^{j}\mathbf{v}_{i}^{j}\approx e^{-\sqrt{-1}j\omega\Delta
t}\sum_{i=1}^{N}c_{i}^{j}\mathbf{v}_{i}^{j}$ hold. Thus, for the problem (14),
by setting $\lambda=e^{-\sqrt{-1}\omega\Delta t}$ and
$\mathbf{u}=\sum_{i=1}^{N}c_{i}^{1}\mathbf{v}_{i}^{1}(\mathbf{X}(s))$ as
initial values and applying an iterative method such as gradient descent
method, we can obtain a solution of the problem (14) satisfying
$\mathbf{u}(\mathbf{X}(s))\approx\mathbf{1}$. For the problem (15), we set
$\mathbf{u}^{j}=\sum_{i=1}^{N}c_{i}^{j}\mathbf{v}_{i}^{j}$ for $j=2,\ldots,M$.
Then, we have $\mathbf{u}^{j}(\mathbf{X}(s))=\mathbf{1}$.

### B.3. Assumption about the eigenvalues of $K$

If all the $N$ elements oscillate and synchronize completely, we expect that
there exist $\omega\in[0,2\pi)$ and eigenvalues of $K$ equal to
$e^{\sqrt{-1}j\omega\Delta t}$ for $j=1,\ldots,M$ and whose geometric
multiplicities are $N$. This is because if the phase variable $\theta_{i}$
increases with frequency $\omega$, then the variable $j\theta_{i}$ increases
with frequency $j\omega$. On the other hand, if the elements do not completely
synchronize and interact weakly, such eigenvalues do not always exist, but we
expect there are $N$ eigenvalues that are nearly equal to
$e^{\sqrt{-1}j\omega\Delta t}$. Therefore, we assume that there exist $M$ sets
of $N$ standard eigenvalues $\\{\lambda_{i}^{j}\\}_{i=1}^{N}$ ($j=1,\ldots,M$)
that are represented as $\lambda_{i}^{j}\approx e^{\sqrt{-1}j\omega\Delta t}$.

## C. Derivation of Eq. (6)

We explain the derivation of the transformation of Eq. (5) into Eq. (6). Let
$\theta_{i}^{j}(t)=\operatorname{arg}({u}_{i}^{j}(X_{i}(t)))$ and let
$r_{i}^{j}(t)=|{u}_{i}^{j}(X_{i}(t))|$. Then, Eq. (5) implies

$\displaystyle\bigg{(}\prod_{j=1}^{M}r_{i}^{j}(t+\Delta
t)\bigg{)}e^{\sqrt{-1}\sum_{j=1}^{M}\theta_{i}^{j}(t+\Delta t)}$
$\displaystyle=e^{\sqrt{-1}\sum_{j=1}^{M}(j\omega\Delta
t+\theta_{i}^{j}(t))}\prod_{j=1}^{M}\sum_{k=1}^{N}{a}_{i,k}^{j}r_{k}^{j}(t)e^{\sqrt{-1}(\theta_{k}^{j}(t)-\theta_{i}^{j}(t))}.$
(16)

Note that by the normalization, we have $r_{i}^{j}(t)\approx 1$. Therefore, by
ignoring the factor $r_{i}^{j}(t)$ and calculating arguments of both sides of
Eq. (16), we obtain

$\displaystyle\sum_{j=1}^{M}\theta_{i}^{j}(t+\Delta t)$
$\displaystyle\approx\sum_{j=1}^{M}\bigg{(}j\omega\Delta
t+\theta_{i}^{j}(t)+\operatorname{arg}\bigg{(}\sum_{k=1}^{N}{a}_{i,k}^{j}e^{\sqrt{-1}(\theta_{k}^{j}(t)-\theta_{i}^{j}(t))}\bigg{)}\bigg{)}.$
(17)

Since $\lambda_{j}$ and $\mathbf{u}^{j}$ for $j=2,\ldots,M$ satisfy Eq. (15)
and $a_{i,k}^{j}\approx 0$, we have $\theta_{i}^{j}(t+\Delta t)\approx
j\omega\Delta t+\theta_{i}^{j}(t)$. Thus, we set $\theta_{i}=\theta_{i}^{1}$
and regard $\theta_{i}^{j}\approx j\theta_{i}$. Therefore, we obtain Eq. (6).

## D. Vector-valued Reproducing Kernel Hilbert Spaces (vvRKHSs)

We review the theory of vvRKHSs [20, 33].

To construct a vvRKHS, we begin by a positive definite kernel. Let
$\mathcal{X}=\mathbb{R}^{d}$ or $[0,2\pi)$. A map
$\Phi:\mathcal{X}^{N}\times\mathcal{X}^{N}\to\mathbb{C}^{N\times N}$ is called
a $\mathbb{C}^{N\times N}$-valued positive definite kernel if it satisfies the
following conditions:

  1. 1.

$\Phi(\mathbf{x},\mathbf{y})=\Phi(\mathbf{y},\mathbf{x})^{*}$ for
$\mathbf{x},\mathbf{y}\in\mathcal{X}^{N}$,

  2. 2.

$\sum_{i,j=1}^{n}\mathbf{d}_{i}^{*}\Phi(\mathbf{x}_{i},\mathbf{x}_{j})\mathbf{d}_{j}\geq
0$ for $n\in\mathbb{N}$, $\mathbf{d}_{i}\in\mathbb{C}^{N}$,
$\mathbf{x}_{i}\in\mathcal{X}^{N}$.

Here, for a matrix $A$, $A^{*}$ represents the Hermitian conjugate of $A$. For
$\mathbf{x}\in\mathcal{X}^{N}$, let $\phi_{\mathbf{x}}$ be a
$\mathbb{C}^{N\times N}$-valued map on $\mathcal{X}^{N}$ defined as
$\phi_{\mathbf{x}}=\Phi(\cdot,\mathbf{x})$. The following vector-valued
function space is constructed:

$\mathcal{H}_{0}:=\bigg{\\{}\sum_{i=1}^{n}\phi_{\mathbf{x}_{i}}\mathbf{d}_{i}\bigg{|}\
n\in\mathbb{N},\ \mathbf{d}_{i}\in\mathbb{C}^{N},\
\mathbf{x}_{i}\in\mathcal{X}^{N}\bigg{\\}}.$

Then, we define a map
$\left\langle\cdot,\cdot\right\rangle:\mathcal{H}_{0}\times\mathcal{H}_{0}\to\mathbb{C}$
as

$\displaystyle\bigg{\langle}\sum_{i=1}^{n}\phi_{\mathbf{x}_{i}}\mathbf{d}_{i},\sum_{j=1}^{l}\phi_{\mathbf{y}_{j}}\mathbf{h}_{j}\bigg{\rangle}:=\sum_{i=1}^{n}\sum_{j=1}^{l}\mathbf{d}_{i}^{*}\Phi(\mathbf{x}_{i},\mathbf{y}_{j})\mathbf{h}_{j}.$

By the above two properties of $\Phi$, $\left\langle\cdot,\cdot\right\rangle$
is well-defined, satisfies the axiom of inner products, and has the
reproducing property, that is,

$\left\langle\phi_{\mathbf{x}}\mathbf{d},\mathbf{u}\right\rangle=\mathbf{d}^{*}\mathbf{u}(\mathbf{x}),$

for $\mathbf{u}\in\mathcal{H}_{0}$, $\mathbf{x}\in\mathcal{X}^{N}$, and
$\mathbf{d}\in\mathbb{C}^{N}$. The completion of $\mathcal{H}_{0}$ is called
the vector-valued reproducing kernel Hilbert space (vvRKHS) associated with
$\Phi$ and denoted by $\mathcal{H}$.

By setting a proper positive definite kernel, we can construct the required
vector-valued function space, that is, the vector-valued function $\mathbf{v}$
having the form $\mathbf{v}(\mathbf{x})=[{v}_{1}(x_{1}),\ldots,v_{N}(x_{N})]$
for $\mathbf{x}=[x_{1},\ldots,x_{N}]\in\mathcal{X}^{N}$, where $v_{i}$ for
$i=1,\ldots,N$ is a complex-valued function on $\mathcal{X}$. More precisely,
let $\Phi$ be a $\mathbb{C}^{N\times N}$-valued function defined as
$[\Phi(\mathbf{x}_{1},\mathbf{x}_{2})]_{i,j}={k}((x_{1,i},i),(x_{2,j},j))$,
where ${k}$ is a complex-valued positive definite kernel (that is, $k$ is a
$\mathbb{C}^{1\times 1}$-valued positive definite kernel). Typical examples of
$k$ are the Gaussian kernel on $\mathbb{R}^{d}$ defined as
$k(x,y)=e^{-c\|x-y\|_{2}}$ for some $c>0$ and Laplacian kernel on
$\mathbb{R}^{d}$ defined as $k(x,y)=e^{-c\|x-y\|_{1}}$ for some $c>0$. Here,
$\|x\|_{p}$ for $x=[x_{1},\ldots,x_{d}]\in\mathbb{R}^{d}$ and $p=1,2$ is
defined as $\|x\|_{p}=(\sum_{i=1}^{d}|x_{i}|^{p})^{1/p}$. By the definition of
$\Phi$, for a function $\mathbf{v}$ in $\mathcal{H}$, $\mathbf{v}(\mathbf{x})$
has the form $[{v}_{1}(x_{1}),\ldots,v_{N}(x_{N})]$, and we can define the
Koopman operator $K$ as Eq. (3).

## E. Estimation of Koopman operators on vvRKHSs

To obtain the frequency $\omega$, transformation of $X_{i}(t)$ into
$\theta_{i}(t)$, and PCF $\Gamma_{i}$ in Eq. (5) only with given data, we need
to estimate the Koopman operator $K$. Here, we show an approach to estimate
the Koopman operator. Let
$[x_{1,0},\ldots,x_{N,0}]=\mathbf{x}_{0},\ldots,[x_{1,T},\ldots,x_{N,T}]=\mathbf{x}_{T}\in\mathcal{X}^{N}$
be given data generated by the dynamical system (1) at time
$t_{0},\ldots,t_{T}$. We assume that $S$ sequences
$\\{\mathbf{x}^{1}_{0},\ldots,\mathbf{x}^{1}_{T}\\},\ldots,\\{\mathbf{x}^{S}_{0},\ldots,\mathbf{x}^{S}_{T}\\}$
of data are given. For the Hilbert space $\mathcal{H}$, we use a vvRKHS.

We now generate a subspace of $\mathcal{H}$ by the given data and estimate the
Koopman operator $K$ on the subspace. For $s=0,\ldots,T$ and $j=1,\ldots,N$,
let
$\boldsymbol{\eta}_{s,i}=1/S\sum_{j=1}^{S}\phi_{\mathbf{x}^{j}_{s}}\mathbf{e}_{i}$
and let $\mathcal{Q}$ be the subspace of $\mathcal{H}$ spanned by
$\\{\boldsymbol{\eta}_{s,i}\,\mid\,s=0,\ldots,T-1,\ i=1,\ldots,N\\}$.
Moreover, to reduce noise and extract crucial information from observed data,
we use principal component analysis (PCA).

We use the kernel PCA to obtain a subspace of $\mathcal{Q}_{1}$ [41]. Let
$\mathbf{G}\in\mathbb{C}^{NT\times NT}$ be the Gram matrix whose
$(Ns+i,Nt+j)$-element is defined as
$\langle\boldsymbol{\eta}_{s,i},\boldsymbol{\eta}_{t,j}\rangle$. Note that
$\mathbf{G}$ is a Hermitian positive definite matrix. Let
$\lambda_{1},\ldots\lambda_{T^{\prime}}>0$ be the largest $T^{\prime}$
eigenvalues of $\mathbf{G}$ and let
$\mathbf{v}_{1},\ldots,\mathbf{v}_{T^{\prime}}\in\mathbb{C}^{NT}$ be the
corresponding orthonormal eigenvectors. Then, $P=QQ^{*}$, where
$Q=H\mathbf{V}$,
$H=[\boldsymbol{\eta}_{0,1},\ldots,\boldsymbol{\eta}_{0,N},\ldots,\boldsymbol{\eta}_{T-1,1},\ldots,\boldsymbol{\eta}_{T-1,N}]$,
and
$\mathbf{V}=[(\lambda_{1})^{-1/2}\mathbf{v}_{1},\ldots,(\lambda_{T^{\prime}})^{-1/2}\mathbf{v}_{T^{\prime}}]$.
The estimation of the Koopman operator $K$ is obtained by projecting a vector
in $\mathcal{H}$ onto $\mathcal{Q}_{1}$, acting $K$ on the projected vector,
and then projecting it back to the original space $\mathcal{H}$. The estimated
operator is represented as $PKP$, and $Q^{*}KQ$ is a matrix representation of
the estimated operator. We denote by $\tilde{\mathbf{K}}$ the matrix
representation of the estimated operator. Let
$\tilde{\mathbf{G}}\in\mathbb{C}^{NT\times NT}$ be the matrix whose
$(Ns+i,Nt+j)$-element is defined as
$\langle\boldsymbol{\eta}_{s+1,i},\boldsymbol{\eta}_{t,j}\rangle$. The matrix
$\tilde{\mathbf{K}}$ is shown to be represented by using $\tilde{\mathbf{G}}$
and $\mathbf{V}$. By the identities
$\langle\boldsymbol{\eta}_{s+1,i},\boldsymbol{\eta}_{t,j}\rangle=\langle\boldsymbol{\eta}_{s,i},K\boldsymbol{\eta}_{t,j}\rangle$
and $Q=H\mathbf{V}$, we have

$\displaystyle\tilde{\mathbf{K}}=\mathbf{V}^{*}H^{*}KH\mathbf{V}=\mathbf{V}^{*}\tilde{\mathbf{G}}\mathbf{V}.$

Estimations of eigenvalues of $K$ are obtained by computing eigenvalues of
$\tilde{\mathbf{K}}$. In addition, estimations of vector-valued eigenfunctions
of $K$ are computed as $Q\mathbf{u}$, where $\mathbf{u}$ is an eigenvector of
$\tilde{\mathbf{K}}$. For all the numerical experiments in this paper, we set
$[\Phi(\mathbf{x}_{1},\mathbf{x}_{2})]_{i,j}=e^{-1/d\|x_{1,i}-{x}_{2,j}\|_{1}}e^{-0.1|i-j|}$
for $\mathbf{x}_{1},\mathbf{x}_{2}\in\mathbb{R}^{d}$. In addition, we set
$T^{\prime}=3/4T$ for all the numerical experiments.

## F. Estimation of the strength of interactions

By computing the coefficients $a_{i,k}^{j}$ in Eq. (6) using the estimated
Koopman operator, we can estimate the strength of interactions between
dynamical elements. Here, we show additional experimental results regarding
the estimation of the strength of interactions.

We consider real-world time-series data describing the body temperature of
cohabiting mice [38]. Five mice were initially housed in a 12 hours light and
12 hours dark (LD) cycle. Two mice (indices 1 and 2) were housed in one of two
opposite LD cycles (L:D), and the other three mice (indices 3, 4, and 5) were
housed in the other LD cycle (D:L). Then, they began cohabiting. According to
the results by Paul [38, Figure 1 (b)], the cohabitation caused the
synchronization of their circadian rhythms. For $i=1,\ldots,5$, we set $y_{i}$
as a variable describing the body temperature of mouse $i$ and set
$z_{i}(t)=y_{i}(t+1)$ and $X_{i}(t)=[y_{i}(t),z_{i}(t)]\in\mathbb{R}^{2}$.
Thus, $N$ and $\mathcal{X}$ are $5$ and $\mathbb{R}^{2}$ in this example. We
estimated the Koopman operator with the data. We first computed standard
eigenvalues of the Koopman operator, which are illustrated in Fig. 5 A. Then,
we estimated the strength of interactions by obtaining the magnitude of
$a_{i,k}=a_{i,k}^{j}$ $(i,k=1,\ldots,5)$ by solving the MEP using the
eigenvalues colored in red in Fig. 5 A and the corresponding eigenfunctions.
The result is shown in Fig. 5 B. We see that the magnitudes of $a_{i,k}$
between mice initially housed in the same LD cycle tends to be large.

Note that we can also estimate the strength of interactions by using the
Granger causality test [14]. Although the Granger causality test only focuses
on interactions, using our framework enables us to estimate the common
frequency $\omega_{0}$ and the strength of interactions simultaneously. We
also remark that for $N\geq 3$, estimating the rigorous form of the PCF is
challenging. Thus, here, we focused on the strength of interactions.

A |  B
---|---
Figure 5: The real-world time-series data describing the body temperature of
cohabiting mice. (A) The standard eigenvalues of the estimated Koopman
operator. (B) The heat map of the magnitude of the sum of interaction
coefficient $|a_{i,k}^{1}|+|a_{k,i}^{1}|$ obtained by the estimated Koopman
operator.

## G. Experimental details

### G.1. Stuart–Landau model with two oscillators

Regarding the generation of data, we generated
$\mathbf{x}_{0}^{i},\mathbf{x}_{1}^{i},\ldots,\mathbf{x}_{2000}^{i}\in(\mathbb{R}^{2})^{2}$
for $i=1,\ldots,S$, where $S=10$ from Eq. (8) with time interval $\Delta
t=0.25$. For generating Fig. 2 C, we set the initial values as
$\mathbf{x}^{i}_{0}=[[0.1\sin(-3/4\pi+3/2\pi\cdot
i/S),0.1\cos(-3/4\pi+3/2\pi\cdot i/S)],[0,0.1]]$ for $i=1,\ldots,S$. For
generating Figs. 2 D–G, we set the initial value as
$x^{i}_{0}=[[0.9\sin(-3/4\pi+6/4\pi\cdot i/S),0.9\cos(-3/4\pi+6/4\pi\cdot
i/S)],[0,0.9]]$.

To obtain Eq. (6), we solved the multiparameter eigenvalue problem according
to Section B. Computation of the solution of the problem (4). There are pairs
of eigenvalues $\lambda^{j}_{1},\lambda^{j}_{2}\in\mathbb{C}$ which satisfy
$|\lambda_{i}^{j}|\approx 1$ and $\lambda_{i}^{j}\approx
e^{\sqrt{-1}j\omega_{0}\Delta t}$ for $\omega_{0}\approx 0.25$ and for $i=1,2$
and $j=1,2,\ldots$, which are colored in red in Fig. 2 A. Let
$\mathbf{v}^{j}_{1},\mathbf{v}^{j}_{2}\in\mathcal{H}$ be the eigenfunctions
with respect to $\lambda^{j}_{1},\lambda^{j}_{2}$, respectively. Since we
expect the orbits of the elements approach as time progresses, we can set
$t_{0}$ as a sufficiently large time, for example, $t_{0}=\Delta tl$ with
$l=2000$. In this example, we set $l=2000$. We computed
$c^{1}_{1},c^{1}_{2}\in\mathbb{C}$ such that
$\sum_{i=1}^{N}c^{1}_{i}v^{1}_{i}(x^{S/2}_{l})=\mathbf{1}$. We solved the
minimization problem (14) by setting $\lambda=\lambda_{1}$,
$\mathbf{u}=\sum_{i=1}^{N}c^{1}_{i}\mathbf{v}^{1}_{i}$, and $a_{i,k}=1$ for
$i=k$ and $a_{i,k}=0.0015$ for $i\neq k$ as initial values. We obtained
$\lambda=e^{\sqrt{-1}\omega\Delta t}$ with $\omega=0.2498$. Then, we set $M=4$
and computed a solution of the minimization problem (15) recursively for
$j=2,3,4$, where $\lambda_{j}=e^{\sqrt{-1}j\omega\Delta t}$, $\mathbf{u}_{1}$
and $a_{i,k}^{1}$ are the solutions of the minimization problem (14), and
$a_{i,k}^{l}$ is the solution of the minimization problem (15) for $l<j$.
Moreover, for $l=2000$, we computed $c_{1}^{j},c_{2}^{j}\in\mathbb{C}$ such
that $\sum_{i=1}^{N}c_{i}^{j}v_{i}^{j}(x^{S/2}_{l})=\mathbf{1}$ and set
$\mathbf{u}_{j}=\sum_{i=1}^{N}c_{i}\mathbf{v}_{i}^{j}$ for $j=2,\ldots M$.

It is noted that appropriate data for obtaining a good global estimation of
the phase function and PCF are generally different. For the phase function, if
we need to estimate it globally (not only the neighborhood of the limit
cycle), we need data which is widely distributed. On the other hand, for
estimating the PCF, the samples which are not in the neighborhood of the limit
cycle make the estimation not accurate.

### G.2. FitzHugh–Nagumo model with two oscillators

We generated
$\mathbf{x}^{i}_{0},\mathbf{x}^{i}_{1},\ldots,\mathbf{x}^{i}_{3000}\in(\mathbb{R}^{2})^{2}$
for $i=1,\ldots,20$ from Eq. (8) with time interval $\Delta t=0.1$, and
estimated the Koopman operator. For generating Fig. 3 C, we set the initial
values as $\mathbf{x}^{i}_{0}=[[0.02\sin(-\pi+2\pi\cdot
i/20),0.02\cos(-\pi+2\pi\cdot i/20)],[0,0.02]]$ for $i=1,\ldots,10$. For
generating Figs. 3 D–G, we set the initial value as
$\mathbf{x}^{i}_{0}=[[0.05+0.35\sin(-\pi+2\pi\cdot
i/20),0.05+0.12\sin(-\pi+2\pi\cdot i/20)],[0.05,0.17]]$ for $i=1,\ldots,20$
and removed the first 200 data from the generated sequence. We obtained Eq.
(6) in the same manner as the case of the Stuart–Landau model. We set
$a_{i,k}=1$ for $i=k$ as an initial value. In fact, the scale of the PCF
$\Gamma_{i}$ estimated by our method depends on the initial value of $a_{i,k}$
for $i\neq k$. Thus, we have to choose an appropriate initial value. In this
example, we set $a_{i,k}=s_{\operatorname{n}}/s_{\operatorname{d}}$, where
$s_{\operatorname{d}}=\sum_{l=1}^{T^{\prime}}(|\tilde{\mathbf{K}}_{1,2l-1}|+|\tilde{\mathbf{K}}_{2,2l}|)$
is the sum of the diagonal parts of $N\times N$ blocks of $\tilde{\mathbf{K}}$
and
$s_{\operatorname{n}}=\sum_{l=1}^{T^{\prime}}(|\tilde{\mathbf{K}}_{1,2l}|+|\tilde{\mathbf{K}}_{2,2l-1}|)$
is the sum of the nondiagonal parts of them.

### G.3. Real-world sleep data

Regarding the data, the original sample frequency was 10Hz. To make the data
smooth and remove the effect of noise, we interpolated observations of both
RIP and airflow linearly to obtain interpolated observations with a sample
frequency of 25Hz. We separate each time series into $5$ series and obtain
$\mathbf{x}^{i}_{0},\mathbf{x}^{i}_{1},\ldots,\mathbf{x}^{i}_{2000}\in\mathcal{X}^{N}=(\mathbb{R}^{2})^{2}$
for $i=1,\ldots S$, where $S=5$, with time interval $\Delta t=0.04$. We
obtained Eq. (6) in the same manner as the case of the Stuart–Landau model.
However, for $c^{j}_{1},c^{j}_{2}\in\mathbb{C}$, we computed the values such
that
$S^{-1}\sum_{k=1}^{S}\sum_{i=1}^{N}c^{j}_{i}v^{j}_{i}(x^{k}_{l})=\mathbf{1}$.
Here, we consider the average of $v^{j}_{i}(x^{k}_{l})$ over $k$ rather than
the original $v^{j}_{i}(x^{k}_{l})$ for some $k$ for making use of the all $S$
sequences of data.

### G.4. Real-world body temperature data

The body temperatures of five mice were recorded every 15 minutes for 85 days
(2040 hours). We averaged every consecutive 4 observations to obtain the
averaged body temperatures in every 1 hour. Then, we separated each time
series into 4 series with 501 averaged observations. We obtained
$\mathbf{x}^{i}_{0},\mathbf{x}^{i}_{1},\ldots,\mathbf{x}^{i}_{500}\in\mathcal{X}^{N}=(\mathbb{R}^{2})^{5}$
for $i=1,\ldots S$, where $S=4$, with time interval $\Delta t=1(h)$.

By solving the multiparameter eigenvalue problem (13), we estimated the
strength of interactions by obtaining the magnitude of $a_{i,k}$
$(i,k=1,\ldots,5)$ in Eq. (6). There are sets of three eigenvalues
$\lambda^{j}_{1},\ldots,\lambda^{j}_{5}\in\mathbb{C}$ of the estimated Koopman
operator that satisfy $|\lambda_{k}^{j}|\approx 1$ and
$\lambda_{1}^{j}\approx\lambda_{2}^{j}\approx\lambda_{3}^{j}\approx
e^{\sqrt{-1}j\omega_{0}\Delta t}$ for some $\omega_{0}\in[0,2\pi)$ and for
$k=1,\ldots,5$ and $j=1,2,\ldots$. Since only the magnitude of $a_{i,k}$ is
required here, we fixed $\lambda$ and $\mathbf{u}$ and simplified the problem.
We set $\lambda$ and $\mathbf{u}$ in the same manner as the initial values for
$N=2$. Then, we obtained $a_{i,k}$ that minimizes the minimization problem
(14).