Koopman Operator and its Approximations for Systems with Symmetries

TL;DR

This paper explores how symmetries in nonlinear dynamical systems influence the structure and computation of the Koopman operator, enabling more efficient approximations and revealing hidden organizational features.

Contribution

It demonstrates how symmetry considerations can simplify the computation of Koopman operator approximations using EDMD and kernel DMD methods, and reveals the block diagonal structure induced by symmetries.

Findings

Symmetries induce a block diagonal structure in Koopman operator approximations.

Symmetry considerations can improve the efficiency of EDMD and kernel DMD methods.

The paper discusses the impact of measurement noise on these methods.

Abstract

Nonlinear dynamical systems with symmetries exhibit a rich variety of behaviors, including complex attractor-basin portraits and enhanced and suppressed bifurcations. Symmetry arguments provide a way to study these collective behaviors and to simplify their analysis. The Koopman operator is an infinite dimensional linear operator that fully captures a system's nonlinear dynamics through the linear evolution of functions of the state space. Importantly, in contrast with local linearization, it preserves a system's global nonlinear features. We demonstrate how the presence of symmetries affects the Koopman operator structure and its spectral properties. In fact, we show that symmetry considerations can also simplify finding the Koopman operator approximations using the extended and kernel dynamic mode decomposition methods (EDMD and kernel DMD). Specifically, representation theory allows us to demonstrate that an isotypic component basis induces block diagonal structure in operator approximations, revealing hidden organization. Practically, if the data is symmetric, the EDMD and kernel DMD methods can be modified to give more efficient computation of the Koopman operator approximation and its eigenvalues, eigenfunctions, and eigenmodes. Rounding out the development, we discuss the effect of measurement noise.

Figures

Equations

$$\dot{x}=g_c(x).$$

$$\Phi (x(t),\Delta t)=x(t)+\underset{t}{\overset{t+\Delta t}{\int }}{g}_c(x(\tau ))d\tau .$$

$$x_{i+1}=g(x_i).$$

$$(\mathcal{K}f)(x)=f(g(x)).$$

$$(\mathcal{K}\varphi )(x)=\lambda \varphi (x).$$

$$x=\sum _i{v}_i{\varphi }_i(x),$$

$$\underset{A}{\int }\mathcal{P}\rho (x)dx=\underset{g^{-1}(A)}{\int }\rho (x)dx.$$

$$K={\Psi }_x^{+}{\Psi }_y.$$

$$\begin{array}{c}K=G^{+}A,\text{ where}\\ G=\sum _m\Psi (x_m{)}^{\ast }\Psi (x_m),A=\sum _m\Psi (x_m{)}^{\ast }\Psi (y_m).\end{array}$$

$${\varphi }_j(x)=\Psi (x)u_j.$$

$$v_i=(w_i^{\ast }B{)}^T.$$

$$k(x_i,y_i)=\Psi (x_i)\Psi (y_i{)}^{\ast }.$$

$$g_c({\gamma }_{\rho }x(t))={\gamma }_{\rho }{g}_c(x(t)).$$

$$g({\gamma }_{\rho }{x}_i)={\gamma }_{\rho }g(x_i).$$

$$({\gamma }_{\rho }\circ f)(x)\equiv f({\gamma }_{\rho }^{-1}x).$$

$$\ddot{x}=-\sigma \dot{x}-x(\beta +{\alpha }^{2}x).$$

$$\begin{array}{c}\dot{x}=y,\\ \dot{y}=-\sigma y-x(\beta +{\alpha }^{2}x).\end{array}$$

$${\gamma }_{s}g(\mathbf{x})=g({\gamma }_s\mathbf{x}).$$

$$\begin{array}{c}\dot{x_i}=y_i,\\ \dot{y_i}=-\sigma \dot{y_i}-x_i(\beta +{\alpha }^2{x}_i)+\sum _{ij}{\eta }_{ij}(x_i-x_j).\end{array}$$

$${\gamma }_{\rho }\circ (\mathcal{K}f)(x)=\mathcal{K}({\gamma }_{\rho }\circ f)(x).$$

$$\begin{array}{c}{\gamma }_{\rho }\circ (\mathcal{K}f)(x)={\gamma }_{\rho }\circ f(g(x))=f({\gamma }_{\rho }^{-1}g(x))\\ =f(g({\gamma }_{\rho }^{-1}x))=\mathcal{K}f({\gamma }_{\rho }^{-1}x)=\mathcal{K}({\gamma }_{\rho }\circ f)(x)\end{array}$$

$$\mathcal{K}({\gamma }_{\rho }\circ \varphi (x))={\gamma }_{\rho }\circ (\mathcal{K}\varphi )(x)=\lambda {\gamma }_{\rho }\circ \varphi (x).$$

$$\psi =\sum _p\sum _{m,n=1}^{d_p}{\alpha }_{mn}^p{\xi }_{mn}^p.$$

$${\mathcal{P}}_{mn}^p=\frac{d_p}{|\Gamma |}\sum _{\gamma \in \Gamma }[R_p(\gamma ){]}_{mn}^{\ast }{\gamma }_{\rho }.$$

$${\xi }_{mn}^p(x)=\frac{1}{c_{np}}{\mathcal{P}}_{mn}^p\circ \psi (x).$$

$${\mathcal{P}}^p=\frac{d_p}{|\Gamma |}\sum _{\gamma \in \Gamma }{\chi }_p(\gamma {)}^{\ast }{\gamma }_{\rho }.$$

$$\{({\gamma }_{\rho }{x}_i,{\gamma }_{\rho }{y}_i)\}=\{(x_i,y_i)\}.$$

$$\begin{array}{c}\Psi (x)=\left(\begin{array}{c}{\gamma }_1\circ \Psi (x)...{\gamma }_{|\Gamma |}\circ \Psi (x)\end{array}\right),\\ {\Psi }_x=\left(\begin{array}{ccc}{\Psi }_1(x_1)& ...& {\Psi }_{N/|\Gamma |}(x_1)\\ ⋮& \ddots & ⋮\\ {\Psi }_1(x_M)& ...& {\Psi }_{N/|\Gamma |}(x_M)\end{array}\right).\end{array}$$

$$P_{{\gamma }_k}({\gamma }_1,...,{\gamma }_{|\Gamma |}{)}^T=({\gamma }_k{\gamma }_1,...,{\gamma }_k{\gamma }_{|\Gamma |}{)}^T.$$

$${\mathbf{P}}_{{\gamma }_k}G=G{\mathbf{P}}_{{\gamma }_k},{\mathbf{P}}_{{\gamma }_k}A=A{\mathbf{P}}_{{\gamma }_k}.$$

Full text

Many natural and engineered dynamical systems — power grid networks and biological regulatory networks, to mention two — exhibit symmetries in their connectivity structure and in their internal dynamics. Some have time-reversal symmetry, others rotational and spatial translation invariance, and others still, combinations. These symmetries are often key for understanding the behavior of systems. For instance, a well-known example that exploits the connection between system state and interconnection symmetries arises in locomotion, where spatio-temporal symmetries of the observed locomotion behavior constrain the structure of neural circuits that generate these patterns. For network systems in particular, symmetries in the connectivity structure are of fundamental importance. For instance, the structural symmetries of a network of identical oscillators can determine its admissible patterns of symmetry-breaking. That said, additional information beyond the knowledge of the network structure is often required to address more detailed questions about a system’s dynamics, such as whether a particular configuration is stable in a given parameter regime. In these cases, the system’s linearization near the steady state can be combined with interconnection symmetry to provide the answer. However, these linearization methods are only valid on a local subset of the state space and therefore are not sufficient for global characteristics of nonlinear dynamical systems, such as their attractors, basins, and transients. The Koopman operator, in contrast, is a linear infinite-dimensional operator that evolves the functions on the state space that is valid on the entire state space. We show how to combine symmetry considerations with the Koopman analysis to study nonlinear dynamical systems with symmetries. We use representation theory to determine the effect of symmetries on the Koopman operator and its approximations, drawing out how local dynamical symmetries interact with symmetries arising from the connectivity of system variables. This, in turn, allows us to modify data-driven Koopman operator approximation algorithms to make them more efficient when applied to dynamical systems with symmetries. We illustrate our findings in a simple network of coupled Duffing oscillators with symmetries in individual oscillator dynamics and in their physical couplings.

I Introduction

Symmetries of dynamical systems manifest themselves in asymptotic dynamics, bifurcations, and attractor basin structure. Symmetries play a crucial role in guiding the emergence of synchronization and pattern formation, behaviors broadly observed in natural and engineered systems. Methods from group theory, representation theory, and equivariant bifurcation theory provide useful tools to study the common features of systems with symmetries golubitsky2003symmetry ; golubitsky2012singularities ; golubitsky2015recent .

