Koopman Operator Based Dynamical Similarity Analysis for Data-driven Quantification of Distance between Dynamics

Authors: shunsuke  kamiya
masafumi  oizumi
jun  kitazono
Venue: ICLR_cc_2024_Workshop_Re-Align
Published: 2024-02-09 03:52:56

Abstract: Quantifying distance between two dynamical neural systems is a fundamental problem in neuroscience and machine learning fields. Neural dynamics are known to possess nonlinear features, which makes comparison between systems difficult. Recently, a promising method to quantify distance between dynamics called Dynamic Similarity Analysis (DSA) is proposed (Ostrow et al., 2023), which measures distance between matrices approximating linear operators defined in time delay embedded space and thus takes nonlinearity into accounts. Although being a strong method, DSA is not free from problems, including obscure interpretability, failure to satisfy the triangle inequality among matrices of different dimensions, and long computational time. To address these problems, we propose a modified
novel version of DSA. Our proposed DSA measures distance between approximated Koopman operators, which has better interpretability as a linear operator that drives dynamics in a mapped space. The distance measure adopted in our method satisfies the triangle inequality even between matrices of different dimensions. This distance measure also allows extremely fast computational time. We applied our method and DSA to Lorenz system (Lorenz, 1963) of various parameters, and found that our method revealed clusters with respect to parameters and dynamical properties, while DSA failed to do so. With theoretical underpinnings of Koopman operators and matrix distance, we propose our method as an effective method to quantify distance between dynamics.

