The Clustering Occupation Kernel (CLOCK) Method for Time Series Arising from Dynamical Systems

# 1Introduction
Source: [https://arxiv.org/html/2610.33354](https://arxiv.org/html/2610.33354)
The Clustering Occupation Kernel \(CLOCK\) Method for Time Series Arising from Dynamical Systems

Michael L\. Wells∗

Kamel Lahouel†Victor Rielly†

Bruno Jedynak∗∗Fariborz Maseeh Department of Mathematics and StatisticsPortland State UniversityPortland, OR 97201†Translational Genomics Research InstitutePhoenix, AZ 85004\{mlwells, bjedyna2\}@pdx\.edu,\{klahouel, vrielly\}@tgen\.org

###### Abstract

We introduce the Clustering Occupation Kernel \(CLOCK\) method, a novel representation\-based framework for clustering time series arising from dynamical systems\. CLOCK embeds observed sequences into a reproducing kernel Hilbert space \(RKHS\) via a time\-averaged occupation kernel, enabling trajectory comparison through their RKHS norm\. For systems governed by differential equations, this formulation rigorously approximates the kernel mean embedding of the invariant Sinai\-Ruelle\-Bowen \(SRB\) measure, meaning our dissimilarity metric approximates the maximum mean discrepancy \(MMD\) between attractors\. To accommodate general time series, we introduce time lags via Takens’ Embedding Theorem to reconstruct the phase space\. Crucially, CLOCK natively handles trajectories of unequal lengths and irregularly spaced observations\. We validate the proposed metric using spectral clustering, demonstrating superior performance against benchmark methods on both synthetic datasets and a real\-world intracranial large\-scale electroencephalogram \(EEG\) dataset\. Finally, ablation studies on the Lorenz\-96 system are consistent with the theoretical convergence rates of the occupation kernels and highlight the method’s strong robustness to zero\-mean Gaussian noise\.

## 1Introduction

Time series, sequences of observations with accompanying timestamps\[[22](https://arxiv.org/html/2610.33354#bib.bib30)\], are ubiquitous objects in data science\. The immense volume of available data makes annotating each time series infeasible, so unsupervised methods are needed to cluster data by shared characteristics\. Time series clustering methods fall roughly into three categories: distance\-based, model\-based, and representation\-based\. Distance\-based methods use a custom distance function defined on time series, so standard methods like K\-means\[[19](https://arxiv.org/html/2610.33354#bib.bib25)\]can be adapted\. The DTW barycenter\-averaging method \(DBA\)\[[25](https://arxiv.org/html/2610.33354#bib.bib3)\]using dynamic time warping\[[30](https://arxiv.org/html/2610.33354#bib.bib2)\]\[[4](https://arxiv.org/html/2610.33354#bib.bib27)\]is an example\. Model\-based methods cluster trajectories based on their similarity in the parameters of a model assumed to govern the data\. For example, time\-series\-based Gaussian mixture models fall into this category\[[36](https://arxiv.org/html/2610.33354#bib.bib26)\]\. Representation\-based methods map the time series to latent feature vectors that capture the essential elements of each time series\. Typical clustering methods in Euclidean space may be applied to these representations\. Autoencoders\[[33](https://arxiv.org/html/2610.33354#bib.bib28)\]\[[17](https://arxiv.org/html/2610.33354#bib.bib29)\]and contrastive\-learning\-based methods\[[40](https://arxiv.org/html/2610.33354#bib.bib12)\]are examples\.

In this work, we adopt a representation\-based perspective that views time series, after preprocessing, as solutions to ordinary differential equations \(ODEs\)\. Specifically, we embed each observed trajectoryy:\[0,T\]→My:\[0,T\]\\to Minto a reproducing kernel Hilbert space \(RKHS\) by considering its time\-averaged occupation embedding\[[28](https://arxiv.org/html/2610.33354#bib.bib6),[27](https://arxiv.org/html/2610.33354#bib.bib4)\]:

y↦1T​∫0TK⁡\(⋅,y⁡\(t\)\)​𝑑t,y\\;\\mapsto\\;\\frac\{1\}\{T\}\\int\_\{0\}^\{T\}K\(\\cdot,y\(t\)\)\\,dt,\(1\)whereKKis a positive\-definite kernel defined on the state spaceMM\.

This representation captures the time\-averaged occupation measure of the trajectory in state space and provides a natural way to compare time series through Hilbert space distances\. Its relevance becomes particularly clear when the observed trajectories are generated by dynamical systems\. According to the Palis conjecture\[[23](https://arxiv.org/html/2610.33354#bib.bib7)\], a generic dynamical system defined as the flow of aCrC^\{r\}\(r≥1r\\geq 1\) map on a compact manifold admits \(i\) finitely many attractors, \(ii\) each of which is supported by a Sinai–Ruelle–Bowen \(SRB\) measure, and \(iii\) the union of the basins of these attractors has full Lebesgue measure\. These SRB measures characterize the long\-term statistical behavior of trajectories\.

The main contribution of this work is to \(1\) show that the embedding defined by \([1](https://arxiv.org/html/2610.33354#S1.E1)\) approximates the kernel mean embedding \(KME\)\[[32](https://arxiv.org/html/2610.33354#bib.bib8)\]of these SRB measures, \(2\) leverage this connection to propose the Clustering Occupation Kernel method \(CLOCK\): a principled clustering method that groups trajectories according to their associated SRB measures, and \(3\) evaluate CLOCK against the state of the art\.

We compute dissimilarities between time series by approximating the maximum mean discrepancy \(MMD\) between their estimated SRB measures\. The authors in\[[5](https://arxiv.org/html/2610.33354#bib.bib37)\]took a related approach, embedding time series in high\-dimensional feature space using sliding windows\. Their goal was to provide a kernel\-based distance measure that estimates distances by matching feature distributions that only partially overlap in their support\. However, for time series generated by a dynamical system, their method clusters trajectories based on transient dynamics governed by the initial distribution\. By contrast, our proposed method clusters trajectories according to their long\-term, asymptotic statistics, making the resulting distances invariant to initial conditions as long as the trajectories are long enough\. Ultimately, the choice between these methods depends on whether the application demands sensitivity to short\-term transients or long\-term systemic invariants\.

Whether the time series is generated by an unknown discrete\-time system or a continuous higher\-order ODE, we can reconstruct its multi\-dimensional phase space by preprocessing the observations using time\-delay embedding before using CLOCK\. Provided the underlying states evolve on a smooth, compact manifold, Takens’ theorem\[[35](https://arxiv.org/html/2610.33354#bib.bib23)\]guarantees that these delay\-coordinate trajectories provide a diffeomorphism to the true orbits of the underlying system\.

We demonstrate the proposed metric’s effectiveness through a series of clustering tasks\. Specifically, we apply spectral clustering to similarity matrices derived from the MMD distances between trajectory embeddings\. We present numerical experiments in which the time series are the solutions of ODEs corrupted by noise\. Different clusters correspond to different vector\-field parameters\. In this context, we compare CLOCK with standard benchmark approaches\. We also evaluate CLOCK on a large medical dataset of electroencephalogram measurements with a gold standard\.

## 2Method

### 2\.1Problem setup

We are given a collection of time series\{yi\}i=1n\\\{y\_\{i\}\\\}\_\{i=1\}^\{n\}, each composed of a finite collection of points\{yi​j\}j=0mi\\\{y\_\{ij\}\\\}\_\{j=0\}^\{m\_\{i\}\}and a corresponding collection of timestamps\{ti​j\}j=0mi\\\{t\_\{ij\}\\\}\_\{j=0\}^\{m\_\{i\}\}\. In order to motivate the method, let the observations be noisy snapshots from the solutions to an initial value problem:

yi​j=xi​\(ti​j\)\+ϵi​j,d​xid​t=fc​\(xi\),xi​\(ti​0\)=xi​0,y\_\{ij\}=x\_\{i\}\(t\_\{ij\}\)\+\\epsilon\_\{ij\},\\frac\{dx\_\{i\}\}\{dt\}=f\_\{c\}\(x\_\{i\}\),\\;\\;x\_\{i\}\(t\_\{i0\}\)=x\_\{i0\},\(2\)wherec∈\{1,…,C\},ϵi​j∼N⁡\(0,Σ\)c\\in\\\{1,\\ldots,C\\\},\\epsilon\_\{ij\}\\sim N\(0,\\Sigma\)are i\.i\.d\. and where\{fc\}c∈\{1,…,C\}\\\{f\_\{c\}\\\}\_\{c\\in\\\{1,\\ldots,C\\\}\}is a finite collection of vector fields, typically within the same family but with different parameters\. The goal is to find a cluster assignment\{1,…,n\}→\{1,…,C\}\\\{1,\\ldots,n\\\}\\to\\\{1,\\ldots,C\\\}that maps each time series in the collection\{\(ti​j,yi​j\)\},1≤i≤n,0≤j≤mi\\\{\(t\_\{ij\},y\_\{ij\}\)\\\},1\\leq i\\leq n,0\\leq j\\leq m\_\{i\}to one ofCCclusters without knowledge of the vector fields or the initial conditions\. In this context, one could consider a model\-based approach jointly estimating the cluster assignment, the vector fields, and the initial conditions using an Expectation\-Maximization algorithm\[[6](https://arxiv.org/html/2610.33354#bib.bib31)\]\. We present below a more direct approach: computing a similarity between trajectories that can be used as input to an existing clustering algorithm\. Note that we do not assume that the time series are of equal length or that the observations are regularly spaced\. In fact, the algorithm described below applies to any time\-series clustering problem\. However, it is justified for data generated by a dynamical system or data that can be preprocessed to exhibit the same statistical characteristics as if it were generated by a dynamical system\.

### 2\.2Background

LetMMbe a compact manifold, and letf:M→Mf\\colon M\\to Mbe aCrC^\{r\}map \(r≥1r\\geq 1\)\. ASinai–Ruelle–Bowen \(SRB\) measure forffis anff\-invariant Borel probability measureμ\\muonMMwith several important properties\. We describe the property of interest as follows:111Another property is the existence of Lebesgue\-absolutely\-continuous conditional measures on unstable manifolds\. See\[[1](https://arxiv.org/html/2610.33354#bib.bib35)\]for details\.Letϕ:M→ℝ\\phi\\colon M\\to\\mathbb\{R\}be a continuous function and lety:\[0,∞\)→My\\colon\[0,\\infty\)\\to Mbe a solution to the differential equation defined byff\. There exists a positive\-Lebesgue\-measure set, denotedB⁡\(μ\)B\(\\mu\)and referred to as thebasin of attraction forμ\\mu, such that ify⁡\(0\)∈B⁡\(μ\)y\(0\)\\in B\(\\mu\), then the following holds:

limT→∞1T​∫0Tϕ⁡\(y⁡\(t\)\)​𝑑t=∫Mϕ⁡\(p\)​𝑑μ​\(p\)\.\\lim\_\{T\\to\\infty\}\\frac\{1\}\{T\}\\int\_\{0\}^\{T\}\\phi\(y\(t\)\)dt=\\int\_\{M\}\\phi\(p\)d\\mu\(p\)\.\(3\)The Palis conjecture states that, for a dense setΦ\\Phiin the space of allCrC^\{r\}\(r≥1r\\geq 1\) mapsM→MM\\to Mwith theCrC^\{r\}topology, eachf∈Φf\\in\\Phihas a finite collection of SRB measuresμ1,…,μn\\mu\_\{1\},\\ldots,\\mu\_\{n\}such that∪iB\(μi\)\\cup\_\{i\}B\(\\mu\_\{i\}\)has total Lebesgue measure\. That is,λ\(∪iB\(μi\)\)=λ\(M\)\\lambda\(\\cup\_\{i\}B\(\\mu\_\{i\}\)\)=\\lambda\(M\)whereλ⁡\(⋅\)\\lambda\(\\cdot\)denotes the Lebesgue measure\.

LetHHbe a reproducing kernel Hilbert space \(RKHS\) of functions defined on a topological space𝒳\\mathcal\{X\}\. Letν\\nube a probability measure on𝒳\\mathcal\{X\}\. Then the KME ofν\\nuis given by the image of the mapping:

ν↦∫MK⁡\(⋅,p\)​𝑑ν​\(p\)\.\\nu\\mapsto\\int\_\{M\}K\(\\cdot,p\)d\\nu\(p\)\.\(4\)The MMD between two probability measuresμ1,μ2\\mu\_\{1\},\\mu\_\{2\}defined on a topological space is the following quantity:

MMD​\(μ1,μ2\)=‖μ^1−μ^2‖H,\\text\{\\rm MMD\}\(\\mu\_\{1\},\\mu\_\{2\}\)=\\\|\\hat\{\\mu\}\_\{1\}\-\\hat\{\\mu\}\_\{2\}\\\|\_\{H\},\(5\)whereμ^1,μ^2\\hat\{\\mu\}\_\{1\},\\hat\{\\mu\}\_\{2\}are the KME’s ofμ1,μ2\\mu\_\{1\},\\mu\_\{2\}, respectively\. The MMD quantifies the dissimilarity of the measuresμ1\\mu\_\{1\}andμ2\\mu\_\{2\}\. We say thatKKischaracteristic\(see\[[32](https://arxiv.org/html/2610.33354#bib.bib8)\]\) when the mapping \([4](https://arxiv.org/html/2610.33354#S2.E4)\) is one\-to\-one\. IfKKis characteristic, then the measuresμ1\\mu\_\{1\}andμ2\\mu\_\{2\}are identical if and only ifMMD​\(μ1,μ2\)=0\\text\{\\rm MMD\}\(\\mu\_\{1\},\\mu\_\{2\}\)=0\.

### 2\.3The clustering occupation kernel method \(CLOCK\)

Assume for now that each trajectoryyiy\_\{i\}is observed continuously over the time interval\[ti​0,ti​mi\]\[t\_\{i0\},t\_\{im\_\{i\}\}\]\. Let us then embedyiy\_\{i\}into an RKHS by computing its time\-averaged occupation kernel\. Thus, we select a continuous, positive\-definite kernelKKwith corresponding RKHSHHand compute the mapping:

yi↦1Ti​∫ti​0ti​miK⁡\(⋅,yi​\(t\)\)​𝑑t=:Li∗,y\_\{i\}\\mapsto\\frac\{1\}\{T\_\{i\}\}\\int\_\{t\_\{i0\}\}^\{t\_\{im\_\{i\}\}\}K\(\\cdot,y\_\{i\}\(t\)\)dt=:L\_\{i\}^\{\*\},\(6\)whereTi:=ti​mi−ti​0T\_\{i\}:=t\_\{im\_\{i\}\}\-t\_\{i0\}\. The motivation for this embedding comes from the following result, proved in Appendix[A](https://arxiv.org/html/2610.33354#A1):

###### Theorem[A\.1](https://arxiv.org/html/2610.33354#A1.Thmtheorem1)\.

LetKKbe a continuous, positive\-definite kernel defined on a compact manifoldMM\. Lety:\[0,∞\)→My\\colon\[0,\\infty\)\\to Mbe a solution to a first\-order ODE satisfying the Palis conjecture\. Then we have that the following holds almost surely:

limT→∞1T​∫0TK⁡\(⋅,y⁡\(t\)\)​𝑑t=∫MK⁡\(⋅,p\)​𝑑μ​\(p\),\\lim\_\{T\\to\\infty\}\\frac\{1\}\{T\}\\int\_\{0\}^\{T\}K\(\\cdot,y\(t\)\)dt=\\int\_\{M\}K\(\\cdot,p\)d\\mu\(p\),\(7\)whereμ\\muis the SRB measure of the attractor of the underlying dynamical system for whichy⁡\(0\)∈B⁡\(μ\)y\(0\)\\in B\(\\mu\)\.

Proof sketch:We prove this theorem in two steps: first we prove weak convergence, and then we prove convergence of the norm, from which strong convergence follows\. The convergence of the norm can be shown using the kernel’s Mercer property\.

Now, according to Theorem[A\.1](https://arxiv.org/html/2610.33354#A1.Thmtheorem1), we can approximate

Li∗≈∫MK⁡\(⋅,p\)​d​μi​\(p\)=:μ^i,L\_\{i\}^\{\*\}\\approx\\int\_\{M\}K\(\\cdot,p\)d\\mu\_\{i\}\(p\)=:\\hat\{\\mu\}\_\{i\},\(8\)whereμi\\mu\_\{i\}is the SRB measure of the attractor such thatyi​\(0\)∈B⁡\(μi\)y\_\{i\}\(0\)\\in B\(\\mu\_\{i\}\)\. Note thatμ^i\\hat\{\\mu\}\_\{i\}is the KME ofμi\\mu\_\{i\}\.

We are now ready to propose the dissimilarityd2​\(⋅,⋅\)d^\{2\}\(\\cdot,\\cdot\)between time series\. It is defined by:

d2​\(yi,yj\)=‖Li∗−Lj∗‖H2=1Ti2​∫ti​0ti​mi∫ti​0ti​miK⁡\(yi​\(t\),yi​\(s\)\)​ds​dt\+1Tj2​∫tj​0tj​mj∫tj​0tj​mjK⁡\(yj​\(t\),yj​\(s\)\)​ds​dt−2Ti​Tj​∫ti​0ti​mi∫tj​0tj​mjK⁡\(yi​\(t\),yj​\(s\)\)​ds​dt\.d^\{2\}\(y\_\{i\},y\_\{j\}\)=\\\|L\_\{i\}^\{\*\}\-L\_\{j\}^\{\*\}\\\|^\{2\}\_\{H\}=\\\\ \\frac\{1\}\{T\_\{i\}^\{2\}\}\\int\_\{t\_\{i0\}\}^\{t\_\{im\_\{i\}\}\}\\int\_\{t\_\{i0\}\}^\{t\_\{im\_\{i\}\}\}K\(y\_\{i\}\(t\),y\_\{i\}\(s\)\)dsdt\+\\\\ \\frac\{1\}\{T\_\{j\}^\{2\}\}\\int\_\{t\_\{j0\}\}^\{t\_\{jm\_\{j\}\}\}\\int\_\{t\_\{j0\}\}^\{t\_\{jm\_\{j\}\}\}K\(y\_\{j\}\(t\),y\_\{j\}\(s\)\)dsdt\-\\\\ \\frac\{2\}\{T\_\{i\}T\_\{j\}\}\\int\_\{t\_\{i0\}\}^\{t\_\{im\_\{i\}\}\}\\int\_\{t\_\{j0\}\}^\{t\_\{jm\_\{j\}\}\}K\(y\_\{i\}\(t\),y\_\{j\}\(s\)\)dsdt\.\(9\)In practice, the time series are observed at a finite number of time points, thus we estimate the integrals with trapezoid\-rule quadratures using the samples\{yi​j\}j=0mi\\\{y\_\{ij\}\\\}\_\{j=0\}^\{m\_\{i\}\}and associated time points\{ti​j\}j=0mi\\\{t\_\{ij\}\\\}\_\{j=0\}^\{m\_\{i\}\}, see Appendix[B](https://arxiv.org/html/2610.33354#A2)\. This dissimilarity is an estimate of‖μ^i−μ^j‖H2\\\|\\hat\{\\mu\}\_\{i\}\-\\hat\{\\mu\}\_\{j\}\\\|^\{2\}\_\{H\}, the squared MMD ofμi\\mu\_\{i\}andμj\\mu\_\{j\}\. Next, we convert dissimilarities to similarities in order to perform spectral clustering\. We define the similarityS⁡\(⋅,⋅\)S\(\\cdot,\\cdot\)according to:

S⁡\(yi,yj\):=exp⁡\(−12​σ2​‖Li∗−Lj∗‖H2\),S\(y\_\{i\},y\_\{j\}\):=\\exp\\left\(\-\\frac\{1\}\{2\\sigma^\{2\}\}\\\|L\_\{i\}^\{\*\}\-L\_\{j\}^\{\*\}\\\|\_\{H\}^\{2\}\\right\),\(10\)whereσ:=median​\{‖Li∗−Lj∗‖H\}i≠jn\\sigma:=\\mbox\{\\rm median\}\\\{\\\|L\_\{i\}^\{\*\}\-L\_\{j\}^\{\*\}\\\|\_\{H\}\\\}\_\{i\\neq j\}^\{n\}\[[10](https://arxiv.org/html/2610.33354#bib.bib9)\]\. We obtain an\(n,n\)\(n,n\)matrix of similarities which we use as input for a standard spectral clustering algorithm\. Note that other similarities may be derived from the given dissimilarity\. All the experiments in this paper are presented with the similarity in \([10](https://arxiv.org/html/2610.33354#S2.E10)\)\.

### 2\.4Rate of convergence

We prove the following result in Appendix[A](https://arxiv.org/html/2610.33354#A1):

###### Theorem[A\.2](https://arxiv.org/html/2610.33354#A1.Thmtheorem2)\.

LetM⊂ℝdM\\subset\\mathbb\{R\}^\{d\}be a compact forward\-invariant region for the flowΦt\\Phi\_\{t\}of aC1C^\{1\}vector field, letμ\\mube an invariant probability measure onMM, and letKKbe a continuous positive\-definite kernel with Hölder feature map\. Suppose that time\-averaged occupation kernels started fromX~∼μ\\tilde\{X\}\\sim\\muconverge toμ^\\hat\{\\mu\}at mean\-square rate𝒪⁡\(T−1\)\\mathcal\{O\}\(T^\{\-1\}\)\(e\.g\. under Lemma[A\.2](https://arxiv.org/html/2610.33354#A1.Thmlemma2)\), and that trajectories started fromX∼νX\\sim\\nuare exponentially asymptotic to trajectories on the attractor, with asymptotic phases distributed absolutely continuously with respect toμ\\mu\(Assumption[A\.1](https://arxiv.org/html/2610.33354#A1.Thmassumption1)\)\. Then

‖1T∫0TK\(⋅,Φt\(X\)\)dt−∫K\(⋅,p\)dμ\(p\)‖H=𝒪ℙ\(T−1/2\)\.\\left\\\|\\frac\{1\}\{T\}\\int\_\{0\}^\{T\}K\(\\cdot,\\Phi\_\{t\}\(X\)\)\\,dt\-\\int K\(\\cdot,p\)\\,d\\mu\(p\)\\right\\\|\_\{H\}=\\mathcal\{O\}\_\{\\mathbb\{P\}\}\(T^\{\-1/2\}\)\.\(11\)

We also provide a visualization of this rate in Section[3\.5](https://arxiv.org/html/2610.33354#S3.SS5)\.

### 2\.5Application to general time series

We use Takens’ embedding theorem\[[35](https://arxiv.org/html/2610.33354#bib.bib23)\]to extend our setting to more general time series\. We augment the time\-series dimension by adding time lags\. That is, we use trajectories of the form\(y⁡\(t\),y⁡\(t\+τ\),…,y⁡\(t\+m​τ\)\)∈ℝm\+1\(y\(t\),y\(t\+\\tau\),\\ldots,y\(t\+m\\tau\)\)\\in\\mathbb\{R\}^\{m\+1\}wherem∈ℕm\\in\\mathbb\{N\}andτ\>0\\tau\>0are fixed\. Ifτ\\tauandmmare chosen appropriately, the time\-lagged trajectories will be the images of integral curves under a smooth embedding of a compact manifold intoℝm\+1\\mathbb\{R\}^\{m\+1\}\. Thus, we may approximate the KMEs of their SRB measures\.

The optimal time\-lagτ\\taumay be found by computing the first local minimum of the average mutual information of\(y⁡\(t\),y⁡\(t\+τ\)\)\(y\(t\),y\(t\+\\tau\)\)\[[9](https://arxiv.org/html/2610.33354#bib.bib5)\], while a large enough embedding dimensionmmmay be selected by minimizing the number of false nearest neighbors of the data points\[[12](https://arxiv.org/html/2610.33354#bib.bib36)\]\.

### 2\.6Scaling the method to large datasets

Like distance\-based methods, this method requires computing dissimilarities for all possible pairs of time series in the dataset, so, by default, its complexity scales quadratically with the number of time seriesnn\. However, by approximating an implicit kernel with an explicit one, e\.g\.,pprandom Fourier features\[[26](https://arxiv.org/html/2610.33354#bib.bib39)\]and using truncated \(rank k\) singular value decomposition, the time and space complexities are reduced to, resp\.,O⁡\(n​p​k\)O\(npk\)andO⁡\(k⁡\(n\+p\)\)O\(k\(n\+p\)\)\.

## 3Experiments

We compared the CLOCK method to a collection of standard methods on several synthetic datasets as well as a real, medical dataset\.222For an implementation of these experiments, see:[https://anonymous\.4open\.science/r/RKHS\-CLOCK\-36BD](https://anonymous.4open.science/r/RKHS-CLOCK-36BD)\. In addition, we performed ablation studies to determine the convergence rate of the time\-averaged occupation kernels and the quality of the clustering as a function of trajectory length\.

All experiments on synthetic data were conducted on a laptop with the Intel Core Ultra 9 Processor 185H, 64 GB RAM, and an NVIDIA GPU with 8GB of memory\. The experiment on real data was conducted on a High\-Performance Computing \(HPC\) cluster using a single NVIDIA A40 GPU with 45 GB RAM\.

### 3\.1Synthetic datasets

For the synthetic datasets, we set the number of clusters to be 2 for simplifying the evaluation and generated 30 trajectories per cluster\. Each cluster was generated by the same differential equation but with different parameters\. Trajectories were integrated numerically using a standard fourth\-order Runge\-Kutta method with a fixed step size\. White noise was then added\. The datasets are visualized in Figure[1](https://arxiv.org/html/2610.33354#S3.F1)\.

1. 1\.FitzHugh\-Nagumo \(FHN\) oscillator :\[[8](https://arxiv.org/html/2610.33354#bib.bib18)\]\[[21](https://arxiv.org/html/2610.33354#bib.bib19)\]This is a 2D system given by the equations: d​x1d​t=x1−x133−x2\+R,d​x2d​t=1τ​\(x1\+a−b​x2\),\\frac\{dx\_\{1\}\}\{dt\}=x\_\{1\}\-\\frac\{x\_\{1\}^\{3\}\}\{3\}\-x\_\{2\}\+R,\\\\ \\frac\{dx\_\{2\}\}\{dt\}=\\frac\{1\}\{\\tau\}\(x\_\{1\}\+a\-bx\_\{2\}\),\(12\)wherea,b,R,τa,b,R,\\tauare constants\. For cluster 1, we useda=0\.7,b=0\.8,R=0\.2,τ=12\.5a=0\.7,b=0\.8,R=0\.2,\\tau=12\.5\. For cluster 2, we useda=3\.2,b=1\.7,R=−0\.4,τ=22\.5a=3\.2,b=1\.7,R=\-0\.4,\\tau=22\.5\. We sampled initial conditions uniformly from the box\[−2,2\]2\[\-2,2\]^\{2\}and integrated trajectories over the time interval\[0,20\]\[0,20\]\. We selected 101 regularly spaced points on each trajectory\. We added zero\-mean Gaussian noise with standard deviation0\.10\.1\.
2. 2\.Van der Pol oscillator:\[[37](https://arxiv.org/html/2610.33354#bib.bib21)\]A 2D system with governing equations: d​x1d​t\\displaystyle\\frac\{dx\_\{1\}\}\{dt\}=x2,d​x2d​t=μ⁡\(1−x12\)​x2−x1,\\displaystyle=x\_\{2\},\\;\\;\\frac\{dx\_\{2\}\}\{dt\}=\\mu\(1\-x\_\{1\}^\{2\}\)x\_\{2\}\-x\_\{1\},\(13\)whereμ\\muis a constant\. We selectedμ=3\\mu=3for cluster 1 andμ=5\\mu=5for cluster 2\. Initial conditions were sampled uniformly in the box\[−2,2\]2\[\-2,2\]^\{2\}\. We integrated trajectories over the time interval\[0,80\]\[0,80\]\. We then sampled 101 points at irregular intervals on each trajectory\. We added Gaussian noise with standard deviation0\.50\.5\.
3. 3\.Lorenz63:\[[14](https://arxiv.org/html/2610.33354#bib.bib22)\]This is a 3D dynamical system\. The governing equations are given by: d​x1d​t=σ\(y−x\),d​x2d​t=x\(ρ−z\)−y,d​x3d​t=x​y−β​z,\\frac\{dx\_\{1\}\}\{dt\}=\\sigma\(y\-x\),\\frac\{dx\_\{2\}\}\{dt\}=x\(\\rho\-z\)\-y,\\\\ \\frac\{dx\_\{3\}\}\{dt\}=xy\-\\beta z,\(14\)whereσ,ρ,β\\sigma,\\rho,\\betaare constants\. For cluster 1, we usedσ=10,ρ=28,β=8/3\\sigma=10,\\rho=28,\\beta=8/3\. For cluster 2, we usedσ=8,ρ=33,β=2/3\\sigma=8,\\rho=33,\\beta=2/3\. Initial conditions were sampled uniformly from the box\[−10,10\]3\[\-10,10\]^\{3\}\. We integrated trajectories over the time interval\[0,60\]\[0,60\]and selected 41 points at irregular intervals on each trajectory\. We then added Gaussian noise with standard deviation1\.51\.5
4. 4\.Kuramoto:\[[13](https://arxiv.org/html/2610.33354#bib.bib17)\]We used a 10D Kuramoto model with governing equations: d​xid​t=ωi\+κ10​∑j=1Nsin⁡\(xj−xi\)\\frac\{dx\_\{i\}\}\{dt\}=\\omega\_\{i\}\+\\frac\{\\kappa\}\{10\}\\sum\_\{j=1\}^\{N\}\\sin\(x\_\{j\}\-x\_\{i\}\)\(15\)fori=1,…,10i=1,\\ldots,10where\{ωi\}i=110\\\{\\omega\_\{i\}\\\}\_\{i=1\}^\{10\}andκ\\kappaare constants\. For both clusters, we setω:=\(1,1,…,1,2\)\\omega:=\(1,1,\\ldots,1,\\sqrt\{2\}\)\. For cluster 1, we setκ=0\.05\\kappa=0\.05, and for cluster 2, we setκ=1\\kappa=1\. Initial conditions were sampled from the box\[−1,1\]10\[\-1,1\]^\{10\}\. We integrated trajectories over the time interval\[0,5\]\[0,5\]and selected 51 regularly\-spaced points\. We then added Gaussian noise with standard deviation0\.50\.5\.
5. 5\.Lorenz96:\[[15](https://arxiv.org/html/2610.33354#bib.bib20)\]We used a 40D Lorenz96 system with governing equations: d​xid​t=\(xi\+1−xi−2\)​xi−1−xi\+F,\\frac\{dx\_\{i\}\}\{dt\}=\(x\_\{i\+1\}\-x\_\{i\-2\}\)x\_\{i\-1\}\-x\_\{i\}\+F,\(16\)where we use the conventionx−1:=x39,x0:=x40x\_\{\-1\}:=x\_\{39\},x\_\{0\}:=x\_\{40\},x41:=x1x\_\{41\}:=x\_\{1\}andFFis a constant\. For cluster 1, we usedF=8F=8, and for cluster 2, we usedF=10F=10\. Initial conditions were sampled uniformly from the box\[−10,10\]40\[\-10,10\]^\{40\}\. We integrated trajectories over the time interval\[0,60\]\[0,60\]and selected 41 irregularly\-spaced points on each trajectory\. We then added Gaussian noise with standard deviation1\.51\.5\.

### 3\.2Real dataset

1. 1\.Bonn EEG time series dataset:\[[3](https://arxiv.org/html/2610.33354#bib.bib38)\]This dataset consists of 1D intracranial electroencephalogram \(EEG\) recordings of the human brain used for epileptic seizure detection research\. There were two clusters in the dataset, each with 100 time series\. Cluster 1 consisted of baseline, seizure\-free readings in the epileptogenic zone, while Cluster 2 contained measurements of seizure activity taken from all regions of the brain where such activity was found\. Automatically identifying seizure activity in EEG recordings has the potential to provide insight into brain function, and to aid in diagnosis and seizure prevention\. Each time series was 23\.6 seconds long with 4097 equally\-spaced samples\. For the CLOCK method, we used Takens’ embedding theorem to augment the time series \(see section[2\.5](https://arxiv.org/html/2610.33354#S2.SS5)\)\.We setm=6m=6andτ=13\\tau=13time\-steps for this experiment\.

\(a\)FHN\.![Refer to caption](https://arxiv.org/html/2610.33354v1/vanderpol_two_clusters.png)\(b\)Van der Pol
![Refer to caption](https://arxiv.org/html/2610.33354v1/lorenz63_two_clusters.png)\(c\)Lorenz63
\(d\)Kuramoto![Refer to caption](https://arxiv.org/html/2610.33354v1/lorenz96_two_clusters.png)\(e\)Lorenz96

Figure 1:Shown are the synthetic datasets in the\(x1,x2\)\(x\_\{1\},x\_\{2\}\)\-plane with coloring given by the true cluster labels\. We present the true trajectories along with the noisy samples\. For higher\-dimensional datapoints, we plot only the first two coordinates\.\.

### 3\.3Evaluation metrics

We used the class labels to define several clustering metrics\. See Appendix[D](https://arxiv.org/html/2610.33354#A4)for explicit formulas\.

1. 1\.Adjusted Rand Index \(ARI\):The adjusted Rand index\[[34](https://arxiv.org/html/2610.33354#bib.bib14)\]counts the proportion of data points that are clustered consistently with the true labels after adjusting for random chance\. It takes values in the range of\[−1,1\]\[\-1,1\]with value11being a perfect clustering,00being entirely random, and negative values being worse than random\.
2. 2\.Normalized Mutual Information \(NMI\):The normalized mutual information\[[18](https://arxiv.org/html/2610.33354#bib.bib15)\]is an information\-theoretic metric which quantifies how much the knowledge of the predicted labels reduces the information content of the true labels\. It is normalized to take values in\[0,1\]\[0,1\]with11being perfect correspondence and00indicating that the two sets of labels are independent\.
3. 3\.Purity:Purity\[[2](https://arxiv.org/html/2610.33354#bib.bib16)\]measures the homogeneity of the predicted clusters in terms of the true cluster assignments\. In the setting of two predicted clusters and two true clusters, the purity is just the accuracy of the clustering\.

### 3\.4Comparators

We tested the method against several benchmark methods\. We implemented each method in a manner similar to that of CLOCK\. Thus, each implementation consisted of the following general steps:

1. \(a\)Map the time series to representative feature vectors in latent space \(using the identity mapping by default\)\.
2. \(b\)Compute pairwise distances between the feature vectors using a distance metric \(using Euclidean distance by default\)\.
3. \(c\)Convert distances to similarities using a Gaussian kernel, choosing the medianσd:=median​\{d⁡\(yi,yj\)\}i≠j\\sigma\_\{d\}:=\\mbox\{\\rm median\}\\\{d\(y\_\{i\},y\_\{j\}\)\\\}\_\{i\\neq j\}as the bandwidth\.
4. \(d\)Obtain cluster assignments via spectral clustering on the matrix of similarities\.
5. \(e\)By default, use linear interpolation to handle irregular spacing\.

We tested against the following collection of methods:

1. 1\.Euclidean distance:This baseline method uses ordinary Euclidean distance to compute distances between time series\. If the time seriesyiy\_\{i\}are given as matrices with shapeT×dT\\times d, withTTthe number of timestamps andddthe dimension, we computed2​\(yi,yj\)=‖yi−yj‖ℝT×d2d^\{2\}\(y\_\{i\},y\_\{j\}\)=\\\|y\_\{i\}\-y\_\{j\}\\\|^\{2\}\_\{\\mathbb\{R\}^\{T\\times d\}\}\. Thus, the similarity metricS⁡\(⋅,⋅\)S\(\\cdot,\\cdot\)is given by a standard Gaussian kernel\.
2. 2\.Dynamic Time Warping \(DTW\):We compute the distance between time seriesyiy\_\{i\}andyjy\_\{j\}with the DTW algorithm\[[30](https://arxiv.org/html/2610.33354#bib.bib2)\]\. This algorithm aligns two time series by stretching and compressing segments with a dynamic one\-to\-many mapping before computing the distance between the aligned points\. Time series with similar features that are out of phase along the time axis will be nearby in terms of DTW distance\.
3. 3\.Shape\-Based Distance \(SBD\):SBD is a distance metric for time series used in the KShape algorithm\[[24](https://arxiv.org/html/2610.33354#bib.bib1)\]\. Like DTW, it aligns the time series before computing distance\. However, unlike DTW, it aligns time series by shifting; not stretching and compressing\. The cross correlation statisticcc, which is a similarity measure in the range\(−1,1\)\(\-1,1\), of the aligned,zz\-normalized time series is computed, and the distance is given by1−c1\-c\. This distance is invariant to scaling and shifting\.
4. 4\.Catch22:The Catch22 method\[[16](https://arxiv.org/html/2610.33354#bib.bib10)\]extracts 22 specially\-chosen features from each time series to use as a low\-dimensional representation\. The 22 selected features were chosen from a pool of over 7000 candidate features and found to be optimal across a broad range of tasks \(such as clustering and classification\) using rigorous evaluation criteria\.
5. 5\.Latent ODE:This is a deep learning method\[[29](https://arxiv.org/html/2610.33354#bib.bib11)\]which models the time series as trajectories generated by a learned vector field and initial conditions in a latent, high\-dimensional space\. Since dynamical information from the trajectories is encoded in the latent initial conditions, we used them as feature vectors for the time series\. This method takes irregular spacing into account, and thus we don’t have to interpolate\.
6. 6\.TS2Vec:TS2Vec\[[40](https://arxiv.org/html/2610.33354#bib.bib12)\]is another deep learning method using a contrastive framework to learn representations of time series\. During training, positive pairs of context views are optimized via a contrastive loss\. The context views are aligned across the time dimension, as well as across spatial dimensions\. The contrastive loss is applied hierarchically at various scales to account for low\- and high\-resolution features\. To handle irregular spacing, we formed a uniform grid of time\-points\[tm​i​n,tm​a​x\]\[t\_\{min\},t\_\{max\}\]and inserted samples at times corresponding to observations\. We masked out the other points\.

Latent ODE and TS2Vec have several hyperparameters to tune\. We used an internal clustering metric, the density\-based clustering validation \(DBCV\) score, for validating them\. See Appendix[E](https://arxiv.org/html/2610.33354#A5)for details\. See Table[1](https://arxiv.org/html/2610.33354#S3.T1)for the results of the numerical experiments\.

### 3\.5Ablation studies

We carried out two ablation studies to determine the the rate of convergence of time\-averaged occupation kernels to the associated KMEs of SRB measures, and the quality of the clustering of the CLOCK method as a function of trajectory length\. For both experiments, we used the Gaussian uniscale kernel \(see section[C](https://arxiv.org/html/2610.33354#A3)for details\)\.

The convergence rate of time\-averaged integrals of observables over trajectories to the spatial integral of the observable with respect to the associated SRB measure is guaranteed to be𝒪\(T−1/2\)\\mathcal\{O\}\(T^\{\-1/2\}\)for systems with sufficient regularity \(see\[[39](https://arxiv.org/html/2610.33354#bib.bib33)\]for details\)\. We computed the empirical convergence rate for the 40D Lorenz\-96 system withF=8F=8\. To do so, we integrated a single trajectory over the interval\[0,300\]\[0,300\]with random initial condition in the box\[−10,10\]40\[\-10,10\]^\{40\}\. We computed a time\-averaged occupation kernel for this trajectory, and used it as a proxy for the KME of the SRB measure of the Lorenz\-96 attractor, denotedμ^\\hat\{\\mu\}\. For each of a collection of equally\-spaced increasing time\-points\{Tj\}j=120,Tj<300\\\{T\_\{j\}\\\}\_\{j=1\}^\{20\},T\_\{j\}<300, we integrated 1000 trajectories over the interval\[0,Tj\]\[0,T\_\{j\}\]with random initial conditions drawn from the box\[−10,10\]40\[\-10,10\]^\{40\}and computed‖Li,j∗−μ^‖H\\\|L\_\{i,j\}^\{\*\}\-\\hat\{\\mu\}\\\|\_\{H\}fori=1,…,1000,j=1,…,20i=1,\\ldots,1000,j=1,\\ldots,20, whereLi,j∗L\_\{i,j\}^\{\*\}is the time\-averaged occupation kernel of trajectoryiiof lengthTjT\_\{j\}\. A log\-log plot of the resulting mean error over all trajectories as a function of trajectory length is plotted in Figure[2](https://arxiv.org/html/2610.33354#S3.F2)\(a\)\. Note that after an initial "burn\-in" period, we observe a straight line with slope≈−0\.48\\approx\-0\.48, as expected\.

For the next experiment, we generated two sets of trajectories for the 40D Lorenz\-96 system, each with 30 trajectories\. The first set corresponded to the parameter choice ofF=8F=8, and the second toF=10F=10\. Each trajectory was integrated over the interval\[0,15\]\[0,15\]with initial conditions drawn randomly from the box\[−10,10\]40\[\-10,10\]^\{40\}\. For each of an equally\-spaced increasing sequence of time\-points\{Tj\}j=1100,0≤Tj<6\\\{T\_\{j\}\\\}\_\{j=1\}^\{100\},0\\leq T\_\{j\}<6, we clustered the trajectories truncated to the interval\[0,Tj\]\[0,T\_\{j\}\]with the CLOCK method\. The mean ARI over 40 runs as a function of trajectory length is plotted in Figure[2](https://arxiv.org/html/2610.33354#S3.F2)\(b\)\.

\(a\)\(b\)
Figure 2:In \(a\), the log\-mean\-squared error between the time\-averaged occupation kernels and the empirical KME of the SRB measure of the Lorenz\-96 attractor is plotted as a function of the log of the trajectory lengthTT, along with a 95% confidence interval \(CI\)\. The line with slope−0\.5\-0\.5and intercept−0\.96\-0\.96is plotted for reference\. In \(b\), the ARI of clusterings with the CLOCK method on the Lorenz\-96 system is plotted with a 95% CI as a function ofTT\.Table 1:Results of Numerical Experiments\. The CLOCK method is listed along with its accompanying kernel\. Information on the kernels used can be found in Appendix[C](https://arxiv.org/html/2610.33354#A3)\. The Gaussian Multiscale kernel \(GM\) with one bandwidth per dimension was the default choice\. The Gaussian Uniscale kernel \(GU\) with a single bandwidth for all dimensions was chosen for high\-dimensional data\. The quadratic kernel \(quad\) and Laplace kernel were also used\. Best results for each metric are in bold\.
Van der PolMethodARINMIPurityCLOCK GM0\.68930\.66150\.9167CLOCK quad\-0\.017200\.5Euclidean\-0\.017200\.5DTW00\.03010\.5167SBD111Catch220\.74690\.70840\.9333Latent ODE\-0\.01240\.00330\.5333TS2Vec111

## 4Discussion

As can be seen from Table[1](https://arxiv.org/html/2610.33354#S3.T1), the CLOCK method achieved superior performance on all but one dataset\. The quadratic kernel failed on the Van der Pol dataset most likely due to the degree mismatch between the underlying system and the kernel\. However, no other method did as well as CLOCK overall, with the closest competitor being TS2Vec\. TS2Vec experienced a significant degradation of performance on the Kuramoto dataset\. This may be due to the trajectories in this dataset being locally similar but globally divergent \(see Figure[1\(d\)](https://arxiv.org/html/2610.33354#S3.F1.sf4)\), so that the hierarchical contrastive framework encoded dissimilar trajectories with similar feature vectors\. We speculate that the overall good performance of the CLOCK method can be attributed to its capacity to accommodate the dynamical information contained in the trajectories\. Although Latent ODE also incorporates dynamical information, it yields a single dynamical model for trajectories with multiple dynamics, which may be problematic\. These results highlight the need for a method specifically\-designed to cluster trajectories generated by differential equations\.

In Appendix[F](https://arxiv.org/html/2610.33354#A6), we demonstrate empirically how the CLOCK algorithm maintains good performance on extremely noisy datasets\.

### 4\.1Limitations

While the proposed embedding and clustering framework can be applied to general time series without explicit assumptions, its theoretical motivation relies on the setting where trajectories are generated by dynamical systems\. In practice, there are datasets for which this assumption is not satisfied\. In such cases, it may be necessary to first transform the data into a representation that better reflects an underlying dynamical structure, for instance via delay\-coordinate embeddings such as Takens’ embedding\.

The method assumes that the time series are of sufficient length to capture the nature of the invariant measures of the attractors\. Moreover, there exist trajectories generated by different dynamical systems that exhibit distinct short\-term \(transient\) behaviors but share the same asymptotic behavior, i\.e\., are associated with the same SRB measure\. In such cases, the proposed method may fail to distinguish between them, as it clusters trajectories based on their long\-term statistical properties rather than their transient dynamics\.

### 4\.2Conclusion

We have presented a clustering method using a novel dissimilarity metric for time series\. The metric is given by the squared RKHS norm between time\-averaged occupation kernels\. The Palis conjecture and Theorem[A\.1](https://arxiv.org/html/2610.33354#A1.Thmtheorem1)assert that this dissimilarity metric measures how distinct the associated SRB measures are, and thus is able to quantify the difference in the dynamics between two time series\.

The method has displayed robust performance on a wide variety of synthetic datasets: low and high dimensions, regular and irregular spacing, and all with a significant amount of noise\. This suggests that the RKHS distance between time\-averaged occupation kernels is an effective metric whenever it is suspected that the data was generated by differential equations\. In addition, this method achieved the best performance on a real dataset consisting of EEG measurements\.

### 4\.3Future work

There are a several directions for future research: Spectral clustering may be replaced with K\-Means or community clustering\. Additionally, different kernels may be selected to compute embeddings, and there are many ways to convert dissimilarities into similarities\. A particularly promising direction is the extension of our framework to stochastic dynamical systems, including stochastic differential equations and continuous\-time Markov chains where clustering could be performed based on invariant measures\.

### AI Use Statement

We were assisted by the Gemini 3\.1 Pro LLM model in this research\. Specifically, the LLM assisted in locating the foundational result that this work is based on: the Palis conjecture justifying the use of time\-averaged occupation kernels as feature vectors for time series\. Additionally, the LLM assisted with the design of the numerical experiments by recommending datasets and evaluation metrics to use\. However, all theoretical proofs and numerical experiments were carried out by the authors themselves, and all theoretical results referenced by the LLM were verified in the original source material\. We take responsibility for the final content of this work\.

## References

- \[1\]\(2020\)Nonuniformly hyperbolic attractors\.Springer\.Cited by:[footnote 1](https://arxiv.org/html/2610.33354#footnote1)\.
- \[2\]E\. Amigó, J\. Gonzalo, J\. Artiles, and F\. Verdejo\(2009\)A comparison of extrinsic clustering evaluation metrics based on formal constraints\.Information Retrieval12\(4\),pp\. 461–486\.Cited by:[item 3](https://arxiv.org/html/2610.33354#S3.I3.i3.p1.1)\.
- \[3\]R\. G\. Andrzejak, K\. Lehnertz, F\. Mormann, C\. Rieke, P\. David, and C\. E\. Elger\(2001\)Indications of nonlinear deterministic and finite\-dimensional structures in time series of brain electrical activity: dependence on recording region and brain state\.Physical Review E64\(6\),pp\. 061907\.Cited by:[item 1](https://arxiv.org/html/2610.33354#S3.I2.i1.p1.1.1)\.
- \[4\]D\. J\. Berndt and J\. Clifford\(1994\)Using dynamic time warping to find patterns in time series\.InProceedings of the 3rd International Conference on Knowledge Discovery and Data Mining,pp\. 359–370\.Cited by:[§1](https://arxiv.org/html/2610.33354#S1.p1.1)\.
- \[5\]S\. Das, H\. N\. Mhaskar, and A\. Cloninger\(2021\)Kernel distance measures for time series, random fields and other structured data\.Frontiers in Applied Mathematics and Statistics7,pp\. 787455\.Cited by:[§1](https://arxiv.org/html/2610.33354#S1.p5.1)\.
- \[6\]A\. P\. Dempster, N\. M\. Laird, and D\. B\. Rubin\(1977\)Maximum likelihood from incomplete data via the em algorithm\.Journal of the Royal Statistical Society: Series B \(Methodological\)39\(1\),pp\. 1–22\.Cited by:[§2\.1](https://arxiv.org/html/2610.33354#S2.SS1.p1.2)\.
- \[7\]D\. Dolgopyat\(1998\)On decay of correlations in anosov flows\.Annals of Mathematics147\(2\),pp\. 357–390\.Cited by:[Appendix A](https://arxiv.org/html/2610.33354#A1.p9.1)\.
- \[8\]R\. FitzHugh\(1961\)Impulses and physiological states in theoretical models of nerve membrane\.Biophysical Journal1\(6\),pp\. 445–466\.Cited by:[item 1](https://arxiv.org/html/2610.33354#S3.I1.i1.p1.1.1)\.
- \[9\]A\. M\. Fraser and H\. L\. Swinney\(1986\)Independent coordinates for strange attractors from mutual information\.Physical Review A33\(2\),pp\. 1134\.Cited by:[§2\.5](https://arxiv.org/html/2610.33354#S2.SS5.p2.1)\.
- \[10\]D\. Garreau, W\. Jitkrittum, and M\. Kanagawa\(2017\)Large sample analysis of the median heuristic\.arXiv preprint arXiv:1707\.07269\.Cited by:[§2\.3](https://arxiv.org/html/2610.33354#S2.SS3.p4.3)\.
- \[11\]J\. Guckenheimer\(1975\)Isochrons and phaseless sets\.Journal of Mathematical Biology1\(3\),pp\. 259–273\.External Links:[Document](https://dx.doi.org/10.1007/BF01273747)Cited by:[§A\.1](https://arxiv.org/html/2610.33354#A1.SS1.SSS0.Px3.p1.1)\.
- \[12\]M\. B\. Kennel, R\. Brown, and H\. D\. Abarbanel\(1992\)Determining embedding dimension for phase\-space reconstruction using a geometrical construction\.Physical Review A45\(6\),pp\. 3403\.Cited by:[§2\.5](https://arxiv.org/html/2610.33354#S2.SS5.p2.1)\.
- \[13\]Y\. Kuramoto\(2005\)Self\-entrainment of a population of coupled non\-linear oscillators\.InInternational Symposium on Mathematical Problems in Theoretical Physics: January 23–29, 1975, Kyoto University/ Kyoto, Japan,pp\. 420–422\.Cited by:[item 4](https://arxiv.org/html/2610.33354#S3.I1.i4.p1.1.1)\.
- \[14\]E\. N\. Lorenz\(1963\)Deterministic nonperiodic flow\.Journal of Atmospheric Sciences20\(2\),pp\. 130–141\.Cited by:[item 3](https://arxiv.org/html/2610.33354#S3.I1.i3.p1.1.1)\.
- \[15\]E\. N\. Lorenz\(1996\)Predictability: a problem partly solved\.InProc\. Seminar on Predictability,Vol\.1,pp\. 1–18\.Cited by:[item 5](https://arxiv.org/html/2610.33354#S3.I1.i5.p1.1.1)\.
- \[16\]C\. H\. Lubba, S\. S\. Sethi, P\. Knaute, S\. R\. Schultz, B\. D\. Fulcher, and N\. S\. Jones\(2019\)Catch22: canonical time\-series characteristics: selected through highly comparative time\-series analysis\.Data Mining and Knowledge Discovery33\(6\),pp\. 1821–1852\.Cited by:[item 4](https://arxiv.org/html/2610.33354#S3.I5.i4.p1.1)\.
- \[17\]P\. Malhotra, V\. TV, L\. Vig, P\. Agarwal, and G\. Shroff\(2017\)TimeNet: pre\-trained deep recurrent neural network for time series classification\.arXiv preprint arXiv:1706\.08838\.Cited by:[§1](https://arxiv.org/html/2610.33354#S1.p1.1)\.
- \[18\]A\. F\. McDaid, D\. Greene, and N\. Hurley\(2011\)Normalized mutual information to evaluate overlapping community finding algorithms\.arXiv preprint arXiv:1110\.2515\.Cited by:[item 2](https://arxiv.org/html/2610.33354#S3.I3.i2.p1.1)\.
- \[19\]J\. B\. McQueen\(1967\)Some methods of classification and analysis of multivariate observations\.InProc\. of 5th Berkeley Symposium on Math\. Stat\. and Prob\.,pp\. 281–297\.Cited by:[§1](https://arxiv.org/html/2610.33354#S1.p1.1)\.
- \[20\]D\. Moulavi, P\. A\. Jaskowiak, R\. J\. G\. B\. Campello, A\. Zimek, and J\. Sander\(2014\)Density\-based clustering validation\.InProceedings of the 2014 SIAM International Conference on Data Mining \(SDM\),pp\. 839–847\.External Links:[Document](https://dx.doi.org/10.1137/1.9781611973440.96)Cited by:[item 1](https://arxiv.org/html/2610.33354#A5.I1.i1.p1.1)\.
- \[21\]J\. Nagumo, S\. Arimoto, and S\. Yoshizawa\(2007\)An active pulse transmission line simulating nerve axon\.Proceedings of the IRE50\(10\),pp\. 2061–2070\.Cited by:[item 1](https://arxiv.org/html/2610.33354#S3.I1.i1.p1.1.1)\.
- \[22\]A\. Nielsen\(2019\)Practical time series analysis: prediction with statistics and machine learning\.O’Reilly Media\.Cited by:[§1](https://arxiv.org/html/2610.33354#S1.p1.1)\.
- \[23\]J\. Palis\(2000\)A global view of dynamics and a conjecture on the denseness of finitude of attractors\.Astérisque261\(xiiixiv\),pp\. 335–347\.Cited by:[§1](https://arxiv.org/html/2610.33354#S1.p3.1)\.
- \[24\]J\. Paparrizos and L\. Gravano\(2016\)K\-shape: efficient and accurate clustering of time series\.SIGMOD Rec\.45\(1\),pp\. 69–76\.External Links:ISSN 0163\-5808,[Link](https://doi-org.proxy.lib.pdx.edu/10.1145/2949741.2949758),[Document](https://dx.doi.org/10.1145/2949741.2949758)Cited by:[item 3](https://arxiv.org/html/2610.33354#S3.I5.i3.p1.1)\.
- \[25\]F\. Petitjean, A\. Ketterlin, and P\. Gançarski\(2011\)A global averaging method for dynamic time warping, with applications to clustering\.Pattern Recognition44\(3\),pp\. 678–693\.Cited by:[§1](https://arxiv.org/html/2610.33354#S1.p1.1)\.
- \[26\]A\. Rahimi and B\. Recht\(2007\)Random features for large\-scale kernel machines\.InAdvances in Neural Information Processing Systems,Vol\.20,pp\. 1177–1184\.Cited by:[§2\.6](https://arxiv.org/html/2610.33354#S2.SS6.p1.1)\.
- \[27\]V\. W\. Rielly, K\. Lahouel, E\. Lew, N\. Fisher, V\. G\. Haney, M\. L\. Wells, and B\. M\. Jedynak\(2025\)MOCK: an algorithm for learning nonparametric differential equations via multivariate occupation kernel functions\.Transactions on Machine Learning Research\.Note:External Links:ISSN 2835\-8856,[Link](https://openreview.net/forum?id=fjVIp2Z9RS)Cited by:[Appendix A](https://arxiv.org/html/2610.33354#A1.p3.1.1),[§1](https://arxiv.org/html/2610.33354#S1.p2.1)\.
- \[28\]J\. A\. Rosenfeld, R\. Kamalapurkar, B\. Russo, and T\. T\. Johnson\(2019\)Occupation kernels and densely defined liouville operators for system identification\.In2019 IEEE 58th Conference on Decision and Control \(CDC\),External Links:[Document](https://dx.doi.org/10.1109/cdc40024.2019.9029337),[Link](https://doi.org/10.1109)Cited by:[§1](https://arxiv.org/html/2610.33354#S1.p2.1)\.
- \[29\]Y\. Rubanova, R\. T\. Chen, and D\. K\. Duvenaud\(2019\)Latent ordinary differential equations for irregularly\-sampled time series\.Advances in Neural Information Processing Systems32\.Cited by:[item 5](https://arxiv.org/html/2610.33354#S3.I5.i5.p1.1)\.
- \[30\]H\. Sakoe and S\. Chiba\(1978\)Dynamic programming algorithm optimization for spoken word recognition\.IEEE Transactions on Acoustics, Speech, and Signal Processing26\(1\),pp\. 43–49\.External Links:[Document](https://dx.doi.org/10.1109/TASSP.1978.1163055)Cited by:[§1](https://arxiv.org/html/2610.33354#S1.p1.1),[item 2](https://arxiv.org/html/2610.33354#S3.I5.i2.p1.1)\.
- \[31\]K\. R\. Shahapure and C\. Nicholas\(2020\)Cluster quality analysis using silhouette score\.In2020 IEEE 7th International Conference on Data Science and Advanced Analytics \(DSAA\),pp\. 747–748\.Cited by:[item 1](https://arxiv.org/html/2610.33354#A5.I1.i1.p1.1)\.
- \[32\]B\. K\. Sriperumbudur, K\. Fukumizu, and G\. R\. Lanckriet\(2011\)Universality, characteristic kernels and rkhs embedding of measures\.\.Journal of Machine Learning Research12\(7\)\.Cited by:[§1](https://arxiv.org/html/2610.33354#S1.p4.1),[§2\.2](https://arxiv.org/html/2610.33354#S2.SS2.p2.3)\.
- \[33\]N\. Srivastava, E\. Mansimov, and R\. Salakhudinov\(2015\)Unsupervised learning of video representations using lstms\.InInternational Conference on Machine Learning,pp\. 843–852\.Cited by:[§1](https://arxiv.org/html/2610.33354#S1.p1.1)\.
- \[34\]D\. Steinley\(2004\)Properties of the hubert\-arable adjusted rand index\.\.Psychological Methods9\(3\),pp\. 386\.Cited by:[item 1](https://arxiv.org/html/2610.33354#S3.I3.i1.p1.1)\.
- \[35\]F\. Takens\(2006\)Detecting strange attractors in turbulence\.InDynamical Systems and Turbulence, Warwick 1980: Proceedings of a Symposium held at the University of Warwick 1979/80,pp\. 366–381\.Cited by:[§1](https://arxiv.org/html/2610.33354#S1.p6.1),[§2\.5](https://arxiv.org/html/2610.33354#S2.SS5.p1.1)\.
- \[36\]R\. Umatani, T\. Imai, K\. Kawamoto, and S\. Kunimasa\(2023\)Time series clustering with an em algorithm for mixtures of linear gaussian state space models\.Pattern Recognition138,pp\. 109375\.External Links:ISSN 0031\-3203,[Document](https://dx.doi.org/https%3A//doi.org/10.1016/j.patcog.2023.109375),[Link](https://www.sciencedirect.com/science/article/pii/S0031320323000766)Cited by:[§1](https://arxiv.org/html/2610.33354#S1.p1.1)\.
- \[37\]B\. Van der Pol\(1926\)LXXXVIII\. on “relaxation\-oscillations”\.The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science2\(11\),pp\. 978–992\.Cited by:[item 2](https://arxiv.org/html/2610.33354#S3.I1.i2.p1.1.1)\.
- \[38\]L\. Young\(1998\)Statistical properties of dynamical systems with some hyperbolicity\.Annals of Mathematics147\(3\),pp\. 585–650\.Cited by:[Appendix A](https://arxiv.org/html/2610.33354#A1.p9.1)\.
- \[39\]L\. Young\(1999\)Recurrence times and rates of mixing\.Israel Journal of Mathematics110\(1\),pp\. 153–188\.Cited by:[Appendix A](https://arxiv.org/html/2610.33354#A1.p6.1),[§3\.5](https://arxiv.org/html/2610.33354#S3.SS5.p2.1)\.
- \[40\]Z\. Yue, Y\. Wang, J\. Duan, T\. Yang, C\. Huang, Y\. Tong, and B\. Xu\(2022\)Ts2vec: towards universal representation of time series\.InProceedings of the AAAI Conference on Artificial Intelligence,Vol\.36,pp\. 8980–8987\.Cited by:[§1](https://arxiv.org/html/2610.33354#S1.p1.1),[item 6](https://arxiv.org/html/2610.33354#S3.I5.i6.p1.1)\.

Supplementary Materials

## Appendix ATheoretical results

An SRB measureμ\\mufor an attractor of an ODE defined on a compact manifoldMMhas the property that

limT→∞1T​∫0Tϕ⁡\(y⁡\(t\)\)​𝑑t=∫Mϕ⁡\(p\)​𝑑μ​\(p\)\\lim\_\{T\\to\\infty\}\\frac\{1\}\{T\}\\int\_\{0\}^\{T\}\\phi\(y\(t\)\)dt=\\int\_\{M\}\\phi\(p\)d\\mu\(p\)\(17\)for every continuous functionϕ:M→ℝ\\phi\\colon M\\to\\mathbb\{R\}and trajectoryy:\[0,∞\)→My\\colon\[0,\\infty\)\\to Msuch thaty⁡\(0\)∈B⁡\(μ\)y\(0\)\\in B\(\\mu\), the basin of attraction forμ\\mu\. The Palis conjecture states that, for a dense set in the set of allCrC^\{r\}maps \(r≥1r\\geq 1\) onMMwith theCrC^\{r\}topology, there exists a finite collection of SRB measuresμ1,…,μn\\mu\_\{1\},\\ldots,\\mu\_\{n\}such that∪iB\(μi\)\\cup\_\{i\}B\(\\mu\_\{i\}\)has total Lebesgue measure\. That is,λ\(∪iB\(μi\)\)=λ\(M\)\\lambda\(\\cup\_\{i\}B\(\\mu\_\{i\}\)\)=\\lambda\(M\)whereλ⁡\(⋅\)\\lambda\(\\cdot\)denotes the Lebesgue measure\.

We wish to show that \([17](https://arxiv.org/html/2610.33354#A1.E17)\) extends to feature vectors in an RKHS\. More specifically, we want to show that

limT→∞1T​∫0TK⁡\(⋅,y⁡\(t\)\)​𝑑t=∫MK⁡\(⋅,p\)​𝑑μ​\(p\)\.\\lim\_\{T\\to\\infty\}\\frac\{1\}\{T\}\\int\_\{0\}^\{T\}K\(\\cdot,y\(t\)\)dt=\\int\_\{M\}K\(\\cdot,p\)d\\mu\(p\)\.\(18\)This is the content of Theorem[A\.1](https://arxiv.org/html/2610.33354#A1.Thmtheorem1)below\. In Corollary ???, we show that the rate of convergence is𝒪\(T−1/2\)\\mathcal\{O\}\(T^\{\-1/2\}\)\.

###### Proposition A\.1\.

LetK:M×M→ℝK\\colon M\\times M\\to\\mathbb\{R\}be a continuous kernel defined on a compact manifoldMM\. Let\{Ti\}i=1∞\\\{T\_\{i\}\\\}\_\{i=1\}^\{\\infty\}be an increasing sequence of positive numbers such thatTi→∞T\_\{i\}\\to\\infty\. Lety:\[0,∞\)→My\\colon\[0,\\infty\)\\to Mdenote a trajectory generated by an ODE such thaty⁡\(0\)∈B⁡\(μ\)y\(0\)\\in B\(\\mu\), the basin of attraction of an SRB measureμ\\mu\. Define the linear functionalLi:H→ℝL\_\{i\}\\colon H\\to\\mathbb\{R\}according to

Li​ϕ=1Ti​∫0Tiϕ⁡\(y⁡\(t\)\)​𝑑t\.L\_\{i\}\\phi=\\frac\{1\}\{T\_\{i\}\}\\int\_\{0\}^\{T\_\{i\}\}\\phi\(y\(t\)\)dt\.\(19\)Then eachLiL\_\{i\}is bounded, and their Riesz representers, denotedLi∗L\_\{i\}^\{\*\}, converge weakly to the functionμ^\\hat\{\\mu\}defined by

μ^:=∫MK⁡\(⋅,p\)​𝑑μ​\(p\)\.\\hat\{\\mu\}:=\\int\_\{M\}K\(\\cdot,p\)d\\mu\(p\)\.\(20\)

###### Proof\.

The boundedness of theLiL\_\{i\}follows from the compactness ofMMand the continuity ofKK\(see\[[27](https://arxiv.org/html/2610.33354#bib.bib4)\]for details\)\. For the same reason, the operator

ϕ↦∫Mϕ⁡\(p\)​𝑑μ​\(p\)\\phi\\mapsto\\int\_\{M\}\\phi\(p\)d\\mu\(p\)\(21\)is bounded, and its Riesz representer isμ^\\hat\{\\mu\}\. The functionsLi∗L\_\{i\}^\{\*\}are given by

Li∗=1Ti​∫0TiK⁡\(⋅,y⁡\(t\)\)​𝑑t\.L\_\{i\}^\{\*\}=\\frac\{1\}\{T\_\{i\}\}\\int\_\{0\}^\{T\_\{i\}\}K\(\\cdot,y\(t\)\)dt\.\(22\)BecauseKKis continuous, every functionϕ\\phiin the corresponding RKHS is continuous\. Thus, sinceμ\\muis an SRB measure, we get

⟨ϕ,Li∗⟩H=Li​ϕ=1Ti​∫0Tiϕ⁡\(y⁡\(t\)\)​𝑑t→∫Mϕ⁡\(p\)​𝑑μ​\(p\)=⟨ϕ,μ^⟩H\.\\langle\\phi,L\_\{i\}^\{\*\}\\rangle\_\{H\}=L\_\{i\}\\phi=\\frac\{1\}\{T\_\{i\}\}\\int\_\{0\}^\{T\_\{i\}\}\\phi\(y\(t\)\)dt\\to\\int\_\{M\}\\phi\(p\)d\\mu\(p\)=\\langle\\phi,\\hat\{\\mu\}\\rangle\_\{H\}\.\(23\)∎

###### Lemma A\.1\.

LetLi∗L\_\{i\}^\{\*\}andμ^\\hat\{\\mu\}be as in Proposition[A\.1](https://arxiv.org/html/2610.33354#A1.Thmproposition1)\. Then the following holds:

‖Li∗‖H→‖μ^‖H\\\|L\_\{i\}^\{\*\}\\\|\_\{H\}\\to\\\|\\hat\{\\mu\}\\\|\_\{H\}\(24\)asi→∞i\\to\\infty\.

###### Proof\.

We have

‖Li∗‖H2=⟨Li∗,Li∗⟩H=1Ti2​∫0Ti∫0TiK⁡\(y⁡\(s\),y⁡\(t\)\)​𝑑s​𝑑t\.\\\|L\_\{i\}^\{\*\}\\\|\_\{H\}^\{2\}=\\langle L\_\{i\}^\{\*\},L\_\{i\}^\{\*\}\\rangle\_\{H\}=\\frac\{1\}\{T\_\{i\}^\{2\}\}\\int\_\{0\}^\{T\_\{i\}\}\\int\_\{0\}^\{T\_\{i\}\}K\(y\(s\),y\(t\)\)dsdt\.\(25\)SinceKKis continuous andMMis compact,KKis a Mercer kernel\. Thus,

1Ti2​∫0Ti∫0TiK⁡\(y⁡\(s\),y⁡\(t\)\)​𝑑s​𝑑t=1Ti2​∫0Ti∫0Ti∑j=1∞λj​ϕj​\(y⁡\(s\)\)​ϕj​\(y⁡\(t\)\)​𝑑s​𝑑t\\frac\{1\}\{T\_\{i\}^\{2\}\}\\int\_\{0\}^\{T\_\{i\}\}\\int\_\{0\}^\{T\_\{i\}\}K\(y\(s\),y\(t\)\)dsdt=\\frac\{1\}\{T\_\{i\}^\{2\}\}\\int\_\{0\}^\{T\_\{i\}\}\\int\_\{0\}^\{T\_\{i\}\}\\sum\_\{j=1\}^\{\\infty\}\\lambda\_\{j\}\\phi\_\{j\}\(y\(s\)\)\\phi\_\{j\}\(y\(t\)\)dsdt\(26\)for some collection of non\-negative eigenvalues\{λj\}\\\{\\lambda\_\{j\}\\\}and eigenfunctions\{ϕj\}\\\{\\phi\_\{j\}\\\}\. Let

hi,N:=1Ti2​∫0Ti∫0Ti∑j=1Nλj​ϕj​\(y⁡\(s\)\)​ϕj​\(y⁡\(t\)\)​𝑑s​𝑑t=∑j=1Nλj​\(1Ti​∫0Tiϕj​\(y⁡\(t\)\)​𝑑t\)2\.h\_\{i,N\}:=\\frac\{1\}\{T\_\{i\}^\{2\}\}\\int\_\{0\}^\{T\_\{i\}\}\\int\_\{0\}^\{T\_\{i\}\}\\sum\_\{j=1\}^\{N\}\\lambda\_\{j\}\\phi\_\{j\}\(y\(s\)\)\\phi\_\{j\}\(y\(t\)\)dsdt=\\sum\_\{j=1\}^\{N\}\\lambda\_\{j\}\\left\(\\frac\{1\}\{T\_\{i\}\}\\int\_\{0\}^\{T\_\{i\}\}\\phi\_\{j\}\(y\(t\)\)dt\\right\)^\{2\}\.\(27\)We have that for any fixedNN, the sequencehi,Nh\_\{i,N\}converges asi→∞i\\to\\inftysinceμ\\muis an SRB measure\. Likewise, for any fixedii, the sequence has a limit asN→∞N\\to\\inftysince the sequence of partial sums converges to the time\-averaged double integral of the kernel over two copies of the trajectory \(using the Dominated Convergence Theorem to pass the limit inside the integrals\)\. We have that

limN→∞limi→∞hi,N=limN→∞∫∫∑j=1Nλj​ϕj​\(p\)​ϕj​\(q\)​𝑑μ​\(p\)​𝑑μ​\(q\)=∫∫⁡K⁡\(p,q\)​𝑑μ​\(p\)​𝑑μ​\(q\)=‖μ^‖H2,\\lim\_\{N\\to\\infty\}\\lim\_\{i\\to\\infty\}h\_\{i,N\}=\\lim\_\{N\\to\\infty\}\\int\\int\\sum\_\{j=1\}^\{N\}\\lambda\_\{j\}\\phi\_\{j\}\(p\)\\phi\_\{j\}\(q\)d\\mu\(p\)d\\mu\(q\)=\\int\\int K\(p,q\)d\\mu\(p\)d\\mu\(q\)=\\\|\\hat\{\\mu\}\\\|^\{2\}\_\{H\},\(28\)where in the second equality we used the Dominated Convergence Theorem with the bound

\|∑j=1Nλj​ϕj​\(p\)​ϕj​\(q\)\|≤K⁡\(p,p\)​K​\(q,q\)\.\\left\\lvert\\sum\_\{j=1\}^\{N\}\\lambda\_\{j\}\\phi\_\{j\}\(p\)\\phi\_\{j\}\(q\)\\right\\rvert\\leq\\sqrt\{K\(p,p\)K\(q,q\)\}\.\(29\)We also have

limi→∞limN→∞hi,N=limi→∞1Ti2​∫0Ti∫0TiK⁡\(y⁡\(s\),y⁡\(t\)\)​𝑑s​𝑑t=limi→∞‖Li∗‖H2,\\lim\_\{i\\to\\infty\}\\lim\_\{N\\to\\infty\}h\_\{i,N\}=\\lim\_\{i\\to\\infty\}\\frac\{1\}\{T\_\{i\}^\{2\}\}\\int\_\{0\}^\{T\_\{i\}\}\\int\_\{0\}^\{T\_\{i\}\}K\(y\(s\),y\(t\)\)dsdt=\\lim\_\{i\\to\\infty\}\\\|L\_\{i\}^\{\*\}\\\|^\{2\}\_\{H\},\(30\)using Dominated Convergence in the first equality\. Thus, we need to show that the order of the limits is exchangeable\. Lethih\_\{i\}denote the limitlimN→∞hi,N\\lim\_\{N\\to\\infty\}h\_\{i,N\}\. We have

\|hi,N−hi\|\\displaystyle\\lvert h\_\{i,N\}\-h\_\{i\}\\rvert=\|1Ti2​∫0Ti∫0Ti∑j=N\+1∞λj​ϕj​\(y⁡\(s\)\)​ϕj​\(y⁡\(t\)\)​𝑑s​𝑑t\|\\displaystyle=\\left\\lvert\\frac\{1\}\{T\_\{i\}^\{2\}\}\\int\_\{0\}^\{T\_\{i\}\}\\int\_\{0\}^\{T\_\{i\}\}\\sum\_\{j=N\+1\}^\{\\infty\}\\lambda\_\{j\}\\phi\_\{j\}\(y\(s\)\)\\phi\_\{j\}\(y\(t\)\)dsdt\\right\\rvert\(31\)≤1Ti2​∫0Ti∫0Timax\(y1,y2\)∈M×M⁡\|∑j=N\+1∞λj​ϕj​\(y1\)​ϕj​\(y2\)\|​𝑑s​𝑑t\\displaystyle\\leq\\frac\{1\}\{T\_\{i\}^\{2\}\}\\int\_\{0\}^\{T\_\{i\}\}\\int\_\{0\}^\{T\_\{i\}\}\\max\_\{\(y\_\{1\},y\_\{2\}\)\\in M\\times M\}\\left\\lvert\\sum\_\{j=N\+1\}^\{\\infty\}\\lambda\_\{j\}\\phi\_\{j\}\(y\_\{1\}\)\\phi\_\{j\}\(y\_\{2\}\)\\right\\rvert dsdt=max\(y1,y2\)∈M×M⁡\|∑j=N\+1∞λj​ϕj​\(y1\)​ϕj​\(y2\)\|→0,\\displaystyle=\\max\_\{\(y\_\{1\},y\_\{2\}\)\\in M\\times M\}\\left\\lvert\\sum\_\{j=N\+1\}^\{\\infty\}\\lambda\_\{j\}\\phi\_\{j\}\(y\_\{1\}\)\\phi\_\{j\}\(y\_\{2\}\)\\right\\rvert\\to 0,where on the last line we used the uniform convergence of the Mercer kernel onM×MM\\times M\. The bound on\|hi,N−hi\|\\lvert h\_\{i,N\}\-h\_\{i\}\\rvertdoes not depend oniiand converges to zero\. Thus, convergence ofhi,Nh\_\{i,N\}inNNis uniform, and we may apply the Moore\-Osgood Theorem to show that the order of the limits is exchangable\. The result follows\. ∎

###### Theorem A\.1\.

Let the setting of Proposition[A\.1](https://arxiv.org/html/2610.33354#A1.Thmproposition1)hold\. Then we have the following:

1Ti​∫0TiK⁡\(⋅,y⁡\(t\)\)​𝑑t→∫MK⁡\(⋅,p\)​𝑑μ​\(p\)\\frac\{1\}\{T\_\{i\}\}\\int\_\{0\}^\{T\_\{i\}\}K\(\\cdot,y\(t\)\)dt\\to\\int\_\{M\}K\(\\cdot,p\)d\\mu\(p\)\(32\)asi→∞i\\to\\infty\.

###### Proof\.

By Proposition[A\.1](https://arxiv.org/html/2610.33354#A1.Thmproposition1)and Lemma[A\.1](https://arxiv.org/html/2610.33354#A1.Thmlemma1), we have

‖Li∗−μ^‖H2=‖Li∗‖H2−2​⟨Li∗,μ^⟩H\+‖μ^‖H2→‖μ^‖H2−2​‖μ^‖H2\+‖μ^‖H2=0\.\\\|L\_\{i\}^\{\*\}\-\\hat\{\\mu\}\\\|^\{2\}\_\{H\}=\\\|L\_\{i\}^\{\*\}\\\|^\{2\}\_\{H\}\-2\\langle L\_\{i\}^\{\*\},\\hat\{\\mu\}\\rangle\_\{H\}\+\\\|\\hat\{\\mu\}\\\|^\{2\}\_\{H\}\\to\\\|\\hat\{\\mu\}\\\|\_\{H\}^\{2\}\-2\\\|\\hat\{\\mu\}\\\|^\{2\}\_\{H\}\+\\\|\\hat\{\\mu\}\\\|^\{2\}\_\{H\}=0\.\(33\)∎

For systems with sufficient regularity \(for example, hyperbolic systems; see\[[39](https://arxiv.org/html/2610.33354#bib.bib33)\]for details\), the time\-average of an observable over a trajectory converges to the spatial average with respect to the associated SRB measure at a rate of𝒪ℙ\(T−1/2\)\\mathcal\{O\}\_\{\\mathbb\{P\}\}\(T^\{\-1/2\}\)\. We use this fact to show that the time\-averaged occupation kernels with random initial conditions drawn from Lebesgue\-absolutely\-continuous distributions converge \(in probability\) to the KMEs of SRB measures at the same rate\.

###### Definition A\.1\.

Letf:M→Mf\\colon M\\to Mbe aCrC^\{r\}\(r≥1r\\geq 1\) map onMM\. LetΦt:M→M\\Phi\_\{t\}\\colon M\\to Mbe the flow offfat timett\. Define the centered\-kernel functionK¯:M×M→ℝ\\overline\{K\}\\colon M\\times M\\to\\mathbb\{R\}according to:

K¯​\(x,y\):=K⁡\(x,y\)−∫MK⁡\(x,p\)​𝑑μ​\(p\)−∫MK⁡\(q,y\)​𝑑μ​\(q\)\+∫M∫MK⁡\(p,q\)​𝑑μ​\(p\)​𝑑μ​\(q\)\.\\overline\{K\}\(x,y\):=K\(x,y\)\-\\int\_\{M\}K\(x,p\)d\\mu\(p\)\-\\int\_\{M\}K\(q,y\)d\\mu\(q\)\+\\int\_\{M\}\\int\_\{M\}K\(p,q\)d\\mu\(p\)d\\mu\(q\)\.\(34\)Define the auto\-correlation function forK¯\\overline\{K\}as:

CK¯​\(τ\):=∫MK¯​\(p,Φτ​\(p\)\)​𝑑μ​\(p\)C\_\{\\overline\{K\}\}\(\\tau\):=\\int\_\{M\}\\overline\{K\}\(p,\\Phi\_\{\\tau\}\(p\)\)d\\mu\(p\)\(35\)

First, we prove a convergence result for the case that the initial condition is distributed according toμ\\mu, an SRB measure forff\. We then extend the result to more general measures in the next result\.

###### Lemma A\.2\.

Assume we are in the setting of Definition[A\.1](https://arxiv.org/html/2610.33354#A1.Thmdefinition1)\. LetX∼μX\\sim\\mu, whereμ\\muis an SRB measure forff\. Suppose that the following regularity condition holds:

σK¯2:=2​∫0∞\|CK¯​\(τ\)\|​𝑑τ<∞\\sigma\_\{\\overline\{K\}\}^\{2\}:=2\\int\_\{0\}^\{\\infty\}\\lvert C\_\{\\overline\{K\}\}\(\\tau\)\\rvert d\\tau<\\infty\(36\)Then we have the following convergence in probability:

‖1T∫0TK\(⋅,Φt\(X\)\)dt−μ^‖H→ℙ0\\left\\\|\\frac\{1\}\{T\}\\int\_\{0\}^\{T\}K\(\\cdot,\\Phi\_\{t\}\(X\)\)dt\-\\hat\{\\mu\}\\right\\\|\_\{H\}\\to\_\{\\mathbb\{P\}\}0\(37\)at the rate𝒪ℙ\(T−1/2\)\\mathcal\{O\}\_\{\\mathbb\{P\}\}\(T^\{\-1/2\}\)\.

###### Proof\.

Denote byL∗​\(μ,T\)L^\{\*\}\(\\mu,T\)the following:

L∗​\(μ,T\):=1T​∫0TK⁡\(⋅,Φt​\(X\)\)​𝑑tL^\{\*\}\(\\mu,T\):=\\frac\{1\}\{T\}\\int\_\{0\}^\{T\}K\(\\cdot,\\Phi\_\{t\}\(X\)\)dt\(38\)whereX∼μX\\sim\\mu\. We have

‖L∗​\(μ,T\)−μ^‖H2=1T2​∫0T∫0TK¯​\(Φt​\(X\),Φs​\(X\)\)​𝑑s​𝑑t\.\\\|L^\{\*\}\(\\mu,T\)\-\\hat\{\\mu\}\\\|^\{2\}\_\{H\}=\\frac\{1\}\{T^\{2\}\}\\int\_\{0\}^\{T\}\\int\_\{0\}^\{T\}\\overline\{K\}\(\\Phi\_\{t\}\(X\),\\Phi\_\{s\}\(X\)\)dsdt\.\(39\)whereK¯\\overline\{K\}is the centered kernel in Definition[A\.1](https://arxiv.org/html/2610.33354#A1.Thmdefinition1)\. We have

𝔼μ​\[‖L∗​\(μ,T\)−μ^‖H2\]=1T2​∫0T∫0TCK¯​\(\|t−s\|\)​𝑑s​𝑑t\.\\mathbb\{E\}\_\{\\mu\}\\left\[\\left\\\|L^\{\*\}\(\\mu,T\)\-\\hat\{\\mu\}\\right\\\|^\{2\}\_\{H\}\\right\]=\\frac\{1\}\{T^\{2\}\}\\int\_\{0\}^\{T\}\\int\_\{0\}^\{T\}C\_\{\\overline\{K\}\}\(\\lvert t\-s\\rvert\)dsdt\.\(40\)by integrating the right side of \([39](https://arxiv.org/html/2610.33354#A1.E39)\) with respect toμ\\muoverMMand using theff\-invariance ofμ\\mu\. Taking the change of variablesτ=t−s\\tau=t\-s, we get

𝔼μ​\[‖L∗​\(μ,T\)−μ^‖H2\]=2T​∫0T\(1−τT\)​CK¯​\(τ\)​𝑑τ≤2T​∫0T\|CK¯​\(τ\)\|​𝑑τ\.\\mathbb\{E\}\_\{\\mu\}\\left\[\\left\\\|L^\{\*\}\(\\mu,T\)\-\\hat\{\\mu\}\\right\\\|^\{2\}\_\{H\}\\right\]=\\frac\{2\}\{T\}\\int\_\{0\}^\{T\}\\left\(1\-\\frac\{\\tau\}\{T\}\\right\)C\_\{\\overline\{K\}\}\(\\tau\)d\\tau\\leq\\frac\{2\}\{T\}\\int\_\{0\}^\{T\}\\lvert C\_\{\\overline\{K\}\}\(\\tau\)\\rvert d\\tau\.\(41\)Using Markov’s inequality in conjunction with the regularity condition, we get that

ℙ⁡\(‖L∗​\(μ,T\)−μ^‖H≥δ\)≤σK¯2δ2​T\.\\mathbb\{P\}\(\\\|L^\{\*\}\(\\mu,T\)\-\\hat\{\\mu\}\\\|\_\{H\}\\geq\\delta\)\\leq\\frac\{\\sigma\_\{\\overline\{K\}\}^\{2\}\}\{\\delta^\{2\}T\}\.\(42\)Setη:=σK¯2δ2​T\\eta:=\\frac\{\\sigma\_\{\\overline\{K\}\}^\{2\}\}\{\\delta^\{2\}T\}\. Then, with probability at least1−η1\-\\eta, we have

‖L∗​\(μ,T\)−μ^‖H≤σK¯η​T,\\\|L^\{\*\}\(\\mu,T\)\-\\hat\{\\mu\}\\\|\_\{H\}\\leq\\frac\{\\sigma\_\{\\overline\{K\}\}\}\{\\sqrt\{\\eta T\}\},\(43\)giving us a convergence rate of𝒪ℙ\(T−1/2\)\\mathcal\{O\}\_\{\\mathbb\{P\}\}\(T^\{\-1/2\}\)\. ∎

Suppose that the systemffis exponentially mixing, meaning that correlations decay exponentially\[[38](https://arxiv.org/html/2610.33354#bib.bib32)\]\. Then the auto\-correlation function of the centered kernel will be integrable for kernels which areα\\alpha\-Hölder continuous, forα∈\(0,1\]\\alpha\\in\(0,1\]\. This holds since, in this case,\|CK¯​\(t\)\|≤C​e−γ​t​‖K‖Cα\\lvert C\_\{\\overline\{K\}\}\(t\)\\rvert\\leq Ce^\{\-\\gamma t\}\\\|K\\\|\_\{C^\{\\alpha\}\}\[[7](https://arxiv.org/html/2610.33354#bib.bib34)\], and so2​∫0∞\|CK¯​\(τ\)\|​𝑑τ<∞2\\int\_\{0\}^\{\\infty\}\\lvert C\_\{\\overline\{K\}\}\(\\tau\)\\rvert d\\tau<\\infty\. Standard kernels such as Gaussian, Laplace, Matérn, and polynomial satisfy this condition\.

#### Setting\.

Letffbe aC1C^\{1\}vector field onℝd\\mathbb\{R\}^\{d\}with flowΦt\\Phi\_\{t\}, and letM⊂ℝdM\\subset\\mathbb\{R\}^\{d\}be compact and forward invariant, i\.e\.Φt​\(M\)⊂M\\Phi\_\{t\}\(M\)\\subset Mfor allt≥0t\\geq 0\(e\.g\. a trapping region containing the attractor\)\. Letμ\\mube aΦ\\Phi\-invariant Borel probability measure onMM, letKKbe a continuous positive\-definite kernel onMMwith RKHSHH, and writeKx:=K⁡\(⋅,x\)K\_\{x\}:=K\(\\cdot,x\)andμ^:=∫MKp​𝑑μ​\(p\)\\hat\{\\mu\}:=\\int\_\{M\}K\_\{p\}\\,d\\mu\(p\)\. Forx∈Mx\\in MandT\>0T\>0define

LT​\(x\):=1T​∫0TKΦt​\(x\)​𝑑t∈H\.L\_\{T\}\(x\):=\\frac\{1\}\{T\}\\int\_\{0\}^\{T\}K\_\{\\Phi\_\{t\}\(x\)\}\\,dt\\in H\.\(44\)SinceΦ\\Phiandx↦Kxx\\mapsto K\_\{x\}are continuous,x↦LT​\(x\)x\\mapsto L\_\{T\}\(x\)is continuous fromMMtoHH; in particular it is Borel measurable\.

#### Hypotheses\.

1. \(K\)*Hölder feature map\.*There existLK\>0L\_\{K\}\>0andκ∈\(0,1\]\\kappa\\in\(0,1\]such that‖Kx−Ky‖H≤LK​\|x−y\|κ\\\|K\_\{x\}\-K\_\{y\}\\\|\_\{H\}\\leq L\_\{K\}\|x\-y\|^\{\\kappa\}for allx,y∈Mx,y\\in M\.
2. \(L\)*Rate under the invariant measure\.*There existc≥0c\\geq 0andT0\>0T\_\{0\}\>0such that, forX~∼μ\\tilde\{X\}\\sim\\mu, 𝔼⁡\[‖LT​\(X~\)−μ^‖H2\]≤cTfor all​T≥T0\.\\mathbb\{E\}\\left\[\\\|L\_\{T\}\(\\tilde\{X\}\)\-\\hat\{\\mu\}\\\|\_\{H\}^\{2\}\\right\]\\leq\\frac\{c\}\{T\}\\qquad\\text\{for all \}T\\geq T\_\{0\}\.\(45\)

###### Assumption A\.1\(Asymptotic phase with absolutely continuous phase distribution\)\.

LetX∼νX\\sim\\nu, whereν\\nuis a Borel probability measure onMM\. There exist a random variableYYwith values inMM, defined on a common probability space withXX, a constantλ\>0\\lambda\>0, and a random variableC⁡\(X,Y\)C\(X,Y\)that is finite almost surely, such that

1. \(A1\)\|Φt​\(X\)−Φt​\(Y\)\|≤C⁡\(X,Y\)​e−λ​t\|\\Phi\_\{t\}\(X\)\-\\Phi\_\{t\}\(Y\)\|\\leq C\(X,Y\)\\,e^\{\-\\lambda t\}for allt≥0t\\geq 0, almost surely;
2. \(A2\)the law ofYYis absolutely continuous with respect toμ\\mu\.

Assumption[A\.1](https://arxiv.org/html/2610.33354#A1.Thmassumption1)says that the orbit of a typical initial condition is exponentially asymptotic to the orbit of a pointYYon the attractor \(an*asymptotic phase*ofXX\), and that the law ofYYis not concentrated on aμ\\mu\-null set\. For hyperbolic attractors this is the mechanism by which SRB measures are physical\.

###### Theorem A\.2\.

Assume \(K\), \(L\), and Assumption[A\.1](https://arxiv.org/html/2610.33354#A1.Thmassumption1)\. Then

‖1T∫0TK\(⋅,Φt\(X\)\)dt−μ^‖H=𝒪ℙ\(T−1/2\)\.\\left\\\|\\frac\{1\}\{T\}\\int\_\{0\}^\{T\}K\(\\cdot,\\Phi\_\{t\}\(X\)\)\\,dt\-\\hat\{\\mu\}\\right\\\|\_\{H\}=\\mathcal\{O\}\_\{\\mathbb\{P\}\}\\\!\\left\(T^\{\-1/2\}\\right\)\.\(46\)More precisely, for everyε\>0\\varepsilon\>0there existsDε<∞D\_\{\\varepsilon\}<\\inftysuch that

ℙ\(∥LT\(X\)−μ^∥H\>DεT−1/2\)≤εfor allT≥max\(1,T0\)\.\\mathbb\{P\}\\left\(\\\|L\_\{T\}\(X\)\-\\hat\{\\mu\}\\\|\_\{H\}\>D\_\{\\varepsilon\}T^\{\-1/2\}\\right\)\\leq\\varepsilon\\qquad\\text\{for all \}T\\geq\\max\(1,T\_\{0\}\)\.\(47\)

###### Proof\.

By the triangle inequality,

‖LT​\(X\)−μ^‖H≤‖LT​\(X\)−LT​\(Y\)‖H⏟\(I\)\+‖LT​\(Y\)−μ^‖H⏟\(II\)\.\\\|L\_\{T\}\(X\)\-\\hat\{\\mu\}\\\|\_\{H\}\\leq\\underbrace\{\\\|L\_\{T\}\(X\)\-L\_\{T\}\(Y\)\\\|\_\{H\}\}\_\{\\text\{\(I\)\}\}\+\\underbrace\{\\\|L\_\{T\}\(Y\)\-\\hat\{\\mu\}\\\|\_\{H\}\}\_\{\\text\{\(II\)\}\}\.\(48\)
*Term \(I\)\.*By \(K\) and \(A1\), almost surely,

\(I\)≤1T​∫0T‖KΦt​\(X\)−KΦt​\(Y\)‖H​𝑑t≤LK​C​\(X,Y\)κT​∫0Te−κ​λ​t​𝑑t\\displaystyle\\leq\\frac\{1\}\{T\}\\int\_\{0\}^\{T\}\\\|K\_\{\\Phi\_\{t\}\(X\)\}\-K\_\{\\Phi\_\{t\}\(Y\)\}\\\|\_\{H\}\\,dt\\leq\\frac\{L\_\{K\}\\,C\(X,Y\)^\{\\kappa\}\}\{T\}\\int\_\{0\}^\{T\}e^\{\-\\kappa\\lambda t\}\\,dt≤LK​C​\(X,Y\)κκ​λ​T\.\\displaystyle\\leq\\frac\{L\_\{K\}\\,C\(X,Y\)^\{\\kappa\}\}\{\\kappa\\lambda T\}\.\(49\)
*Term \(II\)\.*Letρ:=d​Law​\(Y\)/d​μ\\rho:=d\\,\\mathrm\{Law\}\(Y\)/d\\mu, which exists by \(A2\) and is finiteμ\\mu\-almost everywhere\. Forδ\>0\\delta\>0letAδ:=\{y∈M:‖LT​\(y\)−μ^‖H\>δ\}A\_\{\\delta\}:=\\\{y\\in M:\\\|L\_\{T\}\(y\)\-\\hat\{\\mu\}\\\|\_\{H\}\>\\delta\\\}, a Borel set\. By Chebyshev’s inequality and \(L\),μ⁡\(Aδ\)≤c/\(δ2​T\)\\mu\(A\_\{\\delta\}\)\\leq c/\(\\delta^\{2\}T\)forT≥T0T\\geq T\_\{0\}\. Hence, for anyB\>0B\>0,

ℙ\(\(II\)\>δ\)=∫Aδρdμ≤∫\{ρ\>B\}ρdμ\+Bμ\(Aδ\)≤ℙ\(ρ\(Y\)\>B\)\+B​cδ2​T\.\\mathbb\{P\}\(\\text\{\(II\)\}\>\\delta\)=\\int\_\{A\_\{\\delta\}\}\\rho\\,d\\mu\\leq\\int\_\{\\\{\\rho\>B\\\}\}\\rho\\,d\\mu\+B\\,\\mu\(A\_\{\\delta\}\)\\leq\\mathbb\{P\}\\big\(\\rho\(Y\)\>B\\big\)\+\\frac\{Bc\}\{\\delta^\{2\}T\}\.\(50\)
*Conclusion\.*Fixε\>0\\varepsilon\>0\. Sinceρ⁡\(Y\)<∞\\rho\(Y\)<\\inftyandC⁡\(X,Y\)<∞C\(X,Y\)<\\inftyalmost surely, chooseBBwithℙ⁡\(ρ⁡\(Y\)\>B\)≤ε/3\\mathbb\{P\}\(\\rho\(Y\)\>B\)\\leq\\varepsilon/3andbbwithℙ⁡\(C⁡\(X,Y\)\>b\)≤ε/3\\mathbb\{P\}\(C\(X,Y\)\>b\)\\leq\\varepsilon/3\. Seta:=3​B​c/εa:=\\sqrt\{3Bc/\\varepsilon\}\. Takingδ=aT−1/2\\delta=aT^\{\-1/2\}above givesℙ\(\(II\)\>aT−1/2\)≤2ε/3\\mathbb\{P\}\(\\text\{\(II\)\}\>aT^\{\-1/2\}\)\\leq 2\\varepsilon/3\. On the event\{C\(X,Y\)≤b\}\\\{C\(X,Y\)\\leq b\\\}and forT≥1T\\geq 1, term \(I\) is at mostLK​bκ/\(κ​λ​T\)≤LK​bκ/\(κ​λ​T\)L\_\{K\}b^\{\\kappa\}/\(\\kappa\\lambda T\)\\leq L\_\{K\}b^\{\\kappa\}/\(\\kappa\\lambda\\sqrt\{T\}\)\. WithDε:=a\+LK​bκ/\(κ​λ\)D\_\{\\varepsilon\}:=a\+L\_\{K\}b^\{\\kappa\}/\(\\kappa\\lambda\), a union bound givesℙ\(∥LT\(X\)−μ^∥H\>DεT−1/2\)≤ε\\mathbb\{P\}\(\\\|L\_\{T\}\(X\)\-\\hat\{\\mu\}\\\|\_\{H\}\>D\_\{\\varepsilon\}T^\{\-1/2\}\)\\leq\\varepsilonfor allT≥max⁡\(1,T0\)T\\geq\\max\(1,T\_\{0\}\)\. ∎

### A\.1Verification for the kernels and systems of Section[3\.1](https://arxiv.org/html/2610.33354#S3.SS1)

#### Kernels\.

All kernels of Appendix[C](https://arxiv.org/html/2610.33354#A3)satisfy \(K\)\. This follows from‖Kx−Ky‖H2=K⁡\(x,x\)\+K⁡\(y,y\)−2​K​\(x,y\)\\\|K\_\{x\}\-K\_\{y\}\\\|\_\{H\}^\{2\}=K\(x,x\)\+K\(y,y\)\-2K\(x,y\)together with1−e−u≤u1\-e^\{\-u\}\\leq uforu≥0u\\geq 0, and, for the quadratic kernel, from its explicit feature mapx↦\(1,2​x,x​x⊤\)x\\mapsto\(1,\\sqrt\{2\}\\,x,xx^\{\\top\}\)\. WithR:=maxx∈M⁡\|x\|R:=\\max\_\{x\\in M\}\|x\|:

- •GU with bandwidthη\\eta:LK=1/ηL\_\{K\}=1/\\eta,κ=1\\kappa=1;
- •GM with bandwidthsσi\\sigma\_\{i\}:LK=1/mini⁡σiL\_\{K\}=1/\\min\_\{i\}\\sigma\_\{i\},κ=1\\kappa=1;
- •Laplace with bandwidthsσi\\sigma\_\{i\}:LK=\(2​d/mini⁡σi\)1/2L\_\{K\}=\(2\\sqrt\{d\}/\\min\_\{i\}\\sigma\_\{i\}\)^\{1/2\},κ=1/2\\kappa=1/2;
- •quadratic:LK=2\+2​RL\_\{K\}=\\sqrt\{2\}\+2R,κ=1\\kappa=1\.

The bandwidths are computed once from the data and then held fixed; the results apply conditionally on them\.

#### FitzHugh–Nagumo\.

For both parameter sets the unique equilibriumppis a hyperbolic sink \(Jacobian eigenvalues−0\.104±0\.28​i\-0\.104\\pm 0\.28iand\{−3\.46,−0\.089\}\\\{\-3\.46,\-0\.089\\\}\), and numerically all initial conditions in\[−2,2\]2\[\-2,2\]^\{2\}converge to it\. Henceμ=δp\\mu=\\delta\_\{p\}andLT​\(p\)=μ^L\_\{T\}\(p\)=\\hat\{\\mu\}, so \(L\) holds withc=0c=0\. Convergence toppis exponential and uniform on compact subsets of the basin, so Assumption[A\.1](https://arxiv.org/html/2610.33354#A1.Thmassumption1)holds withY:=pY:=p\.

#### Van der Pol\.

Forμ\>0\\mu\>0there is a unique hyperbolic attracting limit cycleΓ\\Gamma, of periodPP, attracting every point except the origin\. Fory∈Γy\\in\\Gamma, each full period contributes exactlyP​μ^P\\hat\{\\mu\}, so‖LT​\(y\)−μ^‖H≤2​P​κ¯/T\\\|L\_\{T\}\(y\)\-\\hat\{\\mu\}\\\|\_\{H\}\\leq 2P\\bar\{\\kappa\}/Twithκ¯:=supx‖Kx‖H\\bar\{\\kappa\}:=\\sup\_\{x\}\\\|K\_\{x\}\\\|\_\{H\}, and \(L\) holds withc=4​P2​κ¯2c=4P^\{2\}\\bar\{\\kappa\}^\{2\}\. \(HereCK¯C\_\{\\overline\{K\}\}is periodic, soσK¯2=∞\\sigma\_\{\\overline\{K\}\}^\{2\}=\\inftyand Lemma[A\.2](https://arxiv.org/html/2610.33354#A1.Thmlemma2)does not apply as stated; cf\. Remark[A\.1](https://arxiv.org/html/2610.33354#A1.Thmremark1)\.\) TakeYYto be the asymptotic phase ofXXonΓ\\Gamma: \(A1\) holds by hyperbolicity ofΓ\\Gamma, and since the asymptotic\-phase map is aC1C^\{1\}submersion on the basin\[[11](https://arxiv.org/html/2610.33354#bib.bib40)\], \(A2\) holds whenν\\nuhas a density\.

## Appendix BApproximating the dissimilarity measure

We use integral quadratures to approximate the integrals in the dissimilarity measure in[9](https://arxiv.org/html/2610.33354#S2.E9)\. Assumingmi=mj=mm\_\{i\}=m\_\{j\}=m,

d^2​\(yi,yj\)=1Ti2​∑k=0m∑l=0mwi​k​wi​l​K​\(yi​\(ti​k\),yi​\(ti​l\)\)\+1Tj2​∑k=0m∑l=0mwj​k​wj​l​K​\(yj​\(tj​k\),yj​\(tj​l\)\)−2Ti​Tj​∑k=0m∑l=0mwi​k​wj​l​K​\(yi​\(ti​k\),yj​\(tj​l\)\),\\hat\{d\}^\{2\}\(y\_\{i\},y\_\{j\}\)=\\frac\{1\}\{T\_\{i\}^\{2\}\}\\sum\_\{k=0\}^\{m\}\\sum\_\{l=0\}^\{m\}w\_\{ik\}w\_\{il\}K\(y\_\{i\}\(t\_\{ik\}\),y\_\{i\}\(t\_\{il\}\)\)\+\\frac\{1\}\{T\_\{j\}^\{2\}\}\\sum\_\{k=0\}^\{m\}\\sum\_\{l=0\}^\{m\}w\_\{jk\}w\_\{jl\}K\(y\_\{j\}\(t\_\{jk\}\),y\_\{j\}\(t\_\{jl\}\)\)\-\\\\ \\frac\{2\}\{T\_\{i\}T\_\{j\}\}\\sum\_\{k=0\}^\{m\}\\sum\_\{l=0\}^\{m\}w\_\{ik\}w\_\{jl\}K\(y\_\{i\}\(t\_\{ik\}\),y\_\{j\}\(t\_\{jl\}\)\),\(52\)The trapezoidal quadrature weightswi​kw\_\{ik\}for theii\-th time series \(and similarly for thejj\-th one\) are defined as:

wi​k=\{ti,1−ti,02if​k=0,ti,k\+1−ti,k−12if​0<k<m,ti,m−ti,m−12if​k=m\.w\_\{ik\}=\\begin\{cases\}\\frac\{t\_\{i,1\}\-t\_\{i,0\}\}\{2\}&\\text\{if \}k=0,\\\\\[6\.0pt\] \\frac\{t\_\{i,k\+1\}\-t\_\{i,k\-1\}\}\{2\}&\\text\{if \}0<k<m,\\\\\[6\.0pt\] \\frac\{t\_\{i,m\}\-t\_\{i,m\-1\}\}\{2\}&\\text\{if \}k=m\.\\end\{cases\}\(53\)If the time series is sampled uniformly such thatΔ​ti=Tim\\Delta t\_\{i\}=\\frac\{T\_\{i\}\}\{m\}, the weights simplify to:

wi​k=\{Δ​ti2if​k=0​or​k=m,Δ​tiif​0<k<m\.w\_\{ik\}=\\begin\{cases\}\\frac\{\\Delta t\_\{i\}\}\{2\}&\\text\{if \}k=0\\text\{ or \}k=m,\\\\\[6\.0pt\] \\Delta t\_\{i\}&\\text\{if \}0<k<m\.\\end\{cases\}\(54\)

## Appendix CKernels used

We used several different kernels in the numerical experiments in Section[3](https://arxiv.org/html/2610.33354#S3)\. We list their definitions here\.

1. 1\.Gaussian uniscale \(GU\): K⁡\(x,y\)=exp⁡\(−12​η2​‖x−y‖2\),K\(x,y\)=\\exp\\left\(\-\\frac\{1\}\{2\\eta^\{2\}\}\\\|x\-y\\\|^\{2\}\\right\),\(55\)where we setη:=median​\{‖xk−xl‖\}k≠l\\eta:=\\mbox\{\\rm median\}\\\{\\\|x\_\{k\}\-x\_\{l\}\\\|\\\}\_\{k\\neq l\}wherexkx\_\{k\}ranged over each individual observation in the dataset\. This kernel was used for high\-dimensional data\.
2. 2\.Gaussian multiscale \(GM\): K\(x,y\)=exp\(−∑i=1d\(xi−yi\)22​σi2\),K\(x,y\)=\\exp\\left\(\-\\sum\_\{i=1\}^\{d\}\\frac\{\(x\_\{i\}\-y\_\{i\}\)^\{2\}\}\{2\\sigma\_\{i\}^\{2\}\}\\right\),\(56\)whereσi\\sigma\_\{i\}was set to be the standard deviation of theit​hi^\{th\}coordinate of the dataset\. This is the default kernel for CLOCK\.
3. 3\.Laplace: K\(x,y\)=exp\(−∑i=1d\|xi−yi\|σi\),K\(x,y\)=\\exp\\left\(\-\\sum\_\{i=1\}^\{d\}\\frac\{\\lvert x\_\{i\}\-y\_\{i\}\\rvert\}\{\\sigma\_\{i\}\}\\right\),\(57\)whereσi\\sigma\_\{i\}was set to be the standard deviation of theit​hi^\{th\}coordinate of the dataset\.
4. 4\.Quadratic \(quad\): K⁡\(x,y\)=\(xT​y\+1\)2\.K\(x,y\)=\(x^\{T\}y\+1\)^\{2\}\.\(58\)

## Appendix DFormulas for evaluation metrics

We provide explicit formulas for the evaluation metrics in Section[3\.3](https://arxiv.org/html/2610.33354#S3.SS3)\.

1. 1\.Adjusted Rand Index \(ARI\):The formula is given by ARI=∑i,j\(ni​j2\)−∑i\(ai2\)​∑j\(bj2\)\(n2\)12​\(∑i\(ai2\)\+∑j\(bj2\)\)−∑i\(ai2\)​∑j\(bj2\)\(n2\),\\text\{ARI\}=\\frac\{\\sum\_\{i,j\}\\binom\{n\_\{ij\}\}\{2\}\-\\frac\{\\sum\_\{i\}\\binom\{a\_\{i\}\}\{2\}\\sum\_\{j\}\\binom\{b\_\{j\}\}\{2\}\}\{\\binom\{n\}\{2\}\}\}\{\\frac\{1\}\{2\}\\left\(\\sum\_\{i\}\\binom\{a\_\{i\}\}\{2\}\+\\sum\_\{j\}\\binom\{b\_\{j\}\}\{2\}\\right\)\-\\frac\{\\sum\_\{i\}\\binom\{a\_\{i\}\}\{2\}\\sum\_\{j\}\\binom\{b\_\{j\}\}\{2\}\}\{\\binom\{n\}\{2\}\}\},\(59\)wherennis the total number of time series in the dataset,ni​jn\_\{ij\}is the number of trajectories that are in both the true clusteriiand the predicted clusterjj,aia\_\{i\}is the total number of samples in true clusterii, andbjb\_\{j\}is the total number of samples in predicted clusterjj\.
2. 2\.Normalized Mutual Information \(NMI\):The formula is given by NMI=2​∑i,jni​j​log⁡\(n⋅ni​jai​bj\)−∑iailog\(ain\)−∑jbjlog\(bjn\),\\text\{NMI\}=\\frac\{2\\sum\_\{i,j\}n\_\{ij\}\\log\\left\(\\frac\{n\\cdot n\_\{ij\}\}\{a\_\{i\}b\_\{j\}\}\\right\)\}\{\-\\sum\_\{i\}a\_\{i\}\\log\\left\(\\frac\{a\_\{i\}\}\{n\}\\right\)\-\\sum\_\{j\}b\_\{j\}\\log\\left\(\\frac\{b\_\{j\}\}\{n\}\\right\)\},\(60\)wheren,ni​j,ai,bjn,n\_\{ij\},a\_\{i\},b\_\{j\}are as above and by convention0​log⁡00\\log 0is set to00\.
3. 3\.Purity:IfΩ\\Omegais a predicted clustering andCCa true clustering, then the purity is defined as purity​\(Ω,C\)=∑ω∈Ωmaxc∈C⁡\|ω∩c\|\.\\mbox\{\\rm purity\}\(\\Omega,C\)=\\sum\_\{\\omega\\in\\Omega\}\\max\_\{c\\in C\}\\lvert\\omega\\cap c\\rvert\.\(61\)In the setting of two predicted clusters and two true cluseters, the purity is just the accuracy of the clustering\.

## Appendix EHyperparameter tuning

We validated the hyperparameters of the deep learning methods in Section[3\.4](https://arxiv.org/html/2610.33354#S3.SS4)using an internal clustering metric, the DBCV\.

1. 1\.Density\-based clustering validation \(DBCV\):The DBCV\[[20](https://arxiv.org/html/2610.33354#bib.bib24)\]is an internal clustering metric, meaning it can be used in the absence of true cluster labels\. This metric assigns a high value to partitions with highly\-dense clusters with sparse between\-cluster regions\. It does so without assuming the clusters are globular, and so can measure the quality of clustering for highly\-irregular collections of trajectories, which occurs frequently with chaotic systems, for example\. The DBCV score is normalized to be in the range\[−1,1\]\[\-1,1\]with 1 being a perfect clustering and \-1 being exceptionally poor\. We found it to be superior to the silhouette score\[[31](https://arxiv.org/html/2610.33354#bib.bib13)\]for the task of validating hyperparameters\.

## Appendix FNoise robustness experiment

We speculate that the CLOCK method is highly resistant to zero\-mean Gaussian noise since the feature vector is obtained by performing an integration along a trajectory\. The integral effectively cancels out the perturbations caused by the noise\. In order to gauge the impact of noise on the performance of CLOCK, we performed an experiment in which we gradually increased the noise level in the Lorenz96 dataset while leaving the true trajectories fixed\. The CLOCK method with Gaussian uniscale and quadratic kernels \(see Appendix[C](https://arxiv.org/html/2610.33354#A3)\) was fit on each dataset, and we evaluated the quality of the clusterings\. We ran this experiment with the baseline method TS2Vec to use as a comparison\. The results are graphically represented in Figure[3](https://arxiv.org/html/2610.33354#A6.F3)\.

The results of the noise robustness experiment show that CLOCK is resistant to noise\. It should be noted that the amount of noise added to the trajectories was extreme, even on the lower settings\. In Figure[3](https://arxiv.org/html/2610.33354#A6.F3), we see that CLOCK with Gaussian kernel fails as the noise level approaches 3\.0, while CLOCK with quadratic kernel is affected at level 2\.0\. Thus, the choice of kernel is an important aspect of this clustering method\. In practice, the kernel may be treated as a hyperparameter and validated using an internal clustering metric\. TS2Vec generally performed worse than CLOCK with either kernel\.

\(a\)ARI\(b\)NMI\(c\)Purity
Figure 3:Plots of evaluation metrics in the noise robustness test of Appendix[F](https://arxiv.org/html/2610.33354#A6)\. The horizontal axis represents the standard deviation of the Gaussian noise added to each observation\. The vertical axis represents one of the evaluation metrics: ARI, NMI or purity\. In Figure[3\(a\)](https://arxiv.org/html/2610.33354#A6.F3.sf1), we have ARI results, in Figure[3\(b\)](https://arxiv.org/html/2610.33354#A6.F3.sf2), we have results for NMI, and in Figure[3\(c\)](https://arxiv.org/html/2610.33354#A6.F3.sf3), we have purity results\.
## Appendix GAssets used

The assets we used can be divided into methods, back\-ends, and datasets\.

### G\.1Methods

1. 1\.Dynamic Time Warping \(DTW\):We used thedtwfunction from the Python packagetslearn\.metrics\. This has the BSD\-2\-Clause license\. We used version 0\.8\.0\.
2. 2\.Shape\-Based Distance \(SBD\):We used the packagescipy\.signaland the functioncorrelate\. This has the BSD 3\-Clause license\. We used version 1\.17\.1\.
3. 3\.Catch\-22:We used the packagepycatch22and the functioncatch22\_all\. This has the GNU General Public License v3 \(GPLv3\) license\. We used version 0\.4\.5\.
4. 4\.Latent ODE:We used thepytorchpackage and thetorchdiffeqpackage\. Thepytorchpackage has the BSD 3\-Clause license\. We used version 2\.10\.0\. Thetorchdiffeqpackage has the MIT license\. We used version 0\.2\.5\.
5. 5\.TS2Vec:Code was downloaded from the GitHub repository at[https://github\.com/zhihanyue/ts2vec](https://github.com/zhihanyue/ts2vec)\. This has the MIT license\. We used the AAAI 2022 Official Release\.
6. 6\.Spectral clustering:We usedsklearn\.clusterpackage and the functionSpectralClustering\. This has the BSD 3\-Clause license\. We used version 1\.8\.0\.

### G\.2Back\-ends

1. 1\.JAX:We used thejaxandjaxlibPython packages\. These have the Apache 2\.0 license\. We used version 4\.30\.0 ofjaxand version 4\.30\.0 ofjaxlib\.

### G\.3Datasets

1. 1\.