Dynamical elements organized into a network are an important class of dynamical systems that often exhibit these behaviors, especially when symmetries appear in both network structure and the dynamics of the individual nodes. Studying the effect of symmetries in network topology of synthetic and real-life systems using computational group theory is an active area of research macarthur2008symmetry ; pecora2014cluster ; pecora2017discovering . Those symmetries lead to phenomena such as full synchronization, cluster synchronization, and formation of exotic steady states such as chimeras sorrentino2016complete ; cho2017stable ; matheny2019exotic ; emenheiser2016patterns . Moreover, topological symmetries underlying cluster synchronization of coupled identical elements assist in analyzing the stability of these fully synchronous cluster states pecora2014cluster ; hart2019topological . For networks of identical coupled oscillators, the form of their limit cycle solutions and the form of their bifurcations can be derived from symmetry considerations ashwin1992dynamics . Symmetries are also key in determining network controllability and observability. For example, Refs. rubin1972controllability ; dellnitz2017sensing explored the effect of explicit network symmetries for linear time-independent and time-dependent networks. Similarly, Refs. letellier2002investigating ; whalen2015observability considered nonlinear network motifs with symmetries and studied how the presence of different types of structural symmetries affect the observability and controllability of the system. Ref. mesbahi2019nonlinear , similar to our approach, uses the Koopman operator formalism (discussed below). They provide analytic results that link the presence of permutational symmetries in dynamical systems to their observability properties.

Many dynamical systems of current interest are high dimensional and nonlinear. For instance, this is the case for many complex networks, such as power grids and biological networks. Complexity there arises from the interaction between network interconnectivity structure and the nonlinearities in the node and edge dynamics. And, this often leads to multistability. Linearization methods can provide insight near the system’s attractors, but they poorly approximate the dynamics on the rest of the state space. In contrast, operator-based methods give access to the global characteristics of nonlinear systems. And, they do so in a linear setting and are therefore more well-suited, for instance, to characterize the attractor basin structure of multistable dynamical systems or the design control interventions. The Perron-Frobenius and Koopman operator are adjoint linear infinite-dimensional operators whose spectra can provide global information about the system. Their approximations using data-driven approaches make operator methods potentially applicable when there is no prior knowledge of the system.

The Perron-Frobenius operator evolves densities on state space. It has been used extensively to access global behavior of nonlinear dynamical systems lasota1985probabilistic ; tantet2018crisis . There are several well developed approaches for obtaining its numerical approximations, such as Ulam’s method that relies on the discretization of physical space to obtain an approximation of the Perron-Frobenius operator ulam1960collection . Since the Koopman operator is adjoint to the Perron-Frobenius operator, numerical approximations of the Koopman operators can be obtained using these methods as well klus2015numerical .

The Koopman operator is an infinite dimensional linear operator that describes the evolution of observables (functions of the state space) koopman1931hamiltonian ; lasota1985probabilistic ; klus2015numerical . Its definition and properties in the context of dynamical systems are provided, for instance, in Ref. budivsic2012applied , which also summarizes its applicability to model reduction, coherency analysis, and ergodic theory. Methods based on the Koopman operator decomposition have been proven useful for applications such as model reduction and control of fluid flows bagheri2013koopman , power system analysis susuki2016applied , and extracting spatio-temporal patterns of neural data brunton2016extracting .

Data driven methods to approximate the Koopman operator rely upon snapshot pairs of measurements of the state of the system at consecutive time steps. Reconstructing the operator from these snapshot pairs requires that a set of functions (called a dictionary of observables) is chosen. The first data driven method introduced, dynamic mode decomposition (DMD), implicitly uses linear monomials as a dictionary and thus is most applicable to systems where the Koopman eigenfunctions are well represented by this basis set schmid2010dynamic . A more recent method called extended DMD (EDMD) introduced in Ref. williams2015data can be more powerful than the standard DMD when applied to nonlinear systems as it allows the choice of more complicated sets of dictionary functions. Applying the EDMD is most computationally feasible if the number of the dictionary functions does not exceed the total number of the snapshot pairs used. That is not necessarily the case if a rich function dictionary (e.g., a dictionary of high order polynomials) is chosen. A modification of EDMD called kernel DMD introduced in Ref. williams2014kernel addresses this issue by providing a way to efficiently calculate the Koopman operator approximation in a case when the number of dictionary functions exceeds the number of measurements. Yet, the principled choice of an underlying dictionary that leads to an accurate approximation of the eigenspectrum corresponding to the leading Koopman modes using EDMD or kernel DMD remains an outstanding challenge. That problem is confronted, for instance, in Ref. li2017extended , where an iterative EDMD dictionary learning method is presented. Although the optimal choice of dictionary functions is often unknown, there are some common choices that are known to produce accurate results for certain classes of systems williams2015data .

Here we study nonlinear dynamical systems with discrete symmetries combining operator-based approaches and linear representation theory. Recently, related methods have been applied to dynamical systems with symmetries. On the one hand, Ref. mehta2006symmetry addresses symmetries of the Perron-Frobenius operator in relation to the admissible symmetry properties of attractors. On the other, Ref. sharma2016correspondence links the spatiotemporal symmetries of the Navier-Stokes equation to the spatial and temporal Koopman operator. Additionally, Ref. brunton2016discovering noted that symmetry considerations play an important role in discovering governing equations. And, Ref. kaiser2018discovering shows how conservation laws can be detected with Koopman operator approximations and then used to control Hamiltonian systems.

In contrast, our focus is on dynamical systems with symmetries described by a finite group. We show how the properties of the associated Koopman operator spectrum can be linked to the properties of the spectrum of the finite dimensional approximations of the Koopman operator obtained from finite data. We further show how the analytic properties of the Koopman operator decomposition can inform the choice of dictionary functions that can be used in the Koopman operator approximations. This gives a practical way to reduce the dimensionality of the approximation problem.

Our development builds as follows. Section II defines the Koopman operator, introduces approximation methods (EDMD and kernel DMD), and defines equivariant dynamical systems as well as useful concepts from group theory and representation theory. Section III draws out the implications of dynamical system symmetries for the structure of the Koopman operator and its eigendecomposition. Section IV connects the properties of the Koopman operator and the structure of its EDMD approximation for symmetric systems. This then allows modifying the EDMD method to exploit the underlying symmetries, resulting in a block-diagonal Koopman operator approximation matrix. We also provide numerical examples, showing how using particular dictionary structures speed up the algorithm. Finally, the last section summarizes our findings and outlines directions for future work.

II Preliminaries

II.1 Koopman operator

In this section, we provide some background in operator theoretic approaches to dynamical systems, in particular, the Koopman operator and its adjoint Perron-Frobenius operator. Since in this manuscript we address the approximations of the Koopman operator where the input is discrete time data, we focus on their definition in the discretized setting. The continuous time definitions are similar budivsic2012applied . Our results regarding the degeneracy of Koopman operator eigenvalues and the properties of its corresponding eigenfunctions presented in Section III hold in both discrete and continuous time settings.

Suppose we are interested in studying continuous time autonomous dynamical systems defined as:

[TABLE]

Here, $x\in\mathbb{R}^{n}$, $g_{c}:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$. Let $\Phi(x(t),\Delta t)$ be a flow map mapping the initial condition $x(t)$ to the solution at time $t+\Delta t$. It is defined in the following way:

[TABLE]

The system can be discretized with a finite time step $\Delta t_{\text{step}}$, so that $x_{i+1}=\Phi(x_{i},\Delta t_{\text{step}})$. We denote the function evolving the dynamics of this discretized system by $g$:

[TABLE]

The Koopman operator is a linear infinite dimensional operator that evolves the functions (referred to as observables) of state space variables $f:\mathbb{R}^{n}\rightarrow\mathbb{C}$. The action of the Koopman operator $\mathcal{K}$ on an observable function $f$ for discrete time systems is defined as:

[TABLE]

Since we consider data-driven Koopman operator approximation methods in this manuscript, the discrete time version of the definition is most applicable.

In general, parts of the Koopman operator spectrum can be continuous korda2018data ; budivsic2012applied , for instance, that can be the case for chaotic systems. However, for the purposes of this manuscript, we focus on the case of a discrete spectrum.

Pairs of eigenvalues $\lambda$ and eigenfunctions $\phi$ of the Koopman operator $\mathcal{K}$ are defined as:

[TABLE]

Of particular interest are the Koopman modes that can be used in model reduction and coherency estimation susuki2011nonlinear ; rowley2017model . The Koopman modes $v_{i}$ of the full state observable are defined by:

[TABLE]

and are projections of the observable onto the span of the eigenfunctions of the Koopman operator $\mathcal{K}$.

The other candidate for studying dynamical systems using an operator based approach is the Perron-Frobenius operator $\mathcal{P}$ defined as follows for deterministic dynamical systems:

[TABLE]

