Koopman Operator and its Approximations for Systems with Symmetries

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.

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 wellknown 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, a) avsalova@ucdavis.eduthese 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 for-mation, 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 [13][14][15] .
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 26,30,31 .Those symmetries lead to phenomena such as full synchronization, cluster synchronization, and formation of exotic steady states such as chimeras 8,12,27,38 .Moreover, topological symmetries underlying cluster synchronization of coupled identical elements assist in analyzing the stability of these fully synchronous cluster states 16,30 .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 1 .Symmetries are also key in determining network controllability and observability.For example, Refs. 11,34explored the effect of explicit network symmetries for linear time-independent and time-dependent networks.Similarly, Refs. 24,44considered 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. 29 , 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 wellsuited, 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 23,42 .There are several well developed approaches for obtaining its numerical approximations, such as Ulam's method that re-lies on the discretization of physical space to obtain an approximation of the Perron-Frobenius operator 43 .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 18 .
The Koopman operator is an infinite dimensional linear operator that describes the evolution of observables (functions of the state space) 18,20,23 .Its definition and properties in the context of dynamical systems are provided, for instance, in Ref. 7 , 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 2 , power system analysis 41 , and extracting spatio-temporal patterns of neural data 5 .
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 35 .A more recent method called extended DMD (EDMD) introduced in Ref. 45 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. 46 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. 25 , 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 45 .
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. 28 addresses symmetries of the Perron-Frobenius operator in relation to the admissible symmetry properties of attractors.On the other, Ref. 36 links the spatiotemporal symmetries of the Navier-Stokes equation to the spatial and temporal Koopman operator.Additionally, Ref. 6 noted that symmetry considerations play an important role in discovering governing equations.And, Ref. 17 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 A. 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 7 .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: ẋ = g c (x). (1) Here, x ∈ R n , g c : R n → R n .Let Φ(x(t), ∆t) be a flow map mapping the initial condition x(t) to the solution at time t + ∆t.It is defined in the following way: The system can be discretized with a finite time step ∆t step , so that x i+1 = Φ(x i , ∆t step ).We denote the function evolving the dynamics of this discretized system by g: The Koopman operator is a linear infinite dimensional operator that evolves the functions (referred to as observables) of state space variables f : R n → C. The action of the Koopman operator K on an observable function f for discrete time systems is defined as: 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 7,22 , 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 λ and eigenfunctions φ of the Koopman operator K are defined as: (Kφ)(x) = λφ(x). ( Of particular interest are the Koopman modes that can be used in model reduction and coherency estimation 33,40 .The Koopman modes v i of the full state observable are defined by: and are projections of the observable onto the span of the eigenfunctions of the Koopman operator K.
The other candidate for studying dynamical systems using an operator based approach is the Perron-Frobenius operator P defined as follows for deterministic dynamical systems: Here, ρ(x) is a density on state space, and A ⊆ R n is a subset of the state space.The Perron-Frobenius operator is an adjoint to the Koopman operator 23 , so an approximation of one of them provides an approximation of the other 18 .

B. Koopman operator approximation methods
Extended dynamic mode decomposition (EDMD) introduced in Ref. 45 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. 21 , and its relation to the Perron-Frobenius operator approximation methods is discussed in Ref. 18 .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 ∆t: y i = Φ(x i , ∆t).Typically, the set of snapshots contains measurements from different trajectories in state space.We define a dictionary of linearly independent observables D = {ψ 1 , ..., ψ N } and form vectors of observables Ψ x and Ψ y .Here, Ψ x ∈ R M ×N , where N is the number of dictionary functions used in the approximation, and M is the number of data snapshots.The elements of Ψ x are obtained from (Ψ x ) ij = ψ j (x i ).We also use the notation Ψ(x m ) = (ψ 1 (x m ), ..., ψ N (x m )) for the dictionary functions evaluated at a particular point on the trajectory.
A finite dimensional approximation of the Koopman operator K that we denote as K can be obtained using: Here, Ψ + x denotes the pseudoinverse of Ψ x .We focus on the case of the Moore-Penrose pseudoinverse for the rest of the manuscript 32 .
If the number of snapshots is much higher than the dimensionality of the function dictionary (M 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: 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 18,45 .
The eigendecomposition of K provides the Koopman eigenvalues, eigenfunctions, and modes that allow an approximate linear representation of the underlying system dynamics.Let λ j and u j be the j th eigenvalue and eigen-vector of K. Then the corresponding Koopman eigenfunction can be approximated by: Let b i be the vectors defined by g(x) i = Ψb i , where g(x) i = e * i x denotes the elements the full state observable discussed in Ref. 45 , and B = (b 1 ... b n ).The Koopman eigenmodes can then be obtained as: Here, w i denotes the i th left eigenvector of K.
A modification of EDMD named kernel DMD 46 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 N .The method relies on evaluating the kernel function: That allows efficient computation of M × M matrices Ĝ, Â, and K, where M is the number of trajectory time steps.The eigendecomposition of 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 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.

C. Point symmetries
In this section, we define the concepts useful to study the structure of the Koopman operator 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 Γ.We define groups by their presentations in a form S|R , 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 r|r n = 1 .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 Γ.Let X ⊂ R n be a vector space with elements x ∈ X.We denote the actions γ ρ on a vector space X by γ ρ x if the set of these actions Γ ρ is isomorphic to Γ.A shorthand γ ρ x = γx is sometimes used in the literature when the action corresponding to the subscript ρ 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 γ ρ 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 ẋ = g c (x) be a continuous time system of differential equations.Here, x ∈ R n , and g c : R n → R n .The system is Γ-equivariant with respect to the actions of Γ ρ if for all x ∈ X and γ ρ ∈ Γ ρ : 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: We note that if a continuous time system is Γ-equivariant, so is its discretization.Moreover, the set of trajectories of a γ-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 {γ ρ x 0 , γ ρ x 1 , ...γ ρ 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. 16,27,30,38) 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 ∈ F are functions f : X → C as: 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 γ ∈ Γ to the elements of the general linear group (a group of matrices of degree n with the operation of matrix multiplication denoted by GL(n, V )) on a vector space V (in this case, we are interested in V = C n ).The characters of a group representation T i (γ) are defined as 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 Γ we can obtain all of its irreducible matrix representations.We denote their elements mapping γ ∈ Γ to p × p-dimensional matrices as R i (γ), 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 χ i (γ) of the unitary irreducible representations.
A vector space, e.g. the space of square integrable functions F, can be uniquely decomposed into components that transform like the i th irreducible representation of Γ under the actions of Γ ρ .These components are called isotypic components 13 .We denote these components by F i .An isotypic decomposition of the square integrable function space with respect to Γ ρ is then defined as F = i 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: We can rewrite the above equation as a system of differential equations to obtain: Let x = x y , and let the dynamics be denoted by ẋ = g(x).Let r s = −1 0 0 −1 = −I 2×2 be the action on the state space that flips the signs of variables.The actions r s and e s = I 2×2 form a group Γ s isomorphic to Γ = Z 2 = r|r 2 = e .Let γ s ∈ Γ s , then: Thus, the Duffing oscillator system is Z 2 -equivariant with respect to the actions γ 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 tr (r) = 1 and the sign representation defined by R sign (r) = −1.Then the space of square integrable functions F can be decomposed into F = F tr ⊕ F sign , where In this case, the sets of functions F tr and 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.2We 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 η ij .Then for each node i in the network, we have the following dynamics: This general coupling scheme is used to describe many systems in literature 30,38 .
We now consider the case of a 3-node network.Depending on what the coupling terms are, the system may be Γ-equivariant with respect to different symmetry groups that act by permuting node indexes.Some examples are: . The isotypic components are defined by 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 F = i,j 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 Γ-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 K.
Theorem III.1 For a Γ-equivariant dynamical system ẋ = g(x) and an arbitrary function f , the Koopman operator commutes with the actions of the elements of Γ: 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.
This result is similar to Theorem 3.1 in 28 , where it is shown that the action of Perron-Frobenius operator commutes with the action of the symmetry group Γ for Γequivariant systems.
Corollary III.1.1The space of eigenfunctions of the Koopman operator K with eigenvalue λ for a Γequivariant system is Γ-invariant.
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 K.The result quoted below is useful for that purpose.
Theorem III.2 (Theorem 3.5 in Chapter XII of Ref. 15 ).Let Γ be a compact Lie group acting on the vector space V decomposed into isotypic components V = W 1 ... W t .Let A : V → V be a linear mapping commuting with Γ.Then A(W k ) ⊂ W k .This result is applicable to finite symmetry groups.Isotypic components of F with respect to Γ induce block diagonal structure of the matrix representation of the Koopman operator.Since K and Γ commute, K(F k ) ⊂ 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 9 , 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 9,39 .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 D ψ = {ψ i }.Each of those functions can be expanded in the isotypic component basis with at least one nonzero coefficient α p mn : Here, ξ p mn is a basis function in the p th isotypic component of 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 Γ, where ξ p mn transforms as (m, n) th element of p th irreducible representation of Γ 9 .We define a projection operator and form a new function basis consisting of functions {ξ p mn } as outlined below.
The projection operator is defined as: Here, [R p (γ)] mn denotes the element in n th row and m th column of the i th unitary irreducible representation of γ ∈ Γ, and γ ρ is the group action.We can form an orthonormal basis D ξ = {ξ i } using that projection operator as follows: Here, c np = P p nn ψ, P p nn ψ 1/2 , where , 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: Here, χ p (γ) is a character of the p th irreducible representation of Γ.If this formula is used, each irreducible representation of degree d p provides a basis function, and d 2 p − 1 other basis functions can be formed using the Gram-Schmidt orthogonalization process 9,30,39 .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 ×d p sized blocks in that matrix 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 A. 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 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 Γ and the definition of its actions Γ ρ .
First, we establish that the Koopman operator approximation K commutes with the actions γ ρ of Γ if the data used in the calculation respects the symmetry, meaning the set of pairs of data points satisfies the condition: 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 {γ 1 , ..., γ |Γ| } and create the vectors Ψ x (and analogously Ψ y ) according to that ordering: Given the ordering of the group elements, we can also construct the permutation representation of the group such that: By Cayley's theorem, such permutations form a group isomorphic to Γ. Determining the actions P γ k of the group generators is sufficient to find the actions of all group elements.Let P γ k = P γ k ⊗ I n×n .We note that (P γ k ) * = (P γ k ) −1 .It can be shown that: By definition, A = Ψ * x Ψ y .We note that: That is the case for symmetric trajectories satisfying Eq. ( 27): Thus, A and G commute with the action of the symmetry group.
If G is invertible and G commutes with γ ρ , G −1 commutes with γ ρ as well.Then: 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 Γ commute.
Since K commutes with the actions of Γ, KF i ⊂ F i .This shows that K can be block-diagonalized in the same way as K.
Suppose we start from a dictionary of observables.Since that dictionary is not necessarily an isotypic component dictionary corresponding to Γ 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: • Each isotypic component is present after the isotypic component decomposition of the original function basis: ∀m, p ∃ψ ∈ D, s.t.P p mn ψ = 0 (35)   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 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. 1a.In that case, the system has Z 2 ⊗D 3 symmetry.Suppose we want to construct an isotypic component basis from a given function dictionary D. As an example, we use an initial dictionary D ψ 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 ψ(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 γ ρ ∈ Γ ρ .Due to the form of these functions, γ ρ • ψ(c, x) = ψ(γ −1 ρ c, x).Since Γ = Z 2 ⊗D 3 is a direct product of two groups, we can write the projection operator in the following form: We first consider the irreducible representations of D 3 .
We implement the EDMD algorithm to obtain the approximation of 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 α = 1., β = −1., δ = .5,and η = 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. 2a illustrates the Koopman operator approximation matrix K calculated using an initial dictionary D ψ and requires performing matrix operations on the full 120 × 120 matrix.Fig. 2b 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 × 10 and 2, 20 × 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 Γ by r Γ .In that case, instead of performing O((mr Γ ) α ) operations of matrix inversion, multiplication, and eigendecomposition, it is sufficient to perform these operations for each of the r Γ blocks, with operations being O(m α ).Here 2 < α < 3 19 .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 Γ-equivariant systems is used.We also note that each d p dimensional irreducible representation results in d p of d p × d p dimensional blocks in K 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 I summarizes the modified EDMD algorithm for Γ-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 K evaluated using a kernel trick commonly applied in machine learning 46 .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

Standard EDMD
EDMD for Γ-equivariant systems • Pick a dictionary of N observables • Evaluate the observables at data points xi and yi • Evaluate the entries of G, A: N 2 elements • Obtain G + : N × N matrix • Pick a dictionary of N observables • Identify the symmetry Γ of the system, find the irreducible representations of Γ • Change the basis to a Γ-symmetric basis using Eq. ( 26) and Eq. ( 25) : multiplying at most N/|Γ| |Γ| × |Γ| matrices by vectors |Γ| × 1.Let Np be the number of functions obtained from applying a projection operator P p corresponding to p th irreducible representation of Γ (e.g., Np = N/|Γ| for cyclic groups).
• Evaluate the observables at data points xi and yi, add trajectories to reflect the symmetries if necessary • To obtain the blocks Kpq of K (each isotypic component corresponds to dp blocks), for each p: - calculations.The method is provided in Appendix D.

B. Consequences of symmetry assumptions in the basis
Assume the data is symmetric as defined by Eq. ( 27) with respect to the symmetry group Γ.A "perfect" basis is the one respecting the isotypic decomposition of Γ. Suppose the basis functions belong to isotypic components of Σ = Γ.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 Γ-equivariant and Σ ⊂ Γ and the set of actions of Σ is a subset of actions of Γ, the system is also Σ-equivariant.Thus, picking a basis respecting the isotypic decomposition of Σ will have the block diagonal structure corresponding to Σ.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 Γ-equivariant and Γ ⊂ Σ, functions belonging to particular isotypic components of Σ 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 Σ = Σ 0 ⊗ Γ.In this case, every action of Σ 0 commutes with every action of Γ.Each isotypic component of F with respect to Σ can be expressed as F pq = F p Σ0 ∩ (F Γ ) q , where F p Σ0 de-notes the p th isotypic component of Σ 0 .In this case, the off-diagonal blocks corresponding to interactions between isotypic components F p1q1 and F p2q2 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 Σ = Z 2 ⊗ 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 F with respect to the symmetry group Σ by F p Σ .We note that if F p Σ ∩ F q1 Γ = ∅ and F p Σ ∩ F q2 Γ = ∅ where q 1 and q 2 index different isotypic components of Γ, then the off-diagonal blocks of K corresponding to interactions between those components are generally nonzero.The condition for F p Σ ∩ F q Γ = ∅ is equivalent to: where f is an arbitrary function, and P p Σ denotes the projection operator onto the p th isotypic component with respect to the symmetry group Σ.
Let H be the set of left cosets of Γ in Σ (defined as H = Σ/Γ = {σΓ : σ ∈ Σ}, where σΓ = {σγ : γ ∈ Γ} 13 ).Thus, the condition of Eq. ( 39) holds if for all h ∈ H: Using Eq. ( 41), the structure of Γ can be determined given the structure of K and Σ 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 .
We consider different coupling schemes of networks of 3 Duffing oscillators shown in Fig. 1b and 1c.We first note that the Z 2 symmetry generated by a sign flip is still present in the system for both cases, so two noninteracting 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.
We note that: Thus, the off block-diagonal structure of K is defined by: • KF tr,D3 ∩ F sign,D3 = ∅ • other off-diagonal blocks corresponding to interactions between node permutation isotypic components are generally nonzero This structure is illustrated on Fig. 3c and differs from that on Fig. 3a and Fig. 3b.This example shows that the structure of the approximation K with maximal symmetries assumed provided information about the actual underlying symmetries of the system.

C. 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 systems 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 FIG. 3: Structure of K with basis functions belonging to the isotypic components of Σ = D 3 for different underlying symmetries of the system 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. 10 .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. 3 .
• 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.
The projection operator is defined as: It acts on f ∈ F to produce sets of projected functions according to: ξ p mn = P p mn • f. (A2) We already know that Kξ p mn = h p , where h p ∈ F p .The subspace F p can be decomposed into d p components F p = F p,1 ⊕...⊕F p,dp , where F p,m = {g|g = P mn •f, f ∈ F, n = 1, ..., d p }.This is a well-defined decomposition since P p mn f, P p kl h = f, P p nm P p kl h = f, δ mk P p nl h 9 can be nonzero only when m = k.
We want to show that Kξ p mn ∈ 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: Here, K • f ≡ h.Let f Γ = {γ • f |γ ∈ Γ}.Any set of linearly independent functions that span f Γ can be transformed into a symmetry respecting basis obtained by calculating all the projections P p mn • f γ , where f γ ∈ f Γ .That corresponds to a block diagonal form of the Koopman operator.
We've already shown that K, the approximation of K, also commutes with the actions of the elements of Γ for Γ-equivariant dynamical systems with Γ-equivariant data.Thus, we can obtain an observable dictionary that block diagonalizes K into |Γ| blocks, where each d p -dimensional irreducible representation results in d p d p × d p -dimensional blocks.
Additionally, suppose KP p mn • f = h, then KP p kn • f = P p km KP p mn • f = P p km • h.This gives us the relation between blocks in 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 ψ ∈ F p,i ) are equal for all i (for data respecting the symmetries of the system and a proper ordering of basis functions).

Standard kernel DMD
Kernel DMD for Γ-equivariant systems • Pick a dictionary of N observables • Evaluate the kernel functions at data points xi and yi • Evaluate the entries of Ĝ, Â: M 2 elements • Obtain Ĝ+ : M × M matrix • Find the eigendecomposition of K: M × M matrix • Pick a dictionary of N observables • Identify the symmetry Γ of the system, find the irreducible representations of Γ • Change the basis to a Γ-symmetric basis using Eq. ( 29) and Eq.(D5) • Evaluate the observables at data points xi and yi, add trajectories to reflect the symmetries if necessary • To obtain the blocks Kpq of K (each isotypic component corresponds to dp blocks), for each p: -Evaluate the entries of Ĝp1, Âp1: (Mp/dp) Thus, the expected values of the off-block-diagonal elements of K n are zero in the isotypic component basis.

a)
If all coupling strengths η ij are equal, the network has D 3 symmetry.This case is shown on Fig.1a.Let the state of the system be defined by x = (x 1 y 1 x 2 y 2 x 3 y 3 )T .Then, the symmetry group is presented by D 3 = r, κ|r 3 = κ 2 = e, κrκ = r −1 and generated by the actions r p = If the coupling strengths obey the conditions η ij = η ji and η ij = η jk , the network has Z 3 symmetry.This case is shown on Fig.1b.The symmetry group is presented by Z 3 = r|r 3 = e and generated by the action r p defined above.c) If the coupling strengths obey the conditions η 12 = η 21 and no other equalities hold, the network has Z 2 symmetry.This case is shown on Fig. 1c.The symmetry group is presented by Z 2 = κ|κ 2 = e and generated by the action κ 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 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