Main paper: 
 introduction: 
 Comparing representations of two or more neural systems has become an important problem in neuroscience and machine learning fields. For instance, comparison between brain activity of different subjects (O'Connell & Chun, 2018; Hebart et al., 2020) , or between brain activity and deep neural network outputs (Yamins et al., 2014; Cichy et al., 2016; Horikawa & Kamitani, 2017) , may provide better understandings of common mechanism of neural information processing, or revealing unknown neural working principles. This information processing is essentially reduced as the activity of neurons -they transmit information to each other and change their activities dynamically. As information processing appears as a time series of the dynamical system, capturing dynamical features of neural activity is vital in comparing foundational computations that systems perform. Comparing dynamical systems, or more specifically, quantifying distance between dynamical systems is thus an important question, possibly opening the door to comparing many different neural systems and information processing machines such as deep neural networks (Sexton & Love, 2022; Fu et al., 2023) or large language models (Kawakita et al., 2023; Wang et al., 2023) .One of the important aspects to consider dynamics of information processing systems is nonlinearity. It has been known that neural systems and machine learning models own nonlinear representational features (Rabinovich & Muezzinoglu, 2010) and that they play important roles in performing com-putations (Breakspear, 2017; Freyer et al., 2012) . Thus, when comparing two information processing systems, methods that take nonlinearity into account are favorable rather than utilizing techniques to compare linear dynamical systems, which have a rather long research history in the field of system identification (Hanzon & Marcus, 1982; Afsari & Vidal, 2013) . A possible approach for comparison taking account of nonlinearity is the use of the Koopman operator method (Koopman, 1931; Mezić, 2005; Brunton et al., 2022) , which has become one of the major data-driven strategies to tackle nonlinear dynamics.Recently, one such method, called Dynamic Similarity Analysis (DSA), has been invented (Ostrow et al., 2023) . DSA is a data-driven method that measures the distance between a pair of time series data. Briefly, this method takes two steps: the first step is Hankel Alternative View of Koopman (HAVOK) analysis (Brunton et al., 2017) , where one finds a matrix for each dynamics as an approximation of a linear operator in time delay coordinate space. The next step is the measuring step, and one quantifies a distance between the obtained pair of matrices. This quantification makes use of what is called Procrustes analysis of vector fields (PAVF), where a change of basis of coordinate is regarded as invariant. The authors claim that DSA succeeded in distinguishing dynamics that are similar in a geometric sense, and identifying dynamical systems of similar dynamical properties that possess different geometrical features (Ostrow et al., 2023) , both of which are unsuccessful with existing methodologies that utilize geometrical characteristics of state space, such as orthogonal Procrustes analysis.Although DSA is a promising and distinctive method to compare dynamics as mentioned above, it is not free from problems. Firstly, its interpretability is not clear. This lack of interpretability of DSA can be attributed to that of HAVOK, whose theoretical foundations and backgrounds have been studied (Arbabi & Mezić, 2017; Hirsh et al., 2021; Kamb et al., 2020) but are yet to be fully explained. Secondly, the distance measure, PAVF, cannot satisfy the triangle inequality when applied to quantify matrices of different dimensions. This may be problematic, as we often apply based on the dissimilarity various machine learning techniques that are based on metric properties of the dissimilarity measure (Williams et al., 2021) . Thirdly, to find distance using PAVF, one has to use optimization algorithms including neural networks (Ostrow et al., 2023) or iterative optimization techniques over Stiefel manifolds (Edelman et al., 1998; Sato & Iwai, 2013) . This can cause a long computational time, limiting the applicability of the method.To address these problems, here we propose an altered version of DSA to quantify distance between two dynamics (Fig. 1 ). We make two alterations. Firstly, our method compares two matrices that are direct approximations of the Koopman operator (Koopman, 1931; Mezić, 2005; Bevanda et al., 2021; Brunton et al., 2022) . Having a clear meaning as linear operators that dominate the dynamics in nonlinearly transformed spaces, the Koopman operator possesses better interpretability and clarity in quantifying distance between dynamics. We can also employ kernel methods to enable flexible representation for the dynamics and at the same time simpler calculations (Baddoo et al., 2022) .Secondly, we modified the distance metric originally proposed in DSA and defined a distance measure called Modified Procrustes analysis of vector fields (M-PAVF). This metric respects the triangle inequality even between matrices of different sizes and successfully evades the problem in the original metric. This modified metric is also advantageous in its computational simplicity, as one can find the metric massively faster than the metric in DSA, preserving similar tendency as the DSA metric.As an application, we employed our version of DSA and the original DSA to time series data of Lorenz system (Lorenz, 1963) of different parameters to see how these methods quantify distances between dynamics. The results show that our method showed clusters of dynamical data of the same parameter, while the original DSA could not. What is more, only did our method succeed in detecting the closeness with respect to dynamics' dynamical properties of different parameters, revealing that our version of DSA can be a favorable measure to quantify distance between dynamics. 
backgrounds: 
 Before describing our method, we first explain some backgrounds about data-driven approaches of nonlinear dynamical analysis. Our study makes use of concepts of the Koopman operator and related methodologies including Extended Dynamic Mode Decomposition (EDMD), kernel EDMD HAVOK Time delay embeddingx 1 x 2 M M d(A 1 , A 2 ) = inf U ∈O(r) A 1 -UA 2 t U F y 1 y 2 x k+1 = F (x k ) xk+1 = A 1 xk y k+1 = G(y k ) K1 EDMD/KEDMD & coordinate transformation K2 Koopman Lifting x 1 x 2 g h M M R m1 R m2 d( K1 , K2 ) = inf U,V ∈St(m2,m1) K1 -U K2 t V F y 1 y 2 g(x 2 ) g(x 1 ) x k+1 = F (x k ) g(x k+1 ) = K 1 g(x k ) h(x k+1 ) = K 2 h(x k ) y k+1 = G(y k ) h(y 2 ) h(y 1 )Original DSA (Ostrow et al. 2023 ) Proposed version of DSAx i x i+1 x i+2 x1 x2 y j y j+1 y j+2 ỹ1 ỹ2 ỹk+1 = A 2 ỹk A 1 r A 2 r m 1 r r m 1 m 2 m 2Procrustes analysis of vector fields Modified Procrustes analysis of vector fields (KEDMD), and Hankel Altered View of Koopman analysis (HAVOK) . In what follows, we are to explain each of these and the original DSA briefly.Initially proposed a century ago in Koopman (1931) and after rediscovery in Mezić (2005) , the Koopman operator has become one of the mainstream approaches in the analysis of nonlinear dynamical systems (Mezic, 2013; Bevanda et al., 2021; Brunton et al., 2022) . In a nutshell, the Koopman operator is a linear operator that drives nonlinearly mapped dynamics. In more detail, assume a discrete-time dynamical systemEQUATIONwith F : M → M. Here F may be a nonlinear transform, and thus the dynamics may be difficult to analyze any further. To get a better perspective, consider a dynamics that is given as a transformation of z n through a function g : M → R,g(z n+1 ) = g • F (z n ), n = 1, 2, • • • . (2)The composition of the functions on the right-hand side can be regarded as an operator K := g • F . This K is called the Koopman operator. When we set g as a member of a Hilbert space F with an inner product ⟨•, •⟩, and assume g • F ∈ F for any g ∈ F, we readily see that the Koopman operator is a linear operator on F.The merit of the Koopman operator is that it transforms a nonlinear dynamics defined as Eq. ( 1) into a linear dynamics. For g ∈ F and {g(zEQUATIONwhich shows that K has a clear mathematical meaning as a linear operator that drives the mapped dynamics {g(z n )}. This mapping g is often called an observable (Williams et al., 2015) . Now with a linear equation Eq.(3) at hand, one may expect to examine properties of the dynamical system in a simple way. For example, one may want to find eigenvalues and eigenfunctions of K to gain an insight into the dynamics, or make a prediction of z n at some unseen time point. Unfortunately, it is usually difficult to do so using Eq. (3), as K is a linear operator in an (generally) infinite dimensional space F. Instead, one has to think of a finite dimensional approximation of K, i.e., one needs to limit the functional space to a finite dimensional subspace. This technique is called Extended Dynamic Mode Decomposition (EDMD) (Williams et al., 2015) and is discussed next.EDMD (Williams et al., 2015) is a data-driven technique to approximate the action of the Koopman operator. As mentioned above, it is often difficult to find the eigenvalues and eigenfunctions of the linear operator K, and one needs to confine the action of the operator in a finite dimensional space.That is, we take an N dimensional space F N (N ∈ N) and consider approximating K with a linear operator K ′ : F N → F N . In EDMD, we take F N as a span:F N := span(ψ 1 , • • • , ψ N ), where the functions ψ 1 , • • • , ψ N ∈ F are arbitrary chosen. Each of the functions ψ i , (i = 1, • • • , N) is called a basis function, and the combination of N functions is called a dictionary. Some typical choices of such a dictionary can be monomial functions, or orthogonal polynomial functions such as Hermite polynomials (under F being an L 2 space with respect to some measure on M).To approximate K with K ′ : F N → F N , we take a function h ∈ F N , and think of the image K ′ h as the projection of Kh onto F N . This projection is necessary as Kh is not always on F N and one has to approximate Kh with a point on F N . Some calculation will show that if {ψ i } N i=1 is linearly independent and h is expressed as h =:N i=1 c i ψ i with c 1 , • • • , c N ∈ R, then K ′ h is written as N i=1 c ′ i ψ i , with t (c ′ 1 , • • • , c ′ N ) = K t (c 1 , • • • , c N ).(4) Here, K is a matrix defined asK = t AG -1 , (5) where (A) ij = ⟨Kψ i , ψ j ⟩ , (G) ij = ⟨ψ i , ψ j ⟩ .(6) Now we obtained the approximated form of K on the finite dimensional subspace F N , which can be conveniently described using linear algebra.How can we find A and G using dynamical data? With the time series data{x n } M n=1 ⊂ M (M ∈ N), we can approximate A and G utilizing Kψ i (x n ) = ψ i (x n+1 ) (n = 1, • • • , M -1)and the Ergodic hypothesis. This approximation, however, can be problematic, with a large number of the basis functions (N ) causing enormous memory load and difficulty in taking the matrix inverse (G -1 ). We often face this large N problem particularly when the dynamics is high-dimensional, since the number of the basis functions grows immensely due to combinatorial explosion. Kernel Extended Dynamic Mode Decomposition (KEDMD) (Williams et al., 2016; Klus et al., 2020) is a method proposed to alleviate such a problem, where kernel methods come in and drastically facilitate the matrix computations. We next look at this method.Basic concepts of KEDMD are almost identical to those of EDMD, except that F is taken to be a reproducing kernel Hilbert space (RKHS) H with respect to a positive definite kernel k : M×M → R and the inner product ⟨•, •⟩ defined using k. In short, k is a generalization of inner products, or similarity measure, of a pair of points in M, and the RKHS is a Hilbert space of functions where functions of the form of k(•, x) (∀ x ∈ M) and their linear sums reside and for any et al., 2020; Baddoo et al., 2022; Rosenfeld et al., 2022) . Note that the subscripts go over 1 to M , not 1 to N . Using the properties of the kernel function, A and G in Eq. ( 6) is written asf ∈ H, ⟨f, k(•, x)⟩ = f (x) holds. In KEDMD, each ψ i is set to be ψ i (•) = k(•, x i ), i = 1, • • • , M (KlusEQUATIONThere are two merits in applying kernel methods to EDMD. First, kernel methods allow us to deal with nonlinear regressions in a tractable manner, thanks to the kernel trick, enabling flexible representation of systems (Shawe-Taylor & Cristianini, 2004) . Second, computational loads lessen because one only has to take M (the number of time points) basis functions at maximum. This drastically improves computational complexity if N is impermissively large, especially due to the matrix inversion process G -1 taking O(M 3 ) steps instead of O(N 3 ). It can still be hard with M basis functions in computing K as in Eq. ( 5) if M is not small enough. But there exist many techniques to evade this problem utilizing essentially row rank property of Gram matrices G, such as Nystrom methods (Drineas & Mahoney, 2005; Sun et al., 2015; Zhang & Kwok, 2010) . Recently, in KEDMD literature, a type of Nystrom method was proposed as a part of an algorithm named linear and nonlinear disambiguation optimization (LANDO) (Baddoo et al., 2022) . We adopt this technique in our computation of K.HAVOK analysis was originally proposed as a method to decompose a chaotic system into a linear system and intermittent forcing (Brunton et al., 2017) . HAVOK makes use of a classical technique called time delay embedding (Kim et al., 1999; Pan & Duraisamy, 2020) , where one stacks time series data {x n } M n=1 ⊂ R d to obtain what is called a Hankel matrix:EQUATIONwhereh ∈ {1, • • • , M }.In HAVOK, one tries to extract a linear relationship between consequent two time steps on this time delay embedded coordinate, but does so using what can be thought as the "basis" for this time series. To compare between two time steps,H (h) 1:M -h and H (h) 2:M -h+1are taken, which we hereafter denote by H 1 and H 2 , respectively, for notational simplicity. One applies singular value decomposition (SVD) of rank r (r ⩽ min(hd, Mh)) to these matrices:H 1 ≃ U 1 Σ 1 t V 1 and H 2 ≃ U 2 Σ 2 t V 2 .These V 1 and V 2 represent dominant structures of the time series {x n } and can be thought as bases for it. Then, using these bases, linear regression is performed to find a matrix A such thatt V 2 = A t V 1 .What does the HAVOK procedure mathematically mean? If H (h) 1:M -h+1 is composed over an infinitesimally short time period (i.e., Mh is small enough), A is a matrix that appears in Frenet-Serret equation, which represents time variation of tangent, normal, and binormal vectors, for the dynamics in time delay embedding space (Hirsh et al., 2021) . In the opposite limit of M → ∞, HAVOK converges to discrete Fourier analysis when {x n } is periodic (Bozzo et al., 2010) . For other general cases, although much work has been done for revealing mathematical properties of HAVOK (Arbabi & Mezić, 2017; Kamb et al., 2020) , a clear-cut mathematical interpretation of HAVOK seems to be missing.DSA is a method to quantify distance between two given systems (Ostrow et al., 2023) . It builds on HAVOK and a matrix distance quantification technique called Procrustes analysis of vector fields (PAVF). Given two dynamical data {x n } M1 n=1 ⊂ R d1 and {y n } M2 n=1 ⊂ R d2 (M 1 , M 2 , d 1 , d 2 ∈ N), the following two steps are taken:1. HAVOK analysis Perform HAVOK to each of the data using the first r ∈ N rows of the Hankel matrices, and obtain matrices A x , A y ∈ R r×r from data {x n } M1 n=1 and {y n } M2 n=1 , respectively. 2. Procrustes analysis of vector fields (PAVF) To see how close A x and A y can be as drift matrices even after a change of basis, solve the next optimization problem:EQUATIONto find the distance between dynamics. Here, O(r) is the set of orthogonal matrices of dimension r.In the original paper (Ostrow et al., 2023) , the authors applied DSA to simulated dynamical data and showed that DSA was able to detect the difference between dynamics that are similar in a geometric sense, and identify dynamical systems of similar dynamical properties that possess different geometrical features, both of which were impossible using existing methods.Being a promising approach that takes nonlinearity into account and owns an empirically favorable ability to distinguish dynamical properties, DSA is not free from problems. We here list three of them. First, its interpretability is not clear. This lack of interpretability is caused by the use of HAVOK, whose theoretical foundations remain unexplained as mentioned above. Secondly, the distance measure, PAVF (Eq. ( 9)), does not always satisfy the triangle inequality when applied to matrices of different dimensions. PAVF was originally proposed to quantify distance between matrices of the same dimension as in Eq. ( 9), and under this condition satisfies the triangle inequality, which is considered favorable as a metric quantifying representational dissimilarity (Williams et al., 2021) . As we are to see later, we can easily extend PAVF for matrices of different dimensions, but that is when a problem occurs; the metric breaches the triangle inequality between matrices of different dimensions. This can be problematic, as after obtaining distance, or dissimilarity, of systems, we often apply machine learning techniques including k-nearest neighbors (Hastie et al., 2009) that are theoretically based on metric properties of the dissimilarity measure. Using metrics that can break the triangle inequality may cause inconsistent results in those later analyses, possibly hindering our understandings. Thirdly, to solve PAVF (Eq. ( 9)), one has to use optimization algorithms including neural networks (Ostrow et al., 2023) or iterative optimization techniques over Stiefel manifolds (Edelman et al., 1998; Sato & Iwai, 2013) . This can cause long computational time, especially when matrix dimensions are large, and thus confines applicability. To address these problems, we propose a novel altered version of DSA, which is to be explained in detail in the next section 3. 
proposed modification of dsa: 
 In this subsection, we explain our modified version of DSA to quantify distance between given two systems. In a nutshell, we measure distance between dynamics by calculating the distance between the approximated Koopman operators. Our method thus consists of the following two steps:1. Apply EDMD or KEDMD to find Koopman approximant matrices for each system.2. Perform a modified version of PAVF and compute the distance between the systems.We are to explain each of the steps in detail in the following. See Fig. 1 for conceptual description.Let us assume that we have access to a pair of dynamical data{x n } M1 n=1 ⊂ M and {y n } M2 n=1 ⊂ M, where M ⊂ R d (M 1 , M 2 , d ∈ N).Here, the dimensions of {x n } and {y n } are set to be identical for a theoretical reason, and practically for data of different dimensions one can zero-pad dynamics as necessary.Apply EDMD or KEDMD First, we perform EDMD or KEDMD to estimate the finite dimensional Koopman approximant for each system. We compute K as in Eq. ( 5) using expressions Eq. ( 6) or Eq. ( 7) for each dynamical data {x n } M1 n=1 and {y n } M2 n=1 , and put them K 1 and K 2 , respectively. Applying EDMD or KEDMD for the first step to quantify distance between dynamics has two main merits. The first is that instead of applying linear methods, EDMD or KEDMD can capture the nonlinear properties of the dynamics since we consider dynamics in the nonlinearly mapped space of observables. As nonlinearity is one of the key factors of dynamical systems including neural circuits (Freyer et al., 2012) , our method may be useful for that purpose. The second is that the matrices used for comparison have clear mathematical meaning as an approximation of the Koopman operator, which is the linear operator in nonlinearly transformed dynamics. This provides better interpretability of what are being compared between dynamics.Here, note that these Koopman approximant matrices K 1 and K 2 are of different dimensions. Readers might think that the dimensions should be the same when applying EDMD, since the number of functions in a dictionary should be the same. However, for simpler and more flexible computation, one can reduce the number of functions in a dictionary in a data-driven manner when using either EDMD or KEDMD. For example, one can apply the almost linearly dependent test in EDMD and KEDMD (Engel et al., 2004) , where one discards redundant bases in terms of almost linear dependency. This reduction technique is recently applied in the Koopman operator literature (Baddoo et al., 2022) . Or one can utilize Nystrom methods for curtailing the row & column numbers of the kernel matrices G and A. After these procedures, the dimension on which we project in the functional space (N in the section 2.2) for estimating K 1 can be different from that of K 2 . We putK 1 ∈ R m1×m1 and K 2 ∈ R m2×m2 .Before calculating the distance between the matrices that are approximation of each dynamics, we have to transform the matrices such that the representation is grounded on orthonormal bases. Eq. (4) implies that K represents transformation of coefficients with respect to the basis {ψ i } m i=1 between adjacent time steps (m = m 1 or m 2 ). To compare K 1 and K 2 in a reasonable manner, we must align the representation such that the basis becomes orthonormal. To get an orthonormal basis {u i } m i=1 of span (ψ 1 , • • • , ψ m ), one can employ QR decomposition in the functional space:(ψ 1 • • • ψ m ) = (u 1 • • • u m ) =: (u 1 • • • u m )R, where the upper triangular matrix R is nonsingular. The vectorsEQUATIONimplying that the representation with respect to the orthonormal basis {u i } is given as RKR -1 =: K. In application, R is explicitly calculated as the Cholesky decomposition of G, since t RR = t R(⟨u i , u j ⟩) ij R = (⟨ψ i , ψ j ⟩) ij = G. We compute K for each dynamics {x n } and {y n } and denote it by K1 and K2 , respectively.Modified Procrustes analysis of vector fields (M-PAVF) Using the matrices K1 ∈ R m1×m1 and K2 ∈ R m2×m2 in the previous section, we move on to quantifying the metric between the two matrices. Note that these matrices are generally not of the same size: we assume m 1 ≥ m 2 without loss of generality. A simple and naive extension of PAVF (Eq. ( 9)) would beEQUATIONor, equivalently,EQUATIONwhere St(m 2 , m 1 ) is a Stiefel manifold of size m 1 ×m 2 , i.e., St(m 2 , m 1 ) = {U ∈ R m1×m2 ; t U U = I m2 } and I m2 is the identity matrix of size m 2 .Although d P might look like a felicitous metric between square matrices of different sizes, it comes with a theoretical and a practical shortcomings. The theoretical problem of d P is that it does not always satisfy the triangle inequality d P (A, B) ≤ d P (B, C) + d P (C, A), which disqualifies d P as a metric (or, more accurately, pseudometric) in a mathematical sense. This can be unfavorable in applying diverse machine learning algorithms after obtaining distances (Williams et al., 2021) .The practical demerit is that to obtain d P one has to resort to optimization algorithms such as neural networks (Ostrow et al., 2023) or iterative optimization over Stiefel manifolds (Edelman et al., 1998; Sato & Iwai, 2013) . This can cause a long computational time, especially when the dimensions of the matrices are large. These problems are discussed in the next subsection 3.2 in detail.Instead, we add a small modification to Eq. ( 13), where we adjust K1 and K2 with two Stiefel elements:EQUATION)We call this metric Modified Procrustes analysis of vector fields (M-PAVF). In fact, this optimization problem has an analytical optimal value, and can be written asEQUATIONwhere {σ( K1) i } m1 i=1 and {σ ( K2) j } m2j=1 are the sets of singular values of K1 and K1 , respectively, and W 2 is the L 2 -Wasserstein distance. Note that when m 1 ̸ = m 2 and for example m 1 > m 2 , one needs to zero-pad the smaller set of singular values. Using this analytical expression, one can immediately prove that d MP defines a proper metric on the space of finite dimensional square matrices.In this subsection, we give our experimental results for comparison of the distance defined by PAVF (d P ) and M-PAVF (d MP ).Triangle Inequality Here we give an example that d P cannot satisfy the triangle inequality while d MP always does. Take A = -2 0 1 0 , B = 2, C = -2. Then computing d P (Eq. ( 13))Table 1 : Values of ρ and property of Lorenz system. Properties of systems were based on Sparrow (1982) . Description of periodic orbits such as [1-2-2] 2 represents that the orbit travels around the first equivalent point for one time and the second two times, repeated two times, then comes back to the same orbit. Orbits of these parameters are shown in Appendix Fig. 4 .5 Periodic [1-1-2-2-1-2] 140 Chaos 152 Periodic [1-1-2-2] 165 Periodic [1-1-2-2] 195 Chaos 220 Periodic [1-1] 250 Periodic [1-1]gives d P (A, B) ≈ 2.920, d P (B, C) = 4, d P (C, A) ≈ 0.7265, which breaches the triangle inequality. On the other hand, computing our new metric d MP (Eq. ( 14)) results ind MP (A, B) ≈ 0.2361, d MP (B, C) = 0, d MP (C, A) ≈ 0.2361, which preserves the triangle inequality as expected.As satisfying the triangle inequality may be preferred in applying various machine learning techniques in later analyses (Williams et al., 2021) , d MP may be a better choice in quantifying distance between matrices. This light computational load can be favorable in examining high dimensional systems, such as neuronal circuits or human brain. Line-like (10) Ring-like (20 We next applied our method to simulated dynamics to see if it can capture dynamical property or dynamical similarity. We did this using the Lorenz system (Lorenz, 1963) : dx dt = σ(yx), dy dt = x(ρz)y, dz dt = xyβz. On simulation, we fixed σ = 10 and β = 8/3 and changed ρ to obtain various dynamical properties of the system. We chose 20 values of ρ according to the properties of dynamics. The values of ρ and the dynamical properties that ρ give are on Table 1 . We performed simulation from the time t = 0 to t = 400 seconds with a time interval of 0.001 second from a common initial value (-8, 8, 27) , and discarded data for the first 200 seconds. Then to collect several data for each parameter, we picked eight periods of the length of five seconds, [200 + 25k, 205 + 25k] , k = 0, • • • , 7, for each ρ. In this way we obtained 160 time series data as a whole (= 20 values of ρ × eight time periods), each of which is described as an R 3×5000 matrix.) [1-1] (220, 250) [1-2-2] (100.5) [1-2-2] 2 (99.7) [1-2-2] 4 (99.6) [1-2-2] 8 (99.537) [1-1-2-2-1-2] (126.5) [1-1-2-2] (152,We then applied our method to every pair of these 160 time series data to find dynamical distances. Before application, we normalized the time series by dividing the whole time series by its average norm. For an approximation of the Koopman operators, we applied KEDMD and set the positive definite kernel as the Gaussian kernel of variance 1 (k(x, y) = exp -|x -y| 2 /2 ). We utilized the almost linear dependency test (Engel et al., 2004; Baddoo et al., 2022) to curtail the number of dimensions of Gram matrix G. For comparison, we also performed DSA in the same manner, with h = 100 points of time delay embedding and the SVD truncation dimension r = 120. This gives two 160-by-160 distance matrices, which are shown in Appendix Fig. 5 (a) and (d). These distance matrix show distinct tendencies, which represents our method and DSA are extracting different features from dynamics. Two-dimensional mappings of dynamics obtained by employing classical multidimensional scaling (cMDS) are shown in Fig. 3 . Each color of points represents each dynamical property described in Table 1 and thus contains one or more than one values of ρ.With our method, we can see clusters of points of the same parameters, while DSA fails to. The mapping using our method also shows that points of similar dynamical properties are placed close, as represented in colors. The results imply that our proposed version of DSA may work better than the original version in identifying and distinguishing dynamics of similar dynamical properties.To check whether KEDMD or M-PAVF was effective, we evaluated distance between dynamics in the following ways: (a) PAVF using HAVOK (i.e., original DSA), (b) M-PAVF using HAVOK, (c) PAVF using KEDMD, and (d) M-PAVF using KEDMD (i.e., proposed method). The results are shown in Appendix Fig. 6 . The panels (a) and (b) show that there are no clusters with respect to parameters or dynamic properties, while the panels in (c) and (d) there exist cluster-like structures. These results imply that applying KEDMD was significant in revealing cluster-like structures of dynamics. 
discussion: 
 In this study, we proposed a new method, a modified version of DSA, which is a data-driven distance quantification between dynamics. Our method quantifies as a distance between dynamics the distance of the Koopman approximant matrices defined in Eq. ( 14). The optimization problem defined in Eq. ( 14) boils down to the Wasserstein distance between singular values of the Koopman operator approximants represented with respect to orthonormal coordinates, and thanks to that analytical expression, we can find distance between dynamics at little computational cost.Although our method reveals clusters when dynamics is periodic as in Fig. 3 , it failed in identifying and distinguishing parameters of chaotic dynamics. Although the prediction of chaotic systems is thought to be difficult due to its sensitivity to initial values (Strogatz, 2000) , it is desirable if our method could extract chaotic dynamical properties better and identify dynamics of the same parameter. We are currently under investigation of hyperparameters and kernel choices that can distinguish chaotic systems well.We only applied our method to Lorenz systems in this paper, but our method has broader applicational directions as its target, and we are on our way to pursuing them. This includes classifying neural dynamics between human subjects performing a behavioral task, revealing learning rules between deep neural networks (DNN), or comparing human and DNN neural systems. Through these applications to various neural data, our method might help examining learning rules of the human brain system, not discarding nonlinear properties of dynamics. 1 0 2 0 2 8 4 5 6 0 7 5 9 9 . 4 9 9 . 5 3 7 9 9 . 6 9 9 . 7 1 0 0 . 5 1 0 0 . 9 3 1 1 5 1 2 6 . 5 1 4 0 1 5 2 1 6 5 1 9 5 2 2 0 2 5 0 1 0 2 0 2 8 4 5 6 0 7 5 9 9 . 4 9 9 . 5 3 7 9 9 . 6 9 9 . 7 1 0 0 . 5 1 0 0 . 9 3 1 1 5 1 2 6 . 5 1 4 0 1 5 2 1 6 5 1 9 5 2 2 0 2 5 0 1 0 2 0 2 8 4 5 6 0 7 5 9 9 . 4 9 9 . 5 3 7 9 9 . 6 9 9 . 7 1 0 0 . 5 1 0 0 . 9 3 1 1 5 1 2 6 . 5 1 4 0 1 5 2 1 6 5 1 9 5 2 2 0 2 5 0 1 0 2 0 2 8 4 5 6 0 7 5 9 9 . 4 9 9 . 5 3 7 9 9 . 6 9 9 . 7 1 0 0 . 5 1 0 0 . 9 3 1 1 5 1 2 6 . 5 1 4 0 1 5 2 1 6 5 1 9 5 2 2 0 2 5 0  
 
 Appendix:  



Reviews: review: title: title: Specific methodological improvements for dynamical similarity analysis show promising initial empiric results in interpretability and tractable implementation, but the paper needs clarifying and restructuring
main_review: main_review: Clarity
* Conceptually clear, and in particular limitations of DSA and how these are addressed by the proposed method are well-explained.
* This structure of DSA problems (1 interpretability, 2 triangle inequality, 3 computational time) is not maintained without. The paper would benefit in terms of clarity and could be made shorter where needed were this the case, including stating in methods which parts of the modified DSA are likely to address which problem. Figure 1 in particular is a great addition in this regard, and the paper should follow that structure.
* Paper is poorly written typographically, linguistically and grammatically: multiple errors of grammar, tense inconsistency, missing commas, missing articles, and other issues limit ease-of-read. I do not list these all as there are numerous, but a careful read-through should be done. Notably, these improve later in the paper (background -- but worse again in section 2.5 -- and results sections).
* I would consider reordering the background section to first discuss DSA, HAVOK, and PAVF as the existing approach, and then the techniques that are employed in the paper -- this would assist with the flow and restating the approach for the reader.

Correctness
* I am not very familiar with the content, but it seems correctly defined mathematically and follows a highly logical flow. E.g. Kernel methods extension is clearly motivated and explained.
* More explanation and intuition of PAVF would help (is it like an alignment transformation then distance?).
* 3.2: should provide a proof that triangle inequality is satisfied rather than example (could be in supplement) -- it is mentioned that m-PAVF is a true metric in section 3.1. I actually think the example can be dropped entirely and the statement made alone.
*3.2: computational time -- what are the costs in terms of accuracy, if any? Given the structure using Cholesky decomposition, you may be able to make theoretical estimates. Whether this is shown empirically or theoretically, I think computational time should be discussed after showing the method works adequately, i.e. after the toy example Lorenz system section 3.3 (and again, to follow the structure of the DSA problem list introduced at the start).
* More discussion of differences in Fig 5 and 6 would be appreciated and interesting; consider restructuring as ablation studies.
* Overall, I think there is too much background and too little discussion of the new method to fully evaluate correctness. It could be interesting to follow the analyses in the Ostrow et al. paper directly.

Novelty
* Methods approach is novel and addresses problems in the existing method literature. Clear acknowledgement of what is new, from where it is adopted, and how it is applied.

Interest
* Relevant new modifications to an existing technique for dynamical systems analysis.
strength_weakness: strength_weakness: Rating:
* Significant new method with initial promising results, that meets criteria in terms of correctness and novelty despite some need for review and improvement.

Fit:
* Clear fit; technical but mathematical details well-explained and of interest.
title: title: A modified version of DSA
main_review: main_review: This paper proposes a modified version of Dynamic Similarity Analysis (DSA) to quantify distance between two dynamical systems. Compared to the original DSA, the proposed method has better interpretability, satisfies the triangle inequality, and has fast computational time. 

**Strengths**
- The modifications on the original DSA are novel and well motivated.
- The experiment on Lorenz systems demonstrates that the proposed method captures dynamical properties more effectively than previous methods.
- The method is able to calculate the distance for Koopman approximate matrices of different dimensions. 

**Weaknesses**
- While having the distance satisfy the triangle inequality is mathematically appealing, it is not clear what the benefits are in applications.
- While the experiments on Lorenz systems seem promising, the authors did not explain why they evaluate their methods on this specific system. Is it a good proxy to real world systems?
strength_weakness: strength_weakness: The proposed method is effective and should be interesting to this community.
title: title: Nice work building upon recent results in the field
main_review: main_review: This paper proposes some tweaks to Dynamical Similarity Analysis (DSA), a recent method introduced by Ostrow et al. (2023; NeurIPS). The bulk of the paper is an explanation of recent developments in Koopman operator theory. The exposition is clear to a general audience; this is an important contribution since the literature being reviewed is recent and I know of no existing authoritative review.

The proposed extensions/reformulations of DSA make sense to me and are generally convincing. I would like the authors to include a derivation of equation 15 (i.e. prove that the modified Procrustes metric has a closed form.

The results / empirical demonstration of this new method need to be fleshed out more, but I think the current submission is sufficient for a workshop.
strength_weakness: strength_weakness: Dynamical similarity analysis is an interesting recent development and this paper gives a thoughtful critique of this method. It is clearly relevant to this workshop and clears the bar for quality.