Here, $\rho(x)$ is a density on state space, and $A\subseteq\mathbb{R}^{n}$ is a subset of the state space. The Perron-Frobenius operator is an adjoint to the Koopman operator lasota1985probabilistic , so an approximation of one of them provides an approximation of the other klus2015numerical .

II.2 Koopman operator approximation methods

Extended dynamic mode decomposition (EDMD) introduced in Ref. williams2015data is a way of approximating the Koopman operator for discretized systems that requires an explicit choice of a dictionary of functions referred to as observables. How to optimally choose those functions remains an open problem for many systems, especially if the form of differential equations describing the governing dynamical system is not known in advance and only finite data on the behavior of the system is available. The method can be very accurate in capturing the dynamics of the system, but its accuracy depends strongly on the choice of an appropriate dictionary of observables. The method’s convergence properties are studied in Ref. korda2018convergence , and its relation to the Perron-Frobenius operator approximation methods is discussed in Ref. klus2015numerical . Here, we summarize the EDMD and its relation to the Koopman operator.

The first requirement for the method is a set pairs of consecutive snapshots $x=[x_{1},x_{2},...,x_{M}]$ and $y=[y_{1},y_{2},...,y_{M}]$, where the measurements $x_{i}$ and $y_{i}$ are performed with a small constant time interval $\Delta t$: $y_{i}=\Phi(x_{i},\Delta t)$. Typically, the set of snapshots contains measurements from different trajectories in state space. We define a dictionary of linearly independent observables $\mathcal{D}=\{\psi_{1},...,\psi_{N}\}$ and form vectors of observables $\Psi_{x}$ and $\Psi_{y}$. Here, $\Psi_{x}\in\mathbb{R}^{M\times N}$, where $N$ is the number of dictionary functions used in the approximation, and $M$ is the number of data snapshots. The elements of $\Psi_{x}$ are obtained from $(\Psi_{x})_{ij}=\psi_{j}(x_{i})$. We also use the notation $\Psi(x_{m})=(\psi_{1}(x_{m}),...,\psi_{N}(x_{m}))$ for the dictionary functions evaluated at a particular point on the trajectory.

A finite dimensional approximation of the Koopman operator $\mathcal{K}$ that we denote as $K$ can be obtained using:

[TABLE]

Here, $\Psi_{x}^{+}$ denotes the pseudoinverse of $\Psi_{x}$. We focus on the case of the Moore-Penrose pseudoinverse for the rest of the manuscript penrose1955generalized .

If the number of snapshots is much higher than the dimensionality of the function dictionary ($M\gg N$), it is more practical instead to define the square matrices $G$ and $A$ as shown below and obtain the approximation in the following way:

[TABLE]

Here, ∗ represents the complex conjugate transpose. If the only observables are the states of the system $x_{1},x_{2},...,x_{n}$, EDMD reduces to DMD williams2015data ; klus2015numerical .

The eigendecomposition of $K$ provides the Koopman eigenvalues, eigenfunctions, and modes that allow an approximate linear representation of the underlying system dynamics. Let $\lambda_{j}$ and $u_{j}$ be the $j$th eigenvalue and eigenvector of $K$. Then the corresponding Koopman eigenfunction can be approximated by:

[TABLE]

Let $b_{i}$ be the vectors defined by $g(x)_{i}=\Psi b_{i}$, where $g(x)_{i}=e^{*}_{i}x$ denotes the elements the full state observable discussed in Ref. williams2015data , and $B=(b_{1}~{}~{}...~{}~{}b_{n})$. The Koopman eigenmodes can then be obtained as:

[TABLE]

Here, $w_{i}$ denotes the $i$th left eigenvector of $K$.

A modification of EDMD named kernel DMD williams2014kernel is better suited for systems with a low number of measurements and a high number of observables (e.g., the full state observable for fluid dynamical systems is very high dimensional, so defining a polynomial dictionary of the full state observable is very computationally expensive), i.e. $M\ll N$. The method relies on evaluating the kernel function:

[TABLE]

That allows efficient computation of $M\times M$ matrices $\hat{G}$, $\hat{A}$, and $\hat{K}$, where $M$ is the number of trajectory time steps. The eigendecomposition of $\hat{K}$ then can be used to obtain the approximations of the Koopman eigenvalues, eigenfunctions, and modes.

In the main body of the manuscript, we focus on the case where the number of measurements is relatively high for each degree of freedom ($M\gg N$), and obtain a way to reduce the dimensionality of the EDMD approximation of the Koopman operator for systems with symmetries in Section IV. A similar modification of the kernel DMD is discussed in Appendix D.

II.3 Point symmetries

In this section, we define the concepts useful to study the structure of the Koopman operator $\mathcal{K}$ and its approximations $K$ for systems with symmetries. Throughout this section and the rest of the manuscript, we use an example of a small network of Duffing oscillators to illustrate the definitions and algorithms.

In this manuscript, we consider dynamical systems (as defined in Eqs. 1 and 3) that respect point symmetries. These systems are called equivariant with respect to the symmetry group $\Gamma$. We define groups by their presentations in a form $\langle S|R\rangle$, where $S$ is a set of generators of the group, and $R$ is a set of relations among these generators that define that group. Every element of the group can be written as a product of powers of some of these generators.

For instance, the cyclic group $Z_{3}$ is presented by $\langle r|r^{n}=1\rangle$. An example of a realization of that group is a set of rotational symmetries of a regular $n$-gon.

To study dynamical systems with symmetries, we need to define the specific actions of the group on a vector space in addition to an abstract presentation of a group $\Gamma$. Let $X\subset\mathbb{R}^{n}$ be a vector space with elements $x\in X$. We denote the actions $\gamma_{\rho}$ on a vector space $X$ by $\gamma_{\rho}x$ if the set of these actions $\Gamma_{\rho}$ is isomorphic to $\Gamma$. A shorthand $\gamma_{\rho}x=\gamma x$ is sometimes used in the literature when the action corresponding to the subscript $\rho$ is clear from the context (for instance, it is defined by a permutation matrix of the same degree as the state space of the system), however, we use the $\gamma_{\rho}$ notation to avoid ambiguity, since the precise definition of group action in particular cases is important in this manuscript, as shown, for instance, in Example II.1 and Example II.2.

Finally, we define what it means for a dynamical system to be symmetric. Let $\dot{x}=g_{c}(x)$ be a continuous time system of differential equations. Here, $x\in\mathbb{R}^{n}$, and $g_{c}:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$. The system is $\Gamma$-equivariant with respect to the actions of $\Gamma_{\rho}$ if for all $x\in\ X$ and $\gamma_{\rho}\in\Gamma_{\rho}$:

[TABLE]

As discussed in Section II, data comes in discretized form, so a discrete form of that definition is useful. For discrete time systems defined by $x_{i+1}=g(x_{i})$, equivariance is defined in a similar manner:

[TABLE]

We note that if a continuous time system is $\Gamma$-equivariant, so is its discretization. Moreover, the set of trajectories of a $\gamma$-equivariant system in state space also respects the symmetries of the system. For discretized systems, it means that if $\{x_{0},x_{1},...x_{n}\}$ form a trajectory in state space, then $\{\gamma_{\rho}x_{0},\gamma_{\rho}x_{1},...\gamma_{\rho}x_{n}\}$ form a trajectory as well.

An important example of equivariant dynamical systems that many of the recent works have focused on (such as Refs.pecora2014cluster ; hart2019topological ; sorrentino2016complete ; matheny2019exotic ) is a system of coupled identical oscillators. In that case, the set (or a subset) of actions under which the system is equivariant is defined by the set of permutational matrices $P$ that commute with the adjacency matrix of that oscillator network. In this case, the action of the group is linear, however, that does not always have to be the case.

We also need to define the action of the group in function space, where $f\in\mathcal{F}$ are functions $f:X\rightarrow\mathbb{C}$ as:

[TABLE]

This definition will be useful since the Koopman operator acts on functions (i.e. observables).

Another concept useful for our work is a linear group representation $T$, which is a mapping from group elements $\gamma\in\Gamma$ to the elements of the general linear group (a group of matrices of degree $n$ with the operation of matrix multiplication denoted by $\text{GL}(n,V)$) on a vector space $V$ (in this case, we are interested in $V=\mathbb{C}^{n}$). The characters of a group representation $T_{i}(\gamma)$ are defined as $\chi_{i}(\gamma)=\text{Tr}(T_{i}(\gamma))$.

A representation is called irreducible if it has no nontrivial invariant subspaces (meaning that the representation matrices corresponding to the group elements can not be simultaneously non-trivially block diagonalized into the same block form). For each $\Gamma$ we can obtain all of its irreducible matrix representations. We denote their elements mapping $\gamma\in\Gamma$ to $p\times p$-dimensional matrices as $R_{i}(\gamma)$, where the index $i$ corresponds to the $i$th irreducible representation. Irreducible representations are defined up to an isomorphism. For the purposes of this manuscript, it is useful to either make use of the unitary irreducible representations or the characters $\chi_{i}(\gamma)$ of the unitary irreducible representations.