FIG. 1 :
FIG.1: Possible symmetries of a network of three identical Duffing oscillators depending on the coupling strength between the oscillators.Different coupling strengths are shown in blue and red.Green arrows correspond to permutational symmetries arising from physical coupling, black arrows correspond to the symmetries of nodal dynamics.
FIG. 2: Structure of K for different choices of dictionary functions Evaluate the entries of Gp1, Ap1: (Np/dp) 2 elements -Obtain G + p1 : (Np/dp) × (Np/dp) matrix -Find Kp1 = G + p1 Ap1: (Np/dp) × (Np/dp) matrices -Find the eigendecomposition of Kp1: (Np/dp) × (Np/dp) matrix -The other Kpq blocks equal to Kp1 • K = p dp q=1 Kpq.Its eigenvalues are the eigenvalues of Kp, and its eigenvectors only have Np nonzero elements.Mathematically, eigenvectors v kl of K are of the form (v kl )i = p δ pk v pl .TABLE I: EDMD vs modified EDMD for Γ-equivariant systems.|Γ| is the order of Γ.The irreducible representations of Γ are indexed by p and are d p -dimensional.

2 elements-
Obtain Ĝ+ p1 : (Mp/dp) × (Mp/dp) matrix -Find Kp1 = Ĝ+ p1 Âp1: (Mp/dp) × (Mp/dp) matrices -Find the eigendecomposition of Kp 1 : (Mp/dp) × (Mp/dp) matrix -The other Kpq blocks equal to Kp1 • K = p dp q=1 Kpq.Its eigenvalues are the eigenvalues of Kp, and its eigenvectors only have Mp nonzero elements.Mathematically, eigenvectors v kl of K are of the form (v kl )i = p δ pk v pl .TABLE II: kernel DMD vs modified kernel DMD for Γ-equivariant systems.|Γ| is the order of Γ.The irreducible representations of Γ are indexed by p and are d p -dimensional.Here, M be the number of data points used by the algorithm, and {(x m , y m )} respect the symmetries of the system.Here, the matrices M i are selected from N Ψ,x/y and Ψ x/y .That follows directly from:E((M 1 P γ k ) * (M 2 P γ k )...(M n−1 P γ k ) * (M n P γ k )) = P −1 γ k E(M * 1 M 2 ...M * n−1 M n )P γ k = E(M 1 M * 2 ...M n−1 M * n ) (E6)