A vector space, e.g. the space of square integrable functions $\mathcal{F}$, can be uniquely decomposed into components that transform like the $i$th irreducible representation of $\Gamma$ under the actions of $\Gamma_{\rho}$. These components are called isotypic components golubitsky2003symmetry . We denote these components by $\mathcal{F}_{i}$. An isotypic decomposition of the square integrable function space with respect to $\Gamma_{\rho}$ is then defined as $\mathcal{F}=\bigoplus\limits_{i}\mathcal{F}_{i}$. We illustrate the construction of an isotypic decomposition using an example of a $Z_{2}$-equivariant system.

Example II.1

Symmetries of a single Duffing oscillator dynamics and isotypic components in function space corresponding to the actions of its symmetry group.

The unforced Duffing oscillator equation has the form:

[TABLE]

We can rewrite the above equation as a system of differential equations to obtain:

[TABLE]

Let $\bm{x}=\begin{pmatrix}x\\ y\end{pmatrix}$, and let the dynamics be denoted by $\dot{\bm{x}}=g(\bm{x})$. Let $r_{s}=\begin{pmatrix}-1&0\\ 0&-1\end{pmatrix}=-I_{2\times 2}$ be the action on the state space that flips the signs of variables. The actions $r_{s}$ and $e_{s}=I_{2\times 2}$ form a group $\Gamma_{s}$ isomorphic to $\Gamma=Z_{2}=\langle r|r^{2}=e\rangle$. Let $\gamma_{s}\in\Gamma_{s}$, then:

[TABLE]

Thus, the Duffing oscillator system is $Z_{2}$-equivariant with respect to the actions $\gamma_{s}$.

We now illustrate the isotypic component decomposition of $Z_{2}$ in function space. $Z_{2}$ has two one-dimensional irreducible representations: the trivial representation defined by $R_{\text{tr}}(r)=1$ and the sign representation defined by $R_{\text{sign}}(r)=-1$. Then the space of square integrable functions $\mathcal{F}$ can be decomposed into $\mathcal{F}=\mathcal{F}_{\text{tr}}\oplus\mathcal{F}_{\text{sign}}$, where $\mathcal{F}_{\text{tr}}=\{f:r_{s}\circ f=f(-x,-y)=f(x,y)\}$ and $\mathcal{F}_{\text{sign}}=\{f:r_{s}\circ f=f(-x,-y)=-f(x,y)\}$. In this case, the sets of functions $\mathcal{F}_{\text{tr}}$ and $\mathcal{F}_{sign}$ consisting of even and odd functions respectively transform like the trivial and sign irreducible representations with respect to sign flip as a group generator action.

We now extend the example to a network of Duffing oscillators and explore additional permutation symmetries.

Example II.2

We now consider the dynamics of a network of Duffing oscillators. Suppose the coupling is linear in $x$ with a coupling coefficient assigned to every edge $\eta_{ij}$. Then for each node $i$ in the network, we have the following dynamics:

[TABLE]

This general coupling scheme is used to describe many systems in literature pecora2014cluster ; sorrentino2016complete .

We now consider the case of a 3-node network. Depending on what the coupling terms are, the system may be $\Gamma$-equivariant with respect to different symmetry groups that act by permuting node indexes. Some examples are:

1. $a)$

If all coupling strengths $\eta_{ij}$ are equal, the network has $D_{3}$ symmetry. This case is shown on Fig. 1(a). Let the state of the system be defined by $x=\left(x_{1}~{}~{}y_{1}~{}~{}x_{2}~{}~{}y_{2}~{}~{}x_{3}~{}~{}y_{3}\right)^{T}$. Then, the symmetry group is presented by $D_{3}=\langle r,\kappa|r^{3}=\kappa^{2}=e,\kappa r\kappa=r^{-1}\rangle$ and generated by the actions $r_{p}=\begin{pmatrix}0&1&0\\ 0&0&1\\ 1&0&0\end{pmatrix}\otimes I_{2\times 2}$ and $\kappa_{p}=\begin{pmatrix}1&0&0\\ 0&0&1\\ 0&1&0\end{pmatrix}\otimes I_{2\times 2}$. 2. $b)$

If the coupling strengths obey the conditions $\eta_{ij}\neq\eta_{ji}$ and $\eta_{ij}=\eta_{jk}$, the network has $Z_{3}$ symmetry. This case is shown on Fig. 1(b). The symmetry group is presented by $Z_{3}=\langle r|r^{3}=e\rangle$ and generated by the action $r_{p}$ defined above. 3. $c)$

If the coupling strengths obey the conditions $\eta_{12}=\eta_{21}$ and no other equalities hold, the network has $Z_{2}$ symmetry. This case is shown on Fig. 1(c). The symmetry group is presented by $Z_{2}=\langle\kappa|\kappa^{2}=e\rangle$ and generated by the action $\kappa_{p}$ defined above.

Even though in case (c) the permutation symmetry is isomorphic to the same group as the sign flip symmetry in Example II.1, the isotypic components in function space $\mathcal{F}$ induced by the group action are different. $Z_{2}$ has two one-dimensional irreducible representations: trivial representation $R_{1}(r)=R_{tr}(r)=1$ and sign representation $R_{2}(r)=R_{sign}(r)=-1$. Let $\bm{x_{i}}=\begin{pmatrix}x_{i}\\ y_{i}\end{pmatrix}$. The isotypic components are defined by $\mathcal{F}_{tr}=\{f:r_{p}\circ f=f(\bm{x}_{1},\bm{x}_{3},\bm{x}_{2})=f(\bm{x}_{1},\bm{x}_{2},\bm{x}_{3})\}$ and $\mathcal{F}_{sign}=\{f:r_{p}\circ f=f(\bm{x}_{1},\bm{x}_{3},\bm{x}_{2})=-f(\bm{x}_{1},\bm{x}_{2},\bm{x}_{3})\}$.

Additionally, each node still has $Z_{2}$ symmetry with respect to the action $r_{s}$ which is not broken since the coupling function is odd. That symmetry is also depicted in Fig. 1. We can find the isotypic components of the whole symmetry group as $\mathcal{F}=\bigoplus\limits_{i,j}\left(\mathcal{F}_{i,s}\otimes\mathcal{F}_{j,p}\right)$.

Any function can be rewritten as a sum of projections into different isotypic components. The procedure is outlined in the following section.

III Properties of the Koopman operator for systems with symmetries

In this section, we consider the structure of the eigenspace of the Koopman operator of $\Gamma$-equivariant systems. We show how to obtain a particular eigenbasis of the system corresponding to the isotypic decomposition in function space and demonstrate that the isotypic decomposition induces a block diagonal structure on the matrix representation of $\mathcal{K}$.

Theorem III.1

For a $\Gamma$-equivariant dynamical system $\dot{x}=g(x)$ and an arbitrary function $f$, the Koopman operator commutes with the actions of the elements of $\Gamma$:

[TABLE]

**Proof:

The commutativity follows from the definitions of the Koopman operator and the definition of the action of the group in state space and function space.**

[TABLE]

This result is similar to Theorem 3.1 in mehta2006symmetry , where it is shown that the action of Perron-Frobenius operator commutes with the action of the symmetry group $\Gamma$ for $\Gamma$-equivariant systems.

Corollary III.1.1

*The space of eigenfunctions of the Koopman operator $\mathcal{K}$ with eigenvalue $\lambda$ for a $\Gamma$-equivariant system is $\Gamma$-invariant.

**Proof:

*Let $S_{\lambda}$ be the set of eigenfunctions of $\mathcal{K}$ with eigenvalue $\lambda$. Let $\phi\in S_{\lambda}$. Then, using the commutativity of $\mathcal{K}$ and $\Gamma_{\rho}$, we can show that $\forall\gamma_{\rho}\in\Gamma_{\rho}$:

[TABLE]

Thus, $\phi_{\gamma,\rho}\in S_{\lambda}$, where $\phi_{\gamma,\rho}$ is also an eigenfunction with an eigenvalue $\lambda$ defined as $\phi_{\gamma,\rho}(x)=\gamma_{\rho}\circ\phi(x)$.

We now consider a particular form of the eigenbasis of the Koopman operator that induces block diagonal structure of the matrix representation of the action of the Koopman operator $\mathcal{K}$. The result quoted below is useful for that purpose.

Theorem III.2

*(Theorem 3.5 in Chapter XII of Ref.golubitsky2012singularities ).

Let $\Gamma$ be a compact Lie group acting on the vector space $V$ decomposed into isotypic components $V=W_{1}\bigoplus...\bigoplus W_{t}$. Let $A:V\rightarrow V$ be a linear mapping commuting with $\Gamma$. Then $A(W_{k})\subset W_{k}$.*

This result is applicable to finite symmetry groups. Isotypic components of $\mathcal{F}$ with respect to $\Gamma$ induce block diagonal structure of the matrix representation of the Koopman operator. Since $\mathcal{K}$ and $\Gamma$ commute, $\mathcal{K}(\mathcal{F}_{k})\subset\mathcal{F}_{k}$. This block structure can be exploited in finding the Koopman operator approximations, as we show in the next section. Thus, we need to be able to obtain an isotypic component basis from an arbitrary function dictionary. This is a well defined procedure cornwell1997group , outlined below. Functions obtained via isotypic decomposition are useful to perform calculations in many areas of physics, for instance, they can simplify finding approximate solutions to Schrodinger equation, or in studying crystallographic point groups cornwell1997group ; stiefel2012group . The construction is also widely applied to dynamical systems, for instance, to study states and their stability using equivariant bifurcation theory.

Suppose we start from an arbitrary basis function dictionary $\mathcal{D}_{\psi}=\{\psi_{i}\}$. Each of those functions can be expanded in the isotypic component basis with at least one nonzero coefficient $\alpha^{p}_{mn}$:

[TABLE]

Here, $\xi^{p}_{mn}$ is a basis function in the $p$th isotypic component of $\mathcal{F}$, and $d_{p}$ is the dimension of that isotypic component. Alternatively, it can be thought of as a sum over all inequivalent (non-isomorphic) irreducible representations of $\Gamma$, where $\xi_{mn}^{p}$ transforms as $(m,n)$th element of $p$th irreducible representation of $\Gamma$ cornwell1997group . We define a projection operator and form a new function basis consisting of functions $\{\xi_{mn}^{p}\}$ as outlined below.

The projection operator is defined as:

[TABLE]

Here, $[R_{p}(\gamma)]_{mn}$ denotes the element in $n$th row and $m$th column of the $i$ th unitary irreducible representation of $\gamma\in\Gamma$, and $\gamma_{\rho}$ is the group action. We can form an orthonormal basis $\mathcal{D}_{\xi}=\{\xi_{i}\}$ using that projection operator as follows:

[TABLE]

Here, $c_{np}=\langle\mathcal{P}^{p}_{nn}\psi,\mathcal{P}^{p}_{nn}\psi\rangle^{1/2}$, where $\langle,\rangle$ denotes the inner product, which can be omitted for our purposes since the scaling of basis functions does not affect the EDMD results.

Equivalently, due to orthogonality of characters of irreducible representations, the projection operator can be obtained using the following expression:

[TABLE]

Here, $\chi_{p}(\gamma)$ is a character of the $p$th irreducible representation of $\Gamma$. If this formula is used, each irreducible representation of degree $d_{p}$ provides a basis function, and $d_{p}^{2}-1$ other basis functions can be formed using the Gram-Schmidt orthogonalization process cornwell1997group ; stiefel2012group ; pecora2014cluster .

Once an isotypic component basis is obtained, the action of the Koopman operator on function space can be presented in the form of a block diagonal matrix. Each irreducible unitary representation of dimension $d_{p}$ in this case corresponds to a number $d_{p}$ of $d_{p}\times d_{p}$ sized blocks in that matrix $\mathcal{K}$. Similar analysis works for the Koopman operator approximation $K$. The reason why this additional decomposition works can be found in Appendix A.

IV Implications for EDMD

IV.1 Constructing a basis for systems with known symmetries

In this section we show that the approximation of $K$ obtained using EDMD can be reduced to the block diagonal structure similar to $\mathcal{K}$ under certain assumption on the data. We provide some examples of constructing an isotypic component basis from a given function dictionary. We highlight that the basis depends on both the structure of $\Gamma$ and the definition of its actions $\Gamma_{\rho}$.

First, we establish that the Koopman operator approximation $K$ commutes with the actions $\gamma_{\rho}$ of $\Gamma$ if the data used in the calculation respects the symmetry, meaning the set of pairs of data points satisfies the condition:

[TABLE]

In other words, the set of trajectories is closed under the action of the symmetry group of the underlying dynamical system.

In order to perform further simplifications, we pick a particular order of group elements $\{\gamma_{1},...,\gamma_{|\Gamma|}\}$ and create the vectors $\Psi_{x}$ (and analogously $\Psi_{y}$) according to that ordering:

[TABLE]

Given the ordering of the group elements, we can also construct the permutation representation of the group such that:

[TABLE]

By Cayley’s theorem, such permutations form a group isomorphic to $\Gamma$. Determining the actions $P_{\gamma_{k}}$ of the group generators is sufficient to find the actions of all group elements. Let $\bm{P}_{\gamma_{k}}=P_{\gamma_{k}}\otimes I_{n\times n}$. We note that $(\bm{P}_{\gamma_{k}})^{*}=(\bm{P}_{\gamma_{k}})^{-1}$. It can be shown that:

[TABLE]

By definition, $A=\Psi_{x}^{*}\Psi_{y}$. We note that:

[TABLE]

That is the case for symmetric trajectories satisfying Eq. 35:

[TABLE]

Thus, $A$ and $G$ commute with the action of the symmetry group.

If $G$ is invertible and $G$ commutes with $\gamma_{\rho}$, $G^{-1}$ commutes with $\gamma_{\rho}$ as well. Then:

[TABLE]

If $G$ is not invertible, the commutativity result still holds for $G^{+}$. $G$ is a normal matrix since it satisfies $GG^{*}=G^{*}G$. In Appendix B, we show that if $G$ is normal, $GG^{+}=G^{+}G$, so $G$ commutes with its Moore-Penrose pseudoinverse, and therefore the actions of $K$ and $\Gamma$ commute.

Since $K$ commutes with the actions of $\Gamma$, $K\mathcal{F}_{i}\subset\mathcal{F}_{i}$. This shows that $K$ can be block-diagonalized in the same way as $\mathcal{K}$.

Suppose we start from a dictionary of observables. Since that dictionary is not necessarily an isotypic component dictionary corresponding to $\Gamma$ and its action, in order to obtain a block diagonal matrix $K$, the dictionary needs to be modified using the procedure outlined in Section III. In the example below, we show explicitly how to perform this transformation into the isotypic component basis.

In order for the basis to faithfully represent the symmetries of the system we require that:

• •

The dictionary is closed under the action of the symmetries of the system:

[TABLE]

• •

Each isotypic component is present after the isotypic component decomposition of the original function basis:

[TABLE]

For instance, using a monomial basis for a $D_{3}$ equivariant system does not satisfy the second requirement.

If these requirements are satisfied, the change of basis does not affect the result obtained by applying the EDMD algorithm as shown in Appendix C. Additionally, we note that the eigenvalues of $K$ do not typically have the same degeneracy properties as the eigenvalues of $\mathcal{K}$, but the symmetries of the underlying dynamical system are preserved in trajectory reconstructions.

Example IV.1

Constructing an isotypic component basis for a network of Duffing oscillators from a given basis.

We consider the case of a system of Duffing oscillators with identical coupling as depicted in Fig. 1(a). In that case, the system has $Z_{2}\otimes D_{3}$ symmetry. Suppose we want to construct an isotypic component basis from a given function dictionary $\mathcal{D}$. As an example, we use an initial dictionary $\mathcal{D}_{\psi}$ of $n$ mesh-free radial basis functions. The radial basis function centers can be obtained by either k-means clustering of the data or sampling from a predetermined distribution. Each function can be presented in a form $\psi(c,x)=r_{c}log(r_{c})$, where $c$ is a 6-dimensional radial basis function center, and $r_{c,x}=||x-c||^{1/2}$. In order to preserve the symmetries of the system, we need to have dictionary elements corresponding to acting on the basis functions by each $\gamma_{\rho}\in\Gamma_{\rho}$. Due to the form of these functions, $\gamma_{\rho}\circ\psi(c,x)=\psi(\gamma_{\rho}^{-1}c,x)$.

Since $\Gamma=Z_{2}\otimes D_{3}$ is a direct product of two groups, we can write the projection operator in the following form:

[TABLE]

We first consider the irreducible representations of $D_{3}$.

• •

Trivial representation $R_{tr}$: $R_{tr}(r)=1$, $R_{tr}(\kappa)=1$

• •

Sign representation $R_{sign}$: $R_{sign}(r)=1$, $R_{sign}(\kappa)=-1$

• •

Standard representation $R_{st}$: $R_{st}(r)=\begin{pmatrix}\omega&0\\ 0&\omega^{2}\end{pmatrix},~{}~{}R_{st}(\kappa)=\begin{pmatrix}0&1\\ 1&0\end{pmatrix}$. Here, $\omega=e^{2\pi i/3}$.

Suppose we form a vector of basis functions in $\mathcal{D}_{\psi}$, $\Psi=(\psi_{1,1}\psi_{2,1},...,\psi_{n,1},...,\psi_{1,|\Gamma|},\psi_{2,|\Gamma|},...,\psi_{n,|\Gamma|})^{T}$, where the first index corresponds to acting on $\psi_{1,i}$ by the $j^{th}$ element of $\Gamma_{\rho}$. Using the equations above, we obtain transformation matrices that we will use to get the isotypic component basis:

[TABLE]

The isotypic component basis then can be obtained by modifying a set of functions in $\mathcal{D}_{\psi}$:

[TABLE]

If we use $\Xi$ as a basis, we obtain $K$ decomposed into 8 blocks, each corresponding to an irreducible representation of $Z_{2}\otimes D_{3}$.

We implement the EDMD algorithm to obtain the approximation of $\mathcal{K}$. Here, the data comes from 500 initial trajectories of length 10 that were then reflected and rotated so that the data respects the symmetries. The parameter values of $\alpha=1.$, $\beta=-1.$, $\delta=.5$, and $\eta=1.$ were used. We plot the approximation matrix $K$ on Fig. 2. In this case, a dictionary of 120 radial basis functions was used. Fig. 2(a) illustrates the Koopman operator approximation matrix $K$ calculated using an initial dictionary $D_{\psi}$ and requires performing matrix operations on the full $120\times 120$ matrix. Fig. 2(b) shows $K$ obtained from the symmetry adapted basis functions. The order of calculations can be reduced significantly since it is only necessary to perform matrix operations on blocks. $K_{symm}$ has $4$, $10\times 10$ and $2$, $20\times 20$ unique blocks.

As shown in examples above, we can construct a basis that diagonalizes the Koopman operator matrix approximation $K$ from the elements of any arbitrary basis. Since the off block-diagonal elements of the matrix are a priori known to be zero, we do not need to compute these elements explicitly. This suggests that for systems with symmetries it is more efficient to perform the EDMD algorithm for isotypic decomposition blocks. We denote the number of conjugacy classes or irreducible representations of $\Gamma$ by $r_{\Gamma}$. In that case, instead of performing $O((mr_{\Gamma})^{\alpha})$ operations of matrix inversion, multiplication, and eigendecomposition, it is sufficient to perform these operations for each of the $r_{\Gamma}$ blocks, with operations being $O(m^{\alpha})$. Here $2<\alpha<3$ knuth1997art . Even though the algorithmic complexity only differs by a factor that scales with the size of the group that is fixed for any given system, in practice, the computation is more efficient when EDMD specific to $\Gamma$-equivariant systems is used. We also note that each $d_{p}$ dimensional irreducible representation results in $d_{p}$ of $d_{p}\times d_{p}$ dimensional blocks in $K_{\text{symm}}$ known to be equal a priori, which further simplifies the calculation. Moreover, in the case of networks of high dimensionality, it allows parallel eigendecomposition computation for isotypic component blocks. Table 1 summarizes the modified EDMD algorithm for $\Gamma$-equivariant systems and highlights that the order of computations can be lowered.

Koopman eigenfunctions and eigenmodes have many applications in dimensionality reduction, finding the basins of attraction, characterizing coherency between oscillatory systems, etc. Block diagonalizing $K$ allows the efficient computation of the Koopman eigenvalues, eigenfunctions, and modes.

The kernel DMD is closely related to the EDMD algorithm. It relies on the calculating the eigentriples associated with $K$ from a dual matrix $\hat{K}$ evaluated using a kernel trick commonly applied in machine learning williams2014kernel . This method can be computationally advantageous for cases when the number of basis functions exceeds the number of available measurements of the state of the system. We find that the kernel DMD can also be modified to include symmetry considerations in order to optimize the calculations. The method is provided in Appendix D.

IV.2 Consequences of symmetry assumptions in the basis

Assume the data is symmetric as defined by Eq. 35 with respect to the symmetry group $\Gamma$. A “perfect” basis is the one respecting the isotypic decomposition of $\Gamma$. Suppose the basis functions belong to isotypic components of $\Sigma\neq\Gamma$. That choice will affect the structure of $K$. We study that structure by evaluating the elements of $A$, since $K$ and $G^{+}$ have the same structure as $A$.

If the system is $\Gamma$-equivariant and $\Sigma\subset\Gamma$ and the set of actions of $\Sigma$ is a subset of actions of $\Gamma$, the system is also $\Sigma$-equivariant. Thus, picking a basis respecting the isotypic decomposition of $\Sigma$ will have the block diagonal structure corresponding to $\Sigma$. This means that choice of basis results in block diagonal $K$, but its structure does not provide any additional information about the symmetries of the system.

If the system is $\Gamma$-equivariant and $\Gamma\subset\Sigma$, functions belonging to particular isotypic components of $\Sigma$ are not preserved by the action of $K$. In the case of symmetric trajectories, that can provide information on what the true symmetries of the system are.

A simple case corresponds to $\Sigma=\Sigma_{0}\otimes\Gamma$. In this case, every action of $\Sigma_{0}$ commutes with every action of $\Gamma$. Each isotypic component of $\mathcal{F}$ with respect to $\Sigma$ can be expressed as $\mathcal{F}_{pq}=\mathcal{F}_{\Sigma_{0}}^{p}\cap(\mathcal{F}_{\Gamma})^{q}$, where $\mathcal{F}_{\Sigma_{0}}^{p}$ denotes the $p$th isotypic component of $\Sigma_{0}$. In this case, the off-diagonal blocks corresponding to interactions between isotypic components $\mathcal{F}_{p_{1}q_{1}}$ and $\mathcal{F}_{p_{2}q_{2}}$ are zero if $q_{1}=q_{2}$, and generally nonzero otherwise. For instance, if a network of three Duffing oscillators similar to the one discussed in, e.g., Example IV.1, has no permutation symmetry and $\Sigma=Z_{2}\otimes D_{3}$, with the action of $Z_{2}$ being a sign flip in nodal dynamics, the isotypic components corresponding to these $Z_{2}$ symmetries will not interact, resulting in two blocks in $K$.

Next, we consider a more general case. We denote the $p$th isotypic component of $\mathcal{F}$ with respect to the symmetry group $\Sigma$ by $\mathcal{F}^{p}_{\Sigma}$. We note that if $\mathcal{F}^{p}_{\Sigma}\cap\mathcal{F}^{q_{1}}_{\Gamma}\neq\varnothing$ and $\mathcal{F}^{p}_{\Sigma}\cap\mathcal{F}^{q_{2}}_{\Gamma}\neq\varnothing$ where $q_{1}$ and $q_{2}$ index different isotypic components of $\Gamma$, then the off-diagonal blocks of $K$ corresponding to interactions between those components are generally nonzero. The condition for $\mathcal{F}^{p}_{\Sigma}\cap\mathcal{F}^{q}_{\Gamma}\neq\varnothing$ is equivalent to:

[TABLE]

where $f$ is an arbitrary function, and $\mathcal{P}^{p}_{\Sigma}$ denotes the projection operator onto the $p$th isotypic component with respect to the symmetry group $\Sigma$.

[TABLE]

Let $H$ be the set of left cosets of $\Gamma$ in $\Sigma$ (defined as $H=\Sigma/\Gamma=\{\sigma\Gamma:\sigma\in\Sigma\}$, where $\sigma\Gamma=\{\sigma\gamma:\gamma\in\Gamma\}$ golubitsky2003symmetry ). Thus, the condition of Eq. 52 holds if for all $h\in H$:

[TABLE]

Using Eq. 58, the structure of $\Gamma$ can be determined given the structure of $K$ and $\Sigma$ used in the calculation. Characters of irreducible representations are available for small order symmetry groups, and scaling up to larger order is possible using computational group theory software. Below is an example for the subgroups of a dihedral group $D_{3}$.

Example IV.2

Coupled Duffing oscillators: $Z_{2}\otimes Z_{2,3}$-equivariant system with $Z_{2}\otimes D_{3}$ basis functions.

We consider different coupling schemes of networks of 3 Duffing oscillators shown in Fig. 1(b) and 1(c). We first note that the $Z_{2}$ symmetry generated by a sign flip is still present in the system for both cases, so two non-interacting blocks corresponding to irreducible representations of that group with respect to that action are still present. Now we focus on the structure of $K$ within each of these non-interacting blocks.

First, let the function dictionary symmetry be $\Sigma=D_{3}=\langle r,\kappa|r^{3}=\kappa^{2}=e,\kappa r\kappa=r^{-1}\rangle$ and the true symmetry of the system be $\Gamma=Z_{3}=\langle r|r^{3}=e\rangle$, where $r_{D_{3}}$ and $r_{Z_{3}}$ have the same action. The isotypic component decomposition of $D_{3}$ is defined in Example IV.1 and can be written as $\mathcal{F}=\mathcal{F}_{tr,D_{3}}\oplus\mathcal{F}_{sign,D_{3}}\oplus\mathcal{F}_{st,D_{3}}$. The isotypic component decomposition of $Z_{3}$ is defined as $\mathcal{F}=\mathcal{F}_{tr,Z_{3}}\oplus\mathcal{F}_{\omega,Z_{3}}\oplus\mathcal{F}_{\omega^{2},Z_{3}}$ ($Z_{3}$ has 3 1-dimensional irreducible representations with $\chi_{tr}(r)=1$, $\chi_{\omega}(r)=\omega$, $\chi_{\omega^{2}}(r)=\omega^{2}$). We note that:

• •

$\mathcal{F}_{tr,Z_{3}}\cap\mathcal{F}_{tr,D_{3}}\neq\varnothing$**

• •

$\mathcal{F}_{tr,Z_{3}}\cap\mathcal{F}_{sign,D_{3}}\neq\varnothing$**

• •

$\mathcal{F}_{\omega,Z_{3}}\cup\mathcal{F}_{\omega^{2},Z_{3}}=\mathcal{F}_{st,D_{3}}$**

Thus, the off block-diagonal structure of $K$ is defined by:

• •

$K\mathcal{F}_{tr,D_{3}}\cap\mathcal{F}_{sign,D_{3}}\neq\varnothing$**

• •

$K\mathcal{F}_{sign,D_{3}}\cap\mathcal{F}_{tr,D_{3}}\neq\varnothing$**

• •

other off-diagonal blocks are zeros

This structure is illustrated on Fig. 3(b) and differs from that on Fig. 3(a).

Now, let $\Sigma=D_{3}$ and $\Gamma=Z_{2}$. Here, $Z_{2}=\langle e,\kappa|\kappa^{2}=e\rangle$, and $\kappa_{D_{3}}$ and $\kappa_{Z_{2}}$ have the same action. The isotypic component decomposition of $Z_{2}$ is defined as $\mathcal{F}=\mathcal{F}_{tr,Z_{2}}\oplus\mathcal{F}_{sign,Z_{2}}$.

We note that:

• •

$\mathcal{F}_{tr,Z_{2}}\cap\mathcal{F}_{tr,D_{3}}\neq\varnothing$**

• •

$\mathcal{F}_{tr,Z_{2}}\cap\mathcal{F}_{st,D_{3}}\neq\varnothing$**

• •

$\mathcal{F}_{sign,Z_{2}}\cap\mathcal{F}_{sign,D_{3}}\neq\varnothing$**

• •

$\mathcal{F}_{sign,Z_{2}}\cap\mathcal{F}_{st,D_{3}}\neq\varnothing$**

Additionally:

• •

$\mathcal{F}_{tr,Z_{2}}\cap\mathcal{F}_{sign,D_{3}}=\varnothing$**

• •

$\mathcal{F}_{sign,Z_{2}}\cap\mathcal{F}_{tr,D_{3}}=\varnothing$**

Thus, the off block-diagonal structure of $K$ is defined by:

• •

$K\mathcal{F}_{tr,D_{3}}\cap\mathcal{F}_{sign,D_{3}}=\varnothing$**

• •

other off-diagonal blocks corresponding to interactions between node permutation isotypic components are generally nonzero

This structure is illustrated on Fig. 3(c) and differs from that on Fig. 3(a) and Fig. 3(b).

This example shows that the structure of the approximation $K$ with maximal symmetries assumed provided information about the actual underlying symmetries of the system.

IV.3 Towards realistic systems

In this manuscript, we provide a general approach for dimensionality reduction in the calculation of Koopman operator approximations by exploiting the underlying symmetries of the system’s dynamics and structure. The exact scaling achieved by the reduction depends on the structure of the symmetry group of the dynamical system, specifically, the number of irreducible representations of the symmetry group and their dimensionality.

The results outlined in this manuscript, similar to most of the other literature related to dynamical systems with symmetries, are immediately applicable in the case of existence of exact symmetries in nonlinear dynamics. That is the case when the system is completely deterministic and the initial conditions respect the symmetries of the system. If the symmetries of the system are known and the available trajectories are deterministic, it is always possible to reconstruct the trajectories that are related via the symmetry group of the system. Then, a full set of trajectories respecting the symmetries of the system can be used to approximate the Koopman operator and its eigendecomposition.

However, in many systems that information is not necessarily available ahead of time and the symmetries are not present in data, even if the initial conditions are symmetric, because of the presence of noise in the system. Some of the examples of not fully symmetric data include the following cases and their combinations:

• •

Deterministic systems with measurement noise. DMD for systems with measurement noise and possible ways to correct for it are presented in Ref. dawson2016characterizing . It is shown in Appendix E that in this case the expected values of off-diagonal elements of $K$ computed using the EDMD are zero, so the block decomposition may still be applicable.

• •

Stochastic systems with symmetric initial conditions and process noise. DMD applied to the systems with process noise is studied, for instance, in Ref. bagheri2014effects .

• •

Systems with imperfect symmetries due to sampling and unknown underlying symmetries.

• •

Systems with imperfect symmetries in dynamics (e.g. slight parameter mismatch).

All these cases require separate treatment, and whether the isotypic component decomposition is still useful in computing the Koopman operator approximation will vary depending on specific characteristics of the data available from the system, such as the strength of noise or the trajectory sampling characteristics.

V Conclusion

In this manuscript, we apply tools from group theory and representation theory to study the structure of the Koopman operator for equivariant dynamical systems. This approach can be applied to systems with permutation symmetries (e.g. networks symmetric under node permutations where the information about the symmetries is contained in the adjacency matrix), systems with intrinsic dynamical symmetries, and systems with both types of symmetries present. We find that the operator itself and its approximations can be block diagonalized using a symmetry basis that respects the isotypic component structure related to the underlying symmetry group and the actions of its elements. For the approximation matrix to be exactly block diagonal, the data must respect the symmetries of the system. That can be readily accomplished if the underlying symmetry is known ahead of time (e.g., the topology of the network is known). Symmetry considerations are applicable to both EDMD and kernel DMD, which means they become useful both in the regime when the number of observables is much larger than that of measurements and vice versa.

Acknowledgements.

This work is supported by the U.S. Army Research Laboratory and the U.S. Army Research Office under contract number W911NF-13-1-0340. The authors thank Jordan Snyder, Mehran Mesbahi, and Afshin Mesbahi for useful discussions.

Appendix A Block diagonalization of isotypic components obtained from $d_{p}$-dimensional irreducible representations

We show that $d$-dimensional irreducible representations of $\Gamma$ yield identical blocks of $\mathcal{K}$ in the isotypic component basis obtained using the unitary irreducible representations of $\Gamma$.

Let the function space be decomposed into isotypic components according to the actions of the symmetry group $\Gamma$ of order $|\Gamma|$, $\gamma_{\rho}\in\Gamma_{\rho}$: $\mathcal{F}=\mathcal{F}_{1}\oplus...\oplus\mathcal{F}_{N}$, where $N$ is the number of irreducible representations of $\Gamma$. Let $\mathcal{F}_{p}$ be one of these isotypic components with a corresponding unitary irreducible representation with elements $R_{p}(\gamma)$ corresponding to $\gamma\in\Gamma$, and let $d_{p}$ be the dimensionality of that representation.

The projection operator is defined as:

[TABLE]

It acts on $f\in\mathcal{F}$ to produce sets of projected functions according to:

[TABLE]

We already know that $\mathcal{K}\xi^{p}_{mn}=h_{p}$, where $h_{p}\in\mathcal{F}_{p}$. The subspace $\mathcal{F}_{p}$ can be decomposed into $d_{p}$ components $\mathcal{F}_{p}=\mathcal{F}_{p,1}\oplus...\oplus\mathcal{F}_{p,d_{p}}$, where $\mathcal{F}_{p,m}=\{g|g=\mathcal{P}_{mn}\circ f,~{}f\in\mathcal{F},n=1,...,d_{p}\}$. This is a well-defined decomposition since $\langle\mathcal{P}^{p}_{mn}f,\mathcal{P}^{p}_{kl}h\rangle=\langle f,\mathcal{P}^{p}_{nm}\mathcal{P}^{p}_{kl}h\rangle=\langle f,\delta_{mk}\mathcal{P}^{p}_{nl}h\rangle$ cornwell1997group can be nonzero only when $m=k$.

We want to show that $\mathcal{K}\xi^{p}_{mn}\in\mathcal{F}_{p,m}$ (also true for any linear operator that commutes with the action of the symmetry group). Since the operator commutes with the actions of the group:

[TABLE]

Here, $\mathcal{K}\circ f\equiv h$.

Let $f_{\Gamma}=\{\gamma\circ f|\gamma\in\Gamma\}$. Any set of linearly independent functions that span $f_{\Gamma}$ can be transformed into a symmetry respecting basis obtained by calculating all the projections $\mathcal{P}^{p}_{mn}\circ f_{\gamma}$, where $f_{\gamma}\in f_{\Gamma}$. That corresponds to a block diagonal form of the Koopman operator.

We’ve already shown that $K$, the approximation of $\mathcal{K}$, also commutes with the actions of the elements of $\Gamma$ for $\Gamma$-equivariant dynamical systems with $\Gamma$-equivariant data. Thus, we can obtain an observable dictionary that block diagonalizes $K$ into $|\Gamma|$ blocks, where each $d_{p}$-dimensional irreducible representation results in $d_{p}$ $d_{p}\times d_{p}$-dimensional blocks.

Additionally, suppose $\mathcal{K}\mathcal{P}^{p}_{mn}\circ f=h$, then $\mathcal{K}\mathcal{P}^{p}_{kn}\circ f=\mathcal{P}^{p}_{km}\mathcal{K}\mathcal{P}^{p}_{mn}\circ f=\mathcal{P}^{p}_{km}\circ h$. This gives us the relation between blocks in $\mathcal{K}$ corresponding to the same irreducible representation $p$. In context of the approximation $K$, it means that we get that $K_{p,i}$ (blocks corresponding to $\psi\in\mathcal{F}_{p,i}$) are equal for all $i$ (for data respecting the symmetries of the system and a proper ordering of basis functions).

Appendix B Commutativity of $K$ and $\gamma_{\rho}$ acting in function space

We show that $G$ and $G^{+}$ (+ denotes the Moore-Penrose pseudoinverse) commute. We note that $G$ is a Hermitian matrix since:

[TABLE]

Thus, $G$ is also normal, i.e. $GG^{*}=G^{*}G$. We show that if $G$ is normal, $GG^{+}=G^{+}G$.

Two of the criteria that define the Moore-Penrose pseudoinverse penrose1955generalized state that $G^{+}=G^{+}GG^{+}$ and $(GG^{+})^{*}=GG^{+}$. It follows that the following relation holds: $G^{+}=G^{+}(GG^{+})^{*}=G^{+}(G^{+})^{*}G^{*}$. Using that relation ben2003generalized and commutativity of + and ∗ operations, we obtain:

[TABLE]

Since the action of $\gamma$ commutes with $A$ and $G$, and since $G$ commutes with $G^{+}$, the action of $\gamma$ commutes with $K=G^{+}A$. which is a Koopman operator approximation.

Appendix C Change of basis and the EDMD approximation

We show that rotating the observable dictionary preserves the symmetries of the reconstructed trajectories.

Suppose we have a basis consisting of dictionary functions $\mathcal{D}_{\psi}$ and a dictionary $\mathcal{D}_{\xi}$ obtained by $\Xi=T\Psi$. Let $\Psi(t)=(\psi_{1}(x(t))~{}~{}\psi_{N}(x(t)))^{T}$ and $\Xi(t)=(\xi_{1}(x(t))~{}~{}\xi_{N}(x(t)))^{T}$. We show that rotating the dictionary function vector does not affect the trajectory reconstruction:

[TABLE]

Next, we show that the state reconstruction preserves the symmetries of the system. Let $P$ be the action of the symmetry group on the basis functions $\Psi$. We aim to show that if $\Psi_{t+1}=K\Psi_{t}$, then $P\Psi_{t+1}=KP\Psi_{t}$. It follows directly from the fact that $K$ and $P$ commute:

[TABLE]

Thus, the trajectories of basis functions reconstructed using the EDMD approximation are $\Gamma$-equivariant, just like the original system. In particular, this is true in case of the evolution of the full state observable.

Appendix D Applicability to kernel methods

Kernel DMD introduced in Ref. williams2014kernel is a variant of approximating the Koopman operator matrix most efficient when the number of measurement points is much smaller than the number of basis functions. Kernel DMD relies on evaluating $\hat{G}$ and $\hat{A}$ using the kernel method. Their elements can be found by indirectly evaluating the inner products in the basis function space: $k(x_{m},y_{n})=\Psi(x_{m})\Psi(y_{n})^{*}$ (e.g., if $k$ is a polynomial kernel, $k(x,y)=(1+xy^{T})^{\alpha}$). We note that $k(\gamma x,\gamma y)=k(x,y)$ due to the properties of inner products.

In kernel DMD, $\hat{G}_{ij}=k(x_{i},x_{j})$ and $\hat{A}_{ij}=k(x_{i},y_{j})$. The eigendecomposition of $\hat{G}=Q\Sigma^{2}Q^{T}$ is then used to find the matrix $\hat{K}$ and use it in computing the eigendecomposition of the Koopman operator approximation matrix $K$:

[TABLE]

Again, we pick a particular order of group elements similarly to Eq. 38:

[TABLE]

We also construct a permutation representation of the group with elements denoted by $P_{\gamma_{i}}$ as defined in Eq. 39.

By Cayley’s theorem, such permutations form a group isomorphic to $\Gamma$. Determining the actions $P_{\gamma_{i}}$ of the group generators is sufficient to find the actions of all the group elements. Let $\bm{P}_{\gamma_{k}}=P_{\gamma_{k}}\otimes I_{n\times n}$. We note that $\left(\bm{P}_{\gamma_{k}}\right)^{*}=(\bm{P}_{\gamma_{k}})^{-1}$. It can be shown that:

[TABLE]

We do so for $\hat{A}$, and the proof for $\hat{G}$ is equivalent. We find that:

[TABLE]

And finally, $k(x_{p},y_{l})=k(\gamma_{i}x_{k},\gamma_{i}y_{q})=k(x_{k},y_{q})$.

Since the relation 78 holds, the same reasoning can be applied to block diagonalize the matrix $\hat{K}$. It is sufficient to apply the projection operator stiefel2012group :

[TABLE]

This projection operator is analogous to the one introduced in equation 32, except the symmetry group in this case acts by permuting the group elements.

We can apply the singular value decomposition of $\mathcal{P}$ to obtain the basis for the projection subspaces of irreducible representations (isotypic components). We form the transformation matrix $T$ by finding the singular value decomposition (SVD) and stacking its eigenvectors as rows of $T$ such that $\bm{T}=T\otimes I_{n\times n}$.

Similarly to EDMD, the isotypic component basis simplifies calculating the approximations of $\hat{K}$.

[TABLE]

The modification is summarized in table 2.

Finally, the approximations of Koopman eigenvalues, eigenfunctions, and eigenmodes can be calculated using $K_{D}$, as shown in Ref. williams2014kernel .

Appendix E Deterministic systems with sensor noise

Transfer operators with process and measurement noise were also studied in Ref. sinha2018robust . Characterizing and correcting for the effect of sensor noise in DMD is discussed in Ref. dawson2016characterizing . We need to extend the results to EDMD to quantify the effect of sensor noise on the structure of the matrix $K$. The main modification that needs to be made is the consideration of the effect of the noise in measuring $X$ and $Y$ on $\Psi_{x}$ and $\Psi_{y}$.

Let $X$ and $Y$ be matrices analogous to $\Psi_{x}$ and $\Psi_{y}$ corresponding to the full-state observable evaluated at discrete time steps. We denote the sensor noise matrices by $N_{x}$ and $N_{y}$, so that the measured $X_{n}$ and $Y_{n}$ can be found from $X_{n}=X+N_{x}$ and $Y_{n}=Y+N_{y}$. We assume that the noise distributions respect the symmetries of the system, which might be the case, for instance, for symmetric networks. Moreover, we assume that the noise is state-independent.

We can form vectors $\Psi_{xn}$ and $\Psi_{yn}$ that can be used to find the approximation $K$ using EDMD:

[TABLE]

Here, $(\Psi_{xn})_{ij}=\psi_{j}((X_{n})_{i})$, $(\Psi_{yn})_{ij}=\psi_{j}((Y_{n})_{i})$, and $N_{\Psi,x}$ and $N_{\Psi,y}$ correspond to noise matrices obtained as:

[TABLE]

We aim to show that $\mathbb{E}(PK_{n})=\mathbb{E}(K_{n}P)$, meaning that the expected value of the Koopman operator $K_{n}$ commutes with the permutation matrix corresponding to an element of the symmetry group. If that is the case, then the expected values of the off-block-diagonal elements of $K_{n}$ in a symmetry adapted basis as defined in Eq. 34 are zero. To do that, we can express $K_{n}$ as:

[TABLE]

If the inverse of the first term exists, it can be expanded into the Taylor series with terms of the form below in a weak noise limit. We need to show that:

[TABLE]

Here, the matrices $M_{i}$ are selected from $N_{\Psi,x/y}$ and $\Psi_{x/y}$. That follows directly from:

[TABLE]

Thus, the expected values of the off-block-diagonal elements of $K_{n}$ are zero in the isotypic component basis.