Time Series Classification through Diffeomorphic Time Warping (DiffTW) Vicky Geneva Haney∗1 , Kamel Lahouel†2 , Victor Rielly‡2 , and Bruno M. Jedynak§1
arXiv:2606.23472v1 [stat.ML] 22 Jun 2026
1 Department of Mathematics and Statistics, Portland State University, USA 2 The Translational Genomics Research Institute (TGen), USA
Abstract Time series classification involves learning a mapping from a continuous, temporally ordered sequence of real-valued observations to a discrete response variable, like class labels. This task is fundamental in domains, including health monitoring, where the temporal structure of data is critical for accurate prediction. Dynamic Time Warping (DTW) is a standard technique for measuring similarity between sequences varying in time or speed. However, DTW is restricted to discrete point matching. To move beyond pairwise alignment, we propose a theoretical framework that learns mappings between real-valued functions. These mappings approximate the flow associated with the characteristic curves of a linear transport equation with a space-dependent velocity field, providing a diffeomorphic transformation between two time series. Using the method of characteristics, we transform this partial differential equation into ordinary differential equations (ODEs) modeling system dynamics. The objective function used to learn these ODEs derives from the fundamental theorem of calculus. To enable flexible, expressive representations of the velocity field, we utilize reproducing kernel Hilbert spaces and optimal control methods. Our method, Diffeomorphic Time Warping (DiffTW), provides a theoretically grounded dissimilarity measure. Using a 1-nearest neighbor classifier, DiffTW outperforms DTW on 60 of 86 datasets. Keywords Backpropagation; DTW; Nearest neighbor classification; RKHS ∗
Corresponding author. Email: [email protected] Corresponding author. Email: [email protected] ‡ Corresponding author. Email: [email protected] § Corresponding author. Email: [email protected] †
1
2
Vicky Geneva Haney, Kamel Lahouel, Victor Rielly, and Bruno M. Jedynak
1
Introduction
Dynamic Time Warping (DTW) is a widely used method for comparing time series that may differ in sampling rate or temporal distortion. DTW computes an optimal alignment path based on a cumulative cost matrix, and the resulting alignment cost serves as a dissimilarity measure Giorgino (2009); Raghavendra et al. (2011). While DTW is effective, it remains a pairwise alignment method and does not model the underlying process generating the sequences. Many extensions, such as modified distances, regularization, and differentiable variants, Liu et al. (2024a,b); Cuturi and Blondel (2017); Buchin et al. (2023), still operate within the alignment framework. Instead of aligning sequences, we learn the dynamics transforming one signal into another. For processes with continuous temporal deformations, we represent these transformations using a first-order partial differential equations (PDE), enabling smooth maps and a dynamical view of time series. Using the method of characteristics, we reduce the PDE to a system of ordinary differential equations (ODEs) and learn the associated vector field in a reproducing kernel Hilbert space (RKHS). Our approach is conceptually related to large deformation diffeomorphic metric mapping (LDDMM) Beg et al. (2005). Our approach differs in two key aspects: we use a time-independent vector field and an explicit kernel composed of random Fourier features Rahimi and Recht (2007). We call the resulting algorithm Diffeomorphic Time Warping (DiffTW). We evaluated DiffTW on electrocardiogram (ECG) classification tasks where DTW has been extensively studied Raghavendra et al. (2011); Middlehurst et al. (2024); Dau et al. (2019). Using standard ECG datasets and subsets of the University of California, Riverside (UCR) archive Dau et al. (2019), we compared DiffTW against a standard Python implementation of DTW Giorgino (2009). While DTW has quadratic complexity in sequence length, DiffTW scales linearly due to its continuous formulation as shown and discussed in Section 3.5 and avoids computation of a full alignment matrix. This property suggests the potential for handling large and high-resolution time series. The rest of this paper is structured as follows: Section 2 describes the methodology we
Diffeomorphic Time Warping (DiffTW)
3
used to learn the vector fields. Section 3 details the experiments we conducted using this approach. Section 4 provides an overview of our results compared to DTW, and finally, Section 5 summarizes our work and discusses its impact on time series classification.
2
Methods
2.1
Defining the Partial Differential Equation
We aim to transition from measuring similarity between sequences as done in DTW to modeling the underlying dynamics and framing the alignment process as a continuous physical flow governed by a first-order advection partial differential equation (PDE). Let u : [0, 1] × [0, 1] 7→ R represent the signal evolving from a source sequence ϕ0 to a target sequence ϕ1 . The boundary states of this evolution are defined as:
u(x, 0) = ϕ0 (x),
u(x, 1) = ϕ1 (x)
This evolution is driven by the advection equation:
α(x)ux + ut = 0
where ut and ux denote the parital derivaties of u with respect to time t and space x, respectively. The velocity field α(x) dictates the rate at which the shape of u propagates. In our formulation, we require α to be a function belonging to a reproducing kernel Hilbert space (RKHS). Using the method of characteristics McOwen (2003), this PDE reduces to a system of ordinary differential equations (ODEs). Let x(t), denote the spatial coordinate along a characteristic curve; for notational simplicity, we will frequently write x in place of x(t). The reduced ODE system is given by:
∂x = ẋ = α(x), ∂t
∂u = u̇ = 0 ∂t
(1)
4
Vicky Geneva Haney, Kamel Lahouel, Victor Rielly, and Bruno M. Jedynak
Equation (1) demonstrates that u is constant along the characteristics. Refer to Figure 2.1 for a description of our characteristics and to explicitly show what we aim to learn. u
ϕ0 (xi (0))
ϕ0 (xi (2/3)) = ϕ0 (xi (0)) ϕ0 (xi (1/3)) = ϕ0 (xi (0))
ϕ1 (xi (1)) = ϕ0 (xi (0))
t=0
t=1
xi (0)
xi (2/3) xi (1/3)
t
xi (1)
ϕ0 x
ϕ1
Figure 1: Visualization of continuous time warping via the method of characteristics. Points xi (0) from the reference signal are transported along characteristic curves xi (t) over time t ∈ [0, 1]. By construction, the signal value u is conserved along these trajectories, meaning ϕ0 (xi (t)) = ϕ0 (xi (0)). The proposed model optimizes the underlying vector field α(x) to construct the specific paths xi (t) that continuously deform ϕ0 into optimal alignment with the target signal ϕ1 .
2.2
Deriving the Objective Function (i)
(i)
(i)
(i)
Given a set of n training pairs {(ϕ0 , ϕ1 )}ni=1 , where we define ϕ0 = ϕ0 (xi (0)) and ϕ1 = ϕ1 (xi (1)), we want to formulate an objective function that encourages each transported point to match its target location at time t = 1. Therefore, we define the following objective function to minimize: n
1X J(α) = [ϕ1 (xi (1)) − ϕ0 (xi (0))]2 + λ∥α∥2 , n i=1
(2)
subject to the constraint ẋi (t) = α(xi (t)), for t ∈ [0, 1], where λ > 0 is a regularization parameter. By the Fundamental Theorem of Calculus, the spatial position at t = 1 is given by: Z 1 xi (1) = xi (0) +
α(xi (t)) dt. 0
Diffeomorphic Time Warping (DiffTW)
5
Because xi (1) depends implicitly on the velocity field α, computing the gradient of J(α) to find the optimal α requires specific care. Sections 2.3 and 2.4 detail the explicit derivation of this gradient.
2.3
Reproducing Kernel Hilbert Spaces
We employ theory of reproducing kernel Hilbert spaces (RKHSs) in order to learn α. Let X be a nonempty set. A Hilbert space of real-valued functions on X that allows the pointwise evaluation at any point x ∈ X to be a continuous linear functional in its norm topology is called an RKHS Pereverzyev (2022). Definition 1 (Positive Definite Kernel). Let X be a non-empty set. A function k : X × X → R is called a positive definite kernel on X if and only if it is symmetric, that is k(x, y) = k(y, x) for any x, y ∈ X , and positive definite, that is n X n X
ci cj k(xi , xj ) ≥ 0
i=1 j=1
for any n > 0, any choice of n points x1 , x2 , . . . , xn ∈ X , and any choice of real numbers c1 , . . . , cn ∈ R. Definition 2 (RKHS). Let X be a set and H ⊂ RX (that is, let H be a subset of functions that take elements from X → R) where H is a class of functions that form a Hilbert space with inner product ⟨·, ·⟩H . The function k : X × X 7→ R is called a reproducing kernel of H if 1. H contains all functions of the form
∀x ∈ X , kx : y 7→ k(x, y). 2. For every x ∈ X and f ∈ H the reproducing property holds:
f (x) = ⟨f, kx ⟩H .
A fundamental property of an RKHS is that there exists a continuous feature map-
6
Vicky Geneva Haney, Kamel Lahouel, Victor Rielly, and Bruno M. Jedynak
ping Ψ : X → H such that the kernel evaluates the inner product in H:
k(x, x′ ) = ⟨Ψ(x), Ψ(x′ )⟩H .
While this feature space is often infinite-dimensional, Bochner’s Theorem allows us to approximate shift-invariant kernels using finite-dimensional mappings. In this work, we utilize the Gaussian kernel,
k(x, x′ ) = exp
−∥x − x′ ∥2 , 2σ 2
where σ is the bandwidth parameter. To approximate this kernel explicitly, we apply random Fourier features Rahimi and Recht (2007). We draw D frequency scalars (where D is the number of features), ω1 , . . . , ωD , independently from a standard normal distribution, ω ∼ N (0, 1). By evaluating the trigonometric transformations scaled by our bandwidth σ, we define our explicit feature map γ(x) ∈ R2D (where p = 2D) as the concatenation of cosine and sine functions:
cos(ω1 x/σ) .. . cos(ω x/σ) D 1 . γ(x) = √ D sin(ω1 x/σ) .. . sin(ωD x/σ)
(3)
By combining Bochner’s Theorem with the strong law of large numbers, the empirical inner product of these random features converges to the exact kernel function as the dimension D approaches infinity. Consequently, for any x, x′ ∈ X , we have:
limD→+∞ γ(x)T γ(x′ ) = k(x, x′ ), a.s. Hence, we can express our target function α ∈ H as a linear combination of these features
Diffeomorphic Time Warping (DiffTW)
7
using learnable weights β ∈ Rp : α(x) = β T γ(x). This reduces the infinite-dimensional representation of α to a low-dimensional, finite parameterization. Consequently, the objective function, Equation (2), becomes: n
1X J(β) = [ϕ1 (xi (1)) − ϕ0 (xi (0))]2 + λβ T β, n i=1
(4)
with the constraint ẋi = β T γ(xi ). The vector β is then learned directly.
2.4
Learning the Characteristics
Computing the gradient of the objective function, ∇β J using Equation (4), presents a computational challenge because the learnable parameters β are implicitly nested within the integration of the final state, xi (1). To overcome this, we frame the gradient computation as an optimal control problem by augmenting the objective function with continuous Lagrange multipliers, pi (t). By applying integration by parts to the system’s dynamical constraint (ẋi = α), we derive a secondary ODE governing the evolution of the multiplier, ṗi (t), subject to a defined terminal condition, pi (1). This formulation enables a continuous backpropagation procedure; by solving the resulting adjoint ODE for pi (t) backwards in time from t = 1 to t = 0, we can efficiently compute the exact gradients in a continuous-time framework that mathematically mirrors backpropagation through time in deep learning Chen et al. (2018). Consider the functional: n 1X T T ψ β γ(·), x0 = ψi β γ(·), x0 , n i=1 with
T
ψi β γ(·), x0
h i2 = ϕ1 (xi (1)) − ϕ0 (xi (0)) ,
8
Vicky Geneva Haney, Kamel Lahouel, Victor Rielly, and Bruno M. Jedynak
for i = {1, . . . , n} under the constraints defined by the ẋ(t) = α(x) = β T γ(x(t)) and x(0) = x0 . Then for each i: ẋi (t) = β T γ(xi (t)), xi (0) = xi (0). Next, we augment our objective function with Lagragian methods to replace this constrained optimization problem and define the Lagrange multipliers pi (t), t ∈ [0, 1] and R1 T T Ji (β) = ψi β γ(·), x0 + 0 pi (t) ẋi (t) − β γ(xi (t))) dt, and compute the derivative of each term with respect to β. Refer to the Appendix B for the completion of the proof. Our objective function gradient with respect to β becomes: n
1X ∇β J(β) = − n i=1
Z 1 pi (t)γ(xi (t))dt + 2λβ, 0
where pi (t) satisfies the adjoint ODE and terminal condition: ∂ ṗi (t) = −pi (t)β T ∂x γ(xi (t)), ∂ ϕ1 (xi (1)). pi (1) = 2[ϕ0 (xi (0) − ϕ1 (xi (1))] ∂x
The algorithm is detailed in Algorithm 1. Algorithm 1 DiffTW using Gradient Descent 1: Initialize β, stopping criteria 2: for iteration = 1 to stopping criteria do 3: Integrate x(t = 1) using β 4: Compute p(1) using x(1) 5: Integrate p(1) backwards in time to obtain p(0) 6: Compute ∇Jβ using quadrature 7: βnew = βold − η∇Jβ 8: end for 9: return βnew
Diffeomorphic Time Warping (DiffTW)
2.5
9
Projection of the Kernel Function
As part of our approach, we require that α vanishes at the boundaries, meaning α(0) = 0 and α(1) = 0. This mirrors the endpoint constraints used in standard DTW. Additionally, it ensures that the endpoints of the signals do not shift independently, allowing the learned dynamics (ẋ = α(x) = β ⊤ γ(x)) to focus on meaningful differences in the body of the signals. To enforce these boundary conditions, we construct a modified RKHS, H0 ⊂ H, that removes any components in the RKHS that are nonzero at the endpoints (x = 0 and x = 1). We do this by projecting out the subspace spanned by the feature vectors evaluated at the boundaries. The resulting kernel ensures that any function drawn from it automatically satisfies the zero-boundary conditions. We define the RKHS H0 ⊂ H as: H0 = {f ∈ H, f (0) = 0 and f (1) = 0}, ⟨f, g⟩H0 = ⟨f, g⟩H ,
g ∈ H0 .
By the reproducing property, f (0) = 0 and f (1) = 0 ⇐⇒ ⟨f, k(·, 0)⟩ = 0 and ⟨f, k(·, 1)⟩ = 0. This implies that H0 = span{k(·, 0), k(·, 1)}⊥ is the orthogonal complement in H. Let k be a positive definite kernel on X with a finite-dimensional feature map γ(x) ∈ Rp from Equation (3), and define Γ = [γ(0), γ(1)] ∈ Rp×2 . Then the projected kernel is: k0 (x, y) = γ(x)⊤ (I − ΓG−1 Γ⊤ )γ(y),
G = Γ⊤ Γ.
Here, A = I − ΓG−1 Γ⊤ is the p × p orthogonal projection matrix onto the subspace of functions vanishing at the boundaries. To see why functions in this space satisfy the boundary constraints, we can verify that A completely annihilates the boundary features:
AΓ = Γ − (ΓG−1 ΓT )Γ = Γ − Γ = 0,
10
Vicky Geneva Haney, Kamel Lahouel, Victor Rielly, and Bruno M. Jedynak
where 0 ∈ p × 2. Since Aγ(0) = 0 and Aγ(1) = 0, it follows that k0 (0, y) = k0 (1, y) = 0 for all y ∈ X , and every α ∈ H0 satisfies α(0) = α(1) = 0. We can now represent the modified kernel as: k0 (x, y) = γ(x)⊤ Aγ(y). Using spectral decomposition, we can write A = U ΛU ⊤ . Here, U ⊤ U = I, U = [u1 , . . . , up−2 ] is the matrix of eigenvectors of A, and Λ is the diagonal matrix of the p−2 positive eigenvalues of A. We can then define the new feature vector for k0 as: γ0 (x) = Λ1/2 U ⊤ γ(x).
(5)
This new feature vector has a dimension of p − 2. Since the projection removes two linearly independent constraints, the feature space of the resulting kernel drops by two dimensions. For details on the derivation, the reproducing property, and the matrix form, see Appendix C.
3
Experiments
We began by evaluating our algorithm on a synthetic toy problem. After confirming its ability to recover known dynamics, we extended the framework to real-world datasets, which we explore later in this section.
3.1
Toy Problem
To test the algorithm, we first constructed a toy problem using random Fourier features. Let the spatial domain be X = [0, 1]. Using random Fourier features, we defined a feature map γ(x) ∈ Rp per Equation (3) and sampled coefficients β0 ∈ Rp to define our initial signal: ϕ0 (x) = β0⊤ γ(x). We then introduced a second set of coefficients β1 ∈ Rp−2 to define the ground-truth vector field governing the temporal dynamics (ẋ = α(x)) and the feature map γ0 (x) ∈ Rp−2 from
Diffeomorphic Time Warping (DiffTW)
11
Equation (5): αtrue (x) = β1⊤ γ0 (x). To ensure smooth, non-trivial temporal dynamics, both β0 and β1 were sampled using p = 2D independent samples from normal distributions scaled proportionally to the feature dimension, where D = 50. For simplicity, we utilized forward Euler integration and evolved ϕ0 through the vector field αtrue to generate the target signal ϕ1 . We trained our algorithm (DiffTW) to recover the vector field using only the start and end signals. Figure 2 visualizes the learned vector field αest alongside the training metrics, and shows the resulting evolved signal ϕ1,est . The hyperparameters used for this toy problem are detailed in Table 1. A standalone interactive demonstration replicating this exact synthetic experiment is available on GitHub.1 After 21 iterations of gradient descent, the algorithm demonstrated fast convergence and successfully reconstructed the underlying dynamics. Encouraged by these findings, we next evaluate our method on real-world datasets, beginning with electrocardiogram (ECG) signals in Section 3.2. Table 1: Parameters of the DiffTW Toy Problem. The algorithm was optimized via gradient descent using the listed parameters, where η is the learning rate and λ is the regularization weight, D number of features, σ for Bandwidth, Iter. for number of iterations of gradient decent, and Step for integration step size. Parameter η Value 0.00675
3.2
λ 1 × 10−4
D 50
σ Iter. Step 0.1 21 1/20
Electrocardiogram Data
After promising results from Section 3.1, we decided to conduct analysis on electrocardiogram (ECG) signals. ECG records the electrical activity of the heart. As shown in Figure 3, a standard cardiac cycle consists of distinct morphological features: the P wave, the sharp QRS complex, and the T wave. For time series modeling, ECGs present signals that contain a complex mixture of smooth variations and sharp peaks. To evaluate our algorithm, we used the MIT-BIH Arrhythmia Database GB and RG (2001), which is one 1
https://github.com/vgeneva/DiffTW/blob/main/demo_toy.py
12
Vicky Geneva Haney, Kamel Lahouel, Victor Rielly, and Bruno M. Jedynak
Original 0 and 1
Amplitude
1 0 1 2 0.0
0.2
0.4
x
0.6
0.8
2 0
1.0
0.0
x(t) Across Iterations
0.4
x
0.6
0.8
1.0
0.8
1.0
_true Iter 0 Iter 5 Iter 10 Iter 15 Iter 20 _est (Final)
0.100 0.075 0.050 0.025
(x)
0.6
0.2
True vs Learned
Warping Path Progression
Iter 0 Iter 5 Iter 10 Iter 15 Iter 20 No Warping
0.8
0.4
0.000 0.025
0.2
0.050 0.075 0.0
0.2
0.4
x(t=0)
0.6
0.8
1.0
Objective Value Progress
0.5 0.4
0.2
0.4
x
0.6
Gradient Norm Progress
5 4
0.3
3
0.2
2
0.1
1
0.0 0.10.0
0.0
Norm Value
Objective Value
4
2
1.0
0.0
1 (Target) offset 0 (Start) 1,est offset
6
Amplitude (with offset)
2
3
The Characteristics: Warping Evolution
1 (Target) 0 (Start)
2.5
5.0
7.5
10.0 12.5 15.0 17.5 20.0
Iteration
00.0
2.5
5.0
7.5
10.0 12.5 15.0 17.5 20.0
Iteration
Figure 2: Top Left: The original function ϕ0 and the target ϕ1 . Top Right: The warping evolution and characteristic curves, illustrating how individual time points migrate from the source ϕ0 to the offset estimated target ϕ1,est across the integration steps. Middle Left: The warping path progression mapping the trajectory dynamics of x(t) across optimization steps, shown for every 5th iteration against the gray dashed identity line. Middle Right: Evolution of the learned velocity field αest , plotted every 5th iteration against a zero reference line. Bottom Left: Objective value across iterations, demonstrating smooth convergence. Bottom Right: Norm of the gradient decreasing during training, indicating optimization stability.
Diffeomorphic Time Warping (DiffTW)
13
R
T
P Q
S
Figure 3: Example of an ECG signal capturing a full cardiac cycle, featuring the P wave, QRS Complex, and T wave. of the most widely used datasets for ECG research. It contains 48 half-hour recordings of two-channel ECG signals from 47 different patients. The dataset includes a variety of heart arrhythmias and normal heartbeats.
3.2.1
An Example: ECG Signal Alignment
Figure 4 presents two processed ECG signals (refer to Section 3.2.2 for processing details), labeled R (right bundle branch block heartbeat) and N (normal heartbeat) along with results of applying our method. For convenience, we denote N as ϕ0 and R as ϕ1 . Optimization of α exhibited decreasing objective functions, settling within 11 iterations. We highlight the changes from ϕ0 (N) evolving to a featured peak in the ϕ1 (R) signal. The hyperparameters used for this ECG problem are detailed in Table 2. The complete code to reproduce this individual ECG signal alignment visualization is available on GitHub.2
Table 2: Parameters of R and N ECG signals for DiffTW. The model was optimized via gradient descent and used the listed parameters where η is the learning rate, λ our regularization term, D number of features, σ for Bandwidth, Iter. for number of iterations of gradient decent, and Step for integration step size. Parameter η Value 1/29 2
λ 1 × 10−4
D 50
σ Iter. Step 0.1 11 1/20
https://github.com/vgeneva/DiffTW/blob/main/demo_ecg.py
14
Vicky Geneva Haney, Kamel Lahouel, Victor Rielly, and Bruno M. Jedynak
Original 0 and 1 5
3 2 1 0
1 (Target) offset 0 (Start) 1,est offset
8
Amplitude (with offset)
Amplitude
4
The Characteristics: Warping Evolution
10
1 (R) 0 (N)
6 4 2 0
1 0.0
0.2
0.4
0.6
x
0.8
1.0
0.0
0.2
Iter 0 Iter 5 Iter 10 No Warping
x(t) Across Iterations
0.8
0.6
x
0.8
Learned Progression
Warping Path Progression
1.0
0.4
1.0 Iter 0 Iter 5 Iter 10 est (Final)
0.08 0.06 0.04
(x)
0.6 0.4
0.02 0.00 0.02
0.2
0.04
0.0
0.0
0.2
0.4
x(t=0)
0.6
0.8
1.0
0.0
0.2
Objective Value Progress
0.4
0.6
x
0.8
1.0
Gradient Norm Progress 12
1.25
10
Objective Value
1.50
Norm Value
1.00 0.75 0.50
8 6 4
0.25
2
0.00 0
2
4
6
Iteration
8
10
00
2
4
6
Iteration
8
10
Figure 4: Top Left: The ECG input signals ϕ0 (N) and the target ϕ1 (R). Top Right: The warping evolution and characteristic curves. Middle Left: Warping path progression mapping the trajectory dynamics of x(t) across 11 optimization iterations. Middle Right: Evolution of the learned velocity field αest over 11 iterations. Bottom Left: Objective value across iterations, showing decreasing values. Bottom Right: Norm of the gradient decreasing during training.
Diffeomorphic Time Warping (DiffTW)
3.2.2
15
ECG Pre-processing and Classification Procedure
To establish a rigorous evaluation framework, we mirror the preprocessing and classification ideas proposed by Raghavendra et al. (2011). The authors proposed a DTW-based approach for arrhythmic ECG beat classification and applied a series of preprocessing steps. Following their protocol to ensure morphological consistency across signals, we first removed low-frequency baseline wanderings using a Butterworth high-pass filter at 0.5 Hz, followed by a low-pass filter at 40 Hz to limit the signal bandwidth. Next, Rpeaks were identified, based on cardiologist annotations, to locate individual heartbeats. Full cardiac cycles were then extracted by taking 100 samples before and 155 samples after each R-peak, which captures the P-wave, QRS complex, and T-wave. This yielded 256 sample points per beat. Finally, each beat was z-normalized so that the average of the signal was zero, keeping the beats comparable. This was slightly different from the authors’ method, which scaled the amplitudes to the range 0–1. To evaluate classification performance, we selected the same eight ECG heartbeat classes as in Raghavendra et al. (2011). The specific classes, along with the distribution of the training and testing sets, are summarized in Table 3. To establish a benchmark Table 3: ECG Labels. Normal (N), Premature ventricular contraction (V), Atrial premature beat (A), Right bundle branch block beat (R), Left bundle branch block beat (L), Ventricular escape beat (E), Junctional (nodal) premature beat (J), and Junctional (Nodal) escape beat (j) for classification. Class label Test Signals Train Signals
N 1000 10
V 982 10
A 1000 10
R 500 10
L 500 10
E 95 10
J 40 10
j 202 10
for our DiffTW algorithm, we reproduced the standard DTW baseline results. While the authors in Raghavendra et al. (2011) utilized a class-wise classification approach and explored sub-sampling to improve detection speed, we opted for a more direct comparison. We performed our analysis without sub-sampling and adopted a standard k-nearest neighbor (k = 1) classification strategy. Specifically, distances are computed between a given test beat and all training beats across the arrhythmia classes, and the test beat is assigned the label of the training beat that yields the minimum distance. To ensure a
16
Vicky Geneva Haney, Kamel Lahouel, Victor Rielly, and Bruno M. Jedynak
straightforward comparison, the distance metrics for both frameworks are derived from their respective alignment processes. For traditional DTW, the distance is defined as the minimum cumulative alignment cost calculated along an optimal discrete warping path. Let ϕ0 and ϕ1 represent two discrete ECG signals. Standard DTW calculates the distance by finding an optimal discrete warping path π that minimizes the accumulated Euclidean distance between the aligned time points:
DDTW (ϕ0 , ϕ1 ) = min π
K X
∥ϕ0 (ik ) − ϕ1 (jk )∥2
(6)
k=1
where K is the total length of the warping path, and (ik , jk ) represents the pair of mapped indices at step k. This minimization is solved via dynamic programming, subject to standard boundary, monotonicity, and step-size continuity constraints Chen et al. (2018). Conversely, DiffTW defines its similarity measure directly through its continuous optimization framework. The distance corresponds to the minimum value achieved by the objective function J(β) from Equation (2) at convergence or stopping criteria:
DDiffTW (ϕ0 , ϕ1 ) = min J(β) β
(7)
This setup provides a direct, distance-based framework to evaluate the accuracy of DiffTW against DTW. Applying these distance metrics within our 1-nearest neighbor framework yielded strong classification performance on the ECG dataset. Specifically, DiffTW outperformed the DTW baseline, achieving an overall accuracy of 98.80% compared to DTW’s 98.12%. Having established the framework’s effectiveness on both synthetic data and specific ECG morphological alignments, we next extended our analysis to broader datasets. To thoroughly assess DiffTW’s performance relative to DTW, we evaluate it across datasets that vary significantly in time-series length and structural complexity, grouping our ECG classification task with the University of California Riverside Time Series Classification Archive, refer to Section 3.3. The standardized hyperparameter tuning and optimization strategy applied across all of these classification tasks is detailed in Section 3.4, and the full benchmark classification results for the ECG dataset,
Diffeomorphic Time Warping (DiffTW)
17
alongside the UCR archive, are presented in Section 4.
3.3
University of California Riverside Time Series Classification Archive
To benchmark our algorithm at scale, we utilized the University of California Riverside (UCR) Time Series Classification Archive Dau et al. (2019), a standard repository containing 128 datasets. Of these, 85 datasets were included in our experimental analysis. We excluded datasets that contained missing values, inconsistent sequence lengths, multivariate data formats, or an insufficient number of training/testing samples required for proper hyperparameter tuning (as detailed in Section 3.4). Of the 85 datasets selected for testing, 70 were provided already z-normalized, 9 required manual z-normalization (annotated in Table 5), and 6 contained exceptionally long time series (length > 1000) or excessively large sample sizes. These were subsampled before applying both DiffTW and DTW to maintain computational feasibility (annotated in Table 5). Detailed results and specific processing annotations for all 85 datasets are provided in Table 5. For additional metadata, including study contexts and dataset origins, we refer readers to Middlehurst et al. (2024).
3.4
Hyperparameter tuning and Broader Testing
While our initial ECG experiments utilized fixed values for λ and bandwidth (Table 2), evaluating the broader UCR archive required a systematic tuning and training strategy. We performed a grid search over 100 log-spaced values for the regularization term λ ∈ [10−2 , 103 ], and 10 uniformly spaced bandwidth values between 1/30 and 1/5. To tune these parameters, we constructed representative validation pairs for each class using only the training set. Specifically, we averaged two training signals to form the initial signal ϕ0 , and averaged a second pair of training signals to form the target signal ϕ1 . By averaging pairs of training signals, we generate stable, representative signals that allow generalization rather than overfitting to an individual possibly noisy signal, a technique supported by sequence averaging frameworks from Petitjean et al. (2011). We then
18
Vicky Geneva Haney, Kamel Lahouel, Victor Rielly, and Bruno M. Jedynak
evaluated each (λ, bandwidth) pair using an index-based 2-fold cross-validation strategy. The temporal indices of ϕ0 and ϕ1 were split equally into even and odd sets. We optimized the weights β for 20 iterations on one subset and computed the objective value on the held-out subset, then swapped the sets. The hyperparameters that yielded the minimum average objective value across both folds were selected to process the full training and test datasets. Finally, we replaced standard gradient descent with the Adaptive Moment Estimation (Adam) optimizer Kingma and Ba (2015) and conducted all the experiments for up to 10 iterations utilizing GPU acceleration. Detailed hardware specifications and execution environments are provided in Appendix A. Refer to Table 4 for summary of parameters. Table 4: Parameters used for training on GPU acceleration where GS stands for grid search, η is the learning rate, λ our regularization term, D number of features, σ for Bandwidth, Iter. for number of iterations of gradient decent, and Step for integration step size. Parameter η Value 1/29
3.5
λ GS
D 50
σ Iter. GS 10
Step 1/20
Computational Complexity
Importantly, Algorithm 1 computes the transport by integrating the ODE from t = 0 to t = 1. This approach yields a computational complexity of O(n · Ksteps · iter), where n is the sequence length (i.e., the number of transported points), Ksteps is the number of integration steps, and iter is the number of optimization steps. In contrast, traditional DTW requires explicitly constructing a pairwise alignment cost matrix, resulting in a quadratic complexity of O(n2 ) with respect to sequence length. By reframing sequence alignment as a continuous dynamical system, our method replaces DTW’s quadratic spatial dependence with a linear dependence on n. While DTW relies on exhaustive discrete path searching to handle non-linear time distortions, DiffTW captures these distortions efficiently through smooth vector fields. This exchange of exact dynamic programming for continuous integration drastically reduces computational overhead, making our method highly scalable for long sequences and large datasets.
Diffeomorphic Time Warping (DiffTW)
4
19
Results
Table 5 reports the overall classification accuracy of DiffTW against the standard DTW baseline across 86 datasets (1 ECG dataset and 85 UCR datasets). Accuracy is defined as the number of correctly classified test signals divided by the total number of test signals. For the 1-Nearest Neighbor (1-NN) classification scheme, standard DTW utilizes the accumulated Euclidean distance along the optimal alignment path as discussed in Equation (6). To ensure a fair comparison, DiffTW classification was based on the optimized objective value between pairs of signals as discussed in Equation (7). Because our objective function measures the sum of squared differences between ϕ0 (xi (0)) and the transformed ϕ1 (xi (1)), it serves as an analogous, continuous counterpart to the discrete DTW distance. Table 5: The classification accuracies for both approaches across all 86 datasets are reported below. This accuracy is the overall accuracy for all the test signals, i.e. total correctly labeled divided by total signals. The bold results represent the higher score.
Datasets
Length
Label
DiffTW
DTW
ACSFOne
1460
10
0.0900
0.8000
Adiac
176
37
0.2046
0.0563
ArrowHead
251
3
0.7943
0.6686
Beef
470
5
0.6000
0.5667
BeetleFly
512
2
0.7000
0.7000
BirdChicken
512
2
0.8500
0.7500
BME
128
3
1.0
0.8867
Car
577
4
0.7833
0.6000
Crop
46
24
0.7188
0.6744
CBF
128
3
0.9700
0.9267
ChlorineConcentration
166
3
0.5029
0.3698
CinCECGTorso
1639
4
0.8688
0.4638
Coffee
286
2
1.0
0.9643
Computers
720
2
0.5320
0.6560
Earthquakes
512
2
0.6547
0.6619
Continued on next page...
20
Vicky Geneva Haney, Kamel Lahouel, Victor Rielly, and Bruno M. Jedynak
Table 5 – Continued from previous page Datasets
Length
Labels
DiffTW
DTW
ECG
256
8
0.9880
0.9812
ECG200
96
2
0.8600
0.7900
ECG5000
140
5
0.9319
0.9228
ECGFiveDays
136
2
0.7666
0.7514
ElectricDevices1
483
7
0.5471
0.2043
EOGHorizontalSignal
1250
12
0.4669
0.4309
EOGVerticalSignal
1250
12
0.4558
0.4337
Ethanollevel1
8752
4
0.2800
0.2800
FaceAll
560
14
0.8041
0.7509
FacesUCR
131
14
0.9093
0.0800
Fish
463
7
0.8857
0.7029
FordA1
2503
2
0.4644
0.5856
FordB1
2503
2
0.4580
0.6111
FreezerRegularTrain
301
2
0.7663
0.8730
FreezerSmallTrain
301
2
0.7193
0.6965
GunPoint
150
2
0.9733
0.8933
GunPointAgeSpan2
150
2
0.9842
0.9146
GunPointMaleVersusFemale2
150
2
0.9937
0.9968
GunPointOldVersusYoung2
150
2
0.9651
0.8159
Ham
431
2
0.6095
0.6095
HandOutlines1
677 3
2
0.8703
0.5649
Haptics
1092
5
0.4448
0.2662
Herring
512
2
0.5938
0.5469
HouseTwenty2
301
2
0.6555
0.9244
InlineSkate
1882
7
0.3964
0.2927
InsectWingbeatSound
256
11
0.5606
0.4333
InsectEPGRegularTrain2
601
3
0.7671
0.7912
ItalyPowerDemand
24
2
0.9631
0.9184
LargeKitchenAppliances
720
3
0.5653
0.8053
Lightning2
637
2
0.7049
0.7049
Lightning7
319
2
0.7123
0.2329
Mallat
1024
8
0.8827
0.2217
Continued on next page...
Diffeomorphic Time Warping (DiffTW)
21
Table 5 – Continued from previous page Datasets
Length
Labels
DiffTW
DTW
Meat
448
3
0.5167
0.9500
MedicalImages
99
10
0.7605
0.0895
MixedShapesRegularTrain
1024
5
0.8973
0.8433
MixedShapesSmallTrain
1024
5
0.8243
0.7790
84
2
0.8850
0.8970
NonInvasiveFetalECGThorax11
3753
42
0.7903
0.0763
NonInvasiveFetalECGThorax21
3753
42
0.8758
0.0830
OliveOil
570
4
0.0667
0.2000
OSULeaf
427
6
0.6364
0.4215
PhalangesOutlinesCorrect
80
3
0.7681
0.5536
Plane
144
7
0.9905
0.2571
PowerCons
144
2
0.9778
0.8778
RefrigerationDevices
720
3
0.4507
0.4613
Rock2
2844
4
0.7800
0.7000
ScreenType
720
3
0.3733
0.4107
SemgHandGenderCh22
1500
2
0.5950
0.6017
SemgHandMovementCh22
1500
6
0.2244
0.3222
SemgHandSubjectCh22
1500
5
0.2889
0.4467
ShapeletSim
500
2
0.5056
0.7167
ShapesAll
512
60
0.8183
0.7633
SmallKitchenAppliances
720
3
0.3787
0.6800
SmoothSubspace
15
3
0.9667
0.8533
SonyAIBORobotSurface1
70
2
0.7787
0.6839
SonyAIBORobotSurface2
65
2
0.8604
0.6590
StarLightCurves1
5123
3
0.8475
0.2580
Strawberry
235
2
0.8865
0.7027
SwedishLeaf
128
15
0.8624
0.2400
SyntheticControl
60
6
0.9700
0.9467
ToeSegmentation1
277
2
0.7807
0.7982
ToeSegmentation2
343
2
0.8846
0.8462
Trace
275
4
0.9900
0.7400
TwoLeadECG
82
2
0.8560
0.9271
MoteStrain
Continued on next page...
22
Vicky Geneva Haney, Kamel Lahouel, Victor Rielly, and Bruno M. Jedynak
Table 5 – Continued from previous page Datasets
4.1
Length
Labels
DiffTW
DTW
TwoPatterns
128
4
0.9982
0.9878
UMD
150
3
0.9861
0.8472
Wafer
152
2
0.8960
0.2052
Wine
234
2
0.4815
0.5926
Worms
900
5
0.4805
0.4675
WormsTwoClass
900
2
0.6753
0.5584
Yoga
426
2
0.8640
0.7853
Discussion
We evaluated DiffTW through a comprehensive benchmarking study of 86 time-series datasets (85 from the UCR archive, plus a specialized ECG heartbeat dataset). To ensure a rigorous comparison, we excluded datasets with missing values, inconsistent lengths, or multivariate components. Classification was performed using a k-nearest neighbor (k = 1) approach, utilizing the respective alignment costs of DiffTW, Equation (7), and standard DTW, Equation (6), as distance metrics. In head-to-head accuracy, DiffTW achieved higher scores on 60 datasets, while DTW performed better on 22, with 4 ties. A two-sided Wilcoxon signed-rank test yielded a p-value of approximately 0.0001, confirming that the performance gains offered by DiffTW are statistically significant across the archive. As illustrated in the boxplots in Figure 5, DiffTW achieves a higher median accuracy and more consistent high-tier performance compared to DTW. The method’s efficacy is particularly evident in tasks requiring precise morphological alignment. In the ECG study, DiffTW successfully aligned key cardiac peaks within 10 iterations by learning an optimal vector field, α. This highlights a core strength of the continuous framework: its ability to capture subtle morphological differences, such as those between arrhythmia 1
Data was subsampled due to large series length or excessively large number of training and testing samples. 2 z-normalized prior to placing in algorithms. 3 Dataset ElectricDevices was subsampled from 96 to 48. Datasets FordA and FordB were subsampled from 500 to 250. Datasets NonInvasiveFetalECGThorax1 and NonInvasiveFetalECGThorax2 were sub-sampled from 750 to 375. Dataset StarLightCurves was subsampled from 1024 to 512. Dataset HandOutlines subsampled from 2709 to 677. Dataset Ethanollevel subsampled from 1751 to 875.
Diffeomorphic Time Warping (DiffTW)
23
Accuracy Comparison: DiffTW vs DTW 1.0
Accuracy
0.8 0.6 0.4 0.2 0.0 DiffTW
DTW
Method
Figure 5: Accuracy comparison across 86 datasets. DiffTW achieves a higher median accuracy and consistently higher classification performance. types, that rigid or purely discrete methods often overlook. This advantage is further demonstrated by the synthetic BME dataset (Figure 6). BME consists of three classes characterized by the timing of a small positive bell: one arising at the initial period (Begin), one with no bell (Middle), and one at the final period (End). DiffTW perfectly captured these localized temporal features, achieving 100% accuracy compared to DTW’s 87%. BME, with z-normalization Class 1
1 0 1
0
20
40
60
80
100
120
140
80
100
120
140
80
100
120
140
Class 2
1.0 0.5 0.0
0
20
40
60
Class 3
1 0 1 0
20
40
60
Figure 6: Performance on the BME dataset, highlighting DiffTW’s ability to precisely align smooth, localized features like the timing of a specific bell curve. Although DiffTW demonstrates strong overall performance, our evaluation reveals specific signal characteristics where discrete DTW retains a distinct advantage. A notable example is the ACSFOne dataset (Figure 7), where DTW significantly outperformed
24
Vicky Geneva Haney, Kamel Lahouel, Victor Rielly, and Bruno M. Jedynak
DiffTW (80% vs. 10% accuracy). ACSFOne signatures are characterized by long idle periods punctuated by sparse, high-energy bursts of power consumption. In such cases, the rigid, step-wise nature of discrete DTW is more robust; the continuous deformation of DiffTW tends to over-smooth or struggle to map isolated, high-magnitude spikes across long stretches of inactivity. Furthermore, in datasets like Ham (Figure 8), which involve ACSFone, with z-normalization 1 0 1 1 0 1
Class 0 0
200
400
600
0
200
400
600
0
200
400
600
0
200
400
600
0
200
400
600
0
200
400
600
0
200
400
600
0
200
400
600
0
200
400
600
0
200
400
600
800
1000
1200
1400
800
1000
1200
1400
800
1000
1200
1400
800
1000
1200
1400
800
1000
1200
1400
800
1000
1200
1400
800
1000
1200
1400
800
1000
1200
1400
800
1000
1200
1400
800
1000
1200
1400
Class 1
Class 2
1 0
10 0 1 0 1
Class 3
Class 4
Class 5
1 0
Class 6
1 0 10 0
2.5 0.0
Class 7
Class 8
Class 9
1 0
Figure 7: ACSFOne dataset signatures, characterized by sparse high-energy bursts that favor discrete DTW alignment. spectrographic measurements for food quality assurance, both methods achieved an equal accuracy of 61%. This parity suggests that the primary challenge lies in the complex, overlapping nature of the dataset itself, proving that the moderate accuracy of 61% is a reflection of the data’s inherent difficulty rather than a flaw in the DiffTW framework. Ultimately, while the 22 datasets where DTW maintained an advantage show that extreme discrete symmetries or sparse bursts still favor traditional warping, the broad success of DiffTW indicates that its continuous vector field framework provides greater
Diffeomorphic Time Warping (DiffTW)
25
Ham, with z-normalization Class 1
5 0 0
50
100
150
200
250
300
350
400
450
300
350
400
450
Class 2 5 0 0
50
100
150
200
250
Figure 8: Ham dataset spectrographic measurements, where both methods achieved identical classification accuracy of 61%. flexibility and accuracy for handling complex temporal distortions across diverse domains.
5
Conclusion
In this work, we introduced DiffTW, a continuous temporal alignment framework driven by the integration of smooth vector fields. Extensive benchmarking across 86 datasets demonstrates that DiffTW matched or outperformed standard discrete DTW on the vast majority of classification tasks, proving particularly effective at capturing and precisely aligning localized morphological features. However, our evaluation also identified specific limitations of the continuous approach. In datasets where DTW outperformed DiffTW, the time series were frequently characterized by high-frequency oscillations. These sharp, rapid fluctuations negatively impact the gradient calculations required to smoothly learn the optimal vector field α. Additionally, on datasets where the two methods tied or exhibited remarkably similar performance, the underlying classes were often visually indistinguishable. This parity suggests that inherent dataset noise or a lack of clear distinguishing features, rather than the alignment mechanism itself, bounded the maximum achievable accuracy in those instances. Future work will explore hybrid approaches to accelerate convergence, such as utilizing discrete DTW to initialize the continuous optimization. By deriving an initial
26
Vicky Geneva Haney, Kamel Lahouel, Victor Rielly, and Bruno M. Jedynak
parameter β from the discrete alignment, we hypothesize that this initialization will significantly reduce training time while preserving the benefits of a continuous alignment path. Furthermore, given that DiffTW already achieves linear theoretical complexity with respect to sequence length, subsequent research will focus on empirical runtime analysis and implementation-level optimizations to further accelerate the framework. Overall, DiffTW offers a highly flexible alternative to traditional dynamic programming that yields statistically significant improvements for time-series classification, provided the underlying signals do not contain severe high-frequency noise that disrupts continuous, gradient-based optimization.
6
Code
The computational framework and algorithm proposed in this paper were implemented in Python. To ensure fast and stable continuous trajectory generation, the model architecture, forward Euler integrations, and gradient descent optimization were built utilizing tensorflow and tensorflow-probability. The complete source code and reproducibility instructions are provided in the Appendix A.
7
Funding
This work was supported in part by the National Science Foundation Research Training Grant (DMS-2136228) and by the National Institute of Health RO1AG021155, R01EY032284, and R01AG027161.
Diffeomorphic Time Warping (DiffTW)
27
References Beg MF, Miller MI, Trouvé A, Younes L (2005). Computing large deformation metric mappings via geodesic flows of diffeomorphisms. International Journal of Computer Vision, 61(2): 139–157. Buchin K, Nusser A, Wong S (2023). Computing continuous dynamic time warping of time series in polynomial time. Chen RTQ, Rubanova Y, Bettencourt J, Duvenaud D (2018). Neural ordinary differential equations. In: Advances in Neural Information Processing Systems, volume 31. Curran Associates, Inc. Cuturi M, Blondel M (2017). Soft-DTW: a differentiable loss function for time-series. In: Proceedings of the 34th International Conference on Machine Learning (D Precup, YW Teh, eds.), volume 70 of Proceedings of Machine Learning Research, 894–903. PMLR. Dau HA, Bagnall A, Kamgar K, Yeh CCM, Yan Zhu SG, Ratanamahatana CA, et al. (2019). The UCR Time Series Archive. GB M, RG M (2001). The Impact of the MIT-BIH Arrhythmia Database. Giorgino T (2009). Computing and visualizing dynamic time warping alignments in r: The dtw package. Version 1.5.3. Kingma DP, Ba J (2015). Adam: A method for stochastic optimization. In: International Conference on Learning Representations (ICLR). San Diego, CA, USA. Liu Y, Zhang YA, Zeng M, Zhao J (2024a). A novel distance measure based on dynamic time warping to improve time series classification. Information Sciences, 656: 119921. Liu Y, Zhang YA, Zeng M, Zhao J (2024b). A novel distance measure based on dynamic time warping to improve time series classification. Information Sciences, 656: 119921. McOwen RC (2003). Partial Differential Equations: Methods and Applications. Prentice Hall, Upper Saddle River, NJ. Middlehurst M, Schafer P, Bagnall A (2024). Bake off redux: a review and experimental evaluation of recent time series classification algorithms. Data Mining and Knowledge Discovery.
28
Vicky Geneva Haney, Kamel Lahouel, Victor Rielly, and Bruno M. Jedynak
Pereverzyev S (2022). An Introduction to Artificial Intelligence Based on Reproducing Kernel Hilbert Spaces. Birkhäuser, Switzerland. Petitjean F, Ketterlin A, Gançarski P (2011). A global averaging method for dynamic time warping, with applications to clustering. Pattern Recognition, 44(3): 678–693. Raghavendra BS, Bera D, Bopardikar AS, Narayanan R (2011). Cardiac arrhythmia detection using dynamic time warping of ecg beats in e-healthcare systems. Technical report, Samsung Advanced Institute of Technology Lab. Rahimi A, Recht B (2007). Random features for large-scale kernel machines. In: Advances in neural information processing systems, volume 20.
Diffeomorphic Time Warping (DiffTW)
A
29
Hardware Environment and Code Execution
To efficiently process the large matrix operations and parallelize the optimization tasks across the 86 evaluated datasets, experiments were executed on a high-performance computing cluster utilizing NVIDIA A40 GPUs (48GB VRAM) running CUDA 12.8. All Python scripts, including demo_toy.py, demo_ecg.py and difftw.py, as well as the accompanying requirements.txt and datasets, are hosted on our public GitHub repository: https://github.com/vgeneva/DiffTW. The repo provides example dataframes to use from the UCR archive. To utilize the algorithm on two timeseries, first install the required packages using pip: pip install -r requirement.txt
B
Derivation of the Objective Function
Proof of the Objective Function. Consider the functional: n 1X T T ψi β γ(·), x0 , ψ β γ(·), x0 = n i=1 with
T
ψi β γ(·), x0
h i2 = ϕ1 (xi (1)) − ϕ0 (xi (0)) ,
for i = {1, . . . , n} under the constraints defined by the ẋ(t) = α(x) = β T γ(x(t)) and x(0) = x0 . Then for each i: ẋi (t) = β T γ(xi (t)), xi (0) = xi (0). Next, we augment our objective function with Lagragian methods to replace this constrained optimization problem and define the Lagrange multipliers pi (t), t ∈ [0, 1] and R1 T T Ji (β) = ψi β γ(·), x0 + 0 pi (t) ẋi (t) − β γ(xi (t))) dt,
(8)
30
Vicky Geneva Haney, Kamel Lahouel, Victor Rielly, and Bruno M. Jedynak
and compute the derivative of each term with respect to β. We begin with the first summand of 8: ∂ ψ ∂β i
∂ T β γ(·), x0 = ∂β [ϕ1 (xi (1)) − ϕ0 (xi (0))]2 = 2 [ϕ1 (xi (1)) − ϕ0 (xi (0))] i h ∂ ∂ ∂ ∂ ϕ1 (xi (1)) ∂β xi (1) − ∂x ϕ0 (xi (0)) ∂β xi (0) · ∂x h i ∂ ∂ = 2 [ϕ1 (xi (1)) − ϕ0 (xi (0))] ∂x ϕ1 (xi (1)) ∂β xi (1) ,
∂ using the fact that xi (0) are constant, thus ∂β xi (0) = 0 for each i. Next,
" ∂ ∂β
R1 0
pi (t)ẋi (t)dt
∂ = ∂β
#
1
R1
pi (t)xi (t) − 0 ṗi (t)xi (t)dt 0
h i R1 pi (1)xi (1) − pi (0)xi (0) − 0 ṗi (t)xi (t)dt R1 ∂ ∂ ∂ = pi (1) ∂β xi (1) − pi (0) ∂β xi (0) − 0 ṗi (t) ∂β xi (t)dt R 1 ∂ ∂ xi (1) − 0 ṗi (t) ∂β xi (t)dt, = pi (1) ∂β ∂ = ∂β
∂ where again we used ∂β xi (0) = 0 and integration by parts in the second line. Lastly,
∂ − ∂β
R1 ∂ T T p (t)β γ(x (t))dt = − p (t) β γ(x (t)) dt i i i i ∂β 0 0 R1 ∂ ∂ γ(xi (t)) ∂β xi (t) dt = − 0 pi (t) γ(xi (t)) + β T ∂x R1 R1 ∂ ∂ = − 0 pi (t)γ(xi (t))dt − 0 pi (t)β T ∂x γ(xi (t)) ∂β xi (t)dt.
R1
Finally, we address the regulation term:
∂ ∂ ∂ ∂ λ ∂β ∥β∥2 = λ( ∂β ⟨β, β⟩) = λ(⟨ ∂β β, β⟩ + ⟨ ∂β β, β⟩) = λ(⟨I, β⟩ + ⟨I, β⟩) = 2λβ.
Combining the results above, we obtain: ∂ J (β) ∂β i
h i ∂ ∂ = 2 [ϕ1 (xi (1)) − ϕ0 (xi (0))] ∂x ϕ1 (xi (1)) ∂β xi (1) R1 ∂ ∂ +pi (1) ∂β xi (1) − 0 ṗi (t) ∂β xi (t)dt R1 − 0 pi (t)γ(xi (t))dt R1 ∂ ∂ γ(xi (t)) ∂β xi (t)dt. − 0 pi (t)β T ∂x
Diffeomorphic Time Warping (DiffTW)
31
We choose: ṗi (t) = −pi (t)β T
∂ γ(xi (t)). ∂x
then we cancel out two terms and obtain: ∂ J (β) ∂β i
= 2 [ϕ1 (xi (1)) − ϕ0 (xi (0))] h i ∂ ∂ · ∂x ϕ1 (xi (1)) ∂β xi (1) R1 ∂ xi (1) − 0 pi (t)γ(xi (t))dt. +pi (1) ∂β
Let pi (1) = 2 (ϕ0 (xi (0)) − ϕ1 (xi (1)))
∂ ϕ1 (xi (1)), ∂x
then our derivative is represented as ∂ Ji (β) = − ∂β
C
Z 1 pi (t)γ(xi (t))dt. 0
Projection Kernel Derivation
Let X ̸= ∅ and k be a positive definite kernel function on X and H be the RKHS of k. We will construct an RKHS H0 ⊂ H such that H0 = {f ∈ H, f (0) = 0 and f (1) = 0}, ⟨f, g⟩H0 = ⟨f, g⟩H . Note that since f ∈ H, we can leverage the reproducing property that shows: f (0) = 0 and f (1) = 0 ⇐⇒ ⟨f, k(·, 0)⟩ = 0 and ⟨f, k(·, 1)⟩ = 0. This implies that H0 is the orthogonal complement in H of V = span{k(·, 0), k(·, 1)}. Then for y ∈ X the orthogonal projection πV k(·, y) of the function k(·, y) on V is char-
32
Vicky Geneva Haney, Kamel Lahouel, Victor Rielly, and Bruno M. Jedynak
acterized by πV k(·, y) ∈ V,
(9)
k(·, y) − πV k(·, y), k(·, u) = 0 for u = 0, 1. We can represent the projection of πV k(·, y) as a linearly combination of the kernel functions of V (since they represent a basis for V ):
πV k(·, y) =
1 X
µv (y)k(·, v) = µT k(·),
v=0
where k(·) = (k(·, 0), k(·, 1))T , µv (y) ∈ R and µ = (µ0 (y), µ1 (y))T . Substituting πV k(·, y) into Equation (9) and using the reproducing property: = k(·, y) − πV k(·, y), k(·, u) P = k(·, y) − 1v=0 µv k(·, v), k(·, u) P = k(·, y), k(·, u) − 1v=0 µv k(·, v), k(·, u) P = k(u, y) − 1v=0 µv k(u, v) P ⇒ k(u, y) = 1v=0 µv k(u, v),
0
for u = 0, 1. We can represent k(u, y) in matrix form for each u = 0, 1. Let k(y) = (k(0, y), k(1, y))T , µ = (µ0 (y), µ1 (y))T , then we can represent k(y) as k(0, 0) k(0, 1) k(y) = Gµ with G = . k(1, 0) k(1, 1) Then if we assume G is invertible we obtain µ = G−1 k(y) and πV k(·, y) = µT k(·) = k(y)T G−1 k(·),
Diffeomorphic Time Warping (DiffTW)
33
since (G−1 )T = (GT )−1 = G−1 . And since k(·, y) − πV k(·, y) ∈ H0 ⇒ πH0 k(·, y) = k(·, y) − πV k(·, y). We can define the kernel k0 as k0 (x, y) = ⟨πH0 k(·, y), πH0 k(·, x)⟩ = ⟨k(·, y) − πV k(·, y), k(·, x) − πV k(·, x)⟩ = k(y, x) − ⟨k(·, y), πV k(·, x)⟩ −⟨πV k(·, y), k(·, x)⟩ + ⟨πV k(·, y), πV k(·, x)⟩ = k(y, x) − ⟨k(·, y), k(x)T G−1 k(·)⟩ −⟨k(y)T G−1 k(·), k(·, x)⟩ +⟨k(y)T G−1 k(·), k(x)T G−1 k(·)⟩.
Let k(x)T G−1 = (α1 (x) α2 (x) = αT (x) ∈ R1×2 , then
⟨k(·, y), k(x)T G−1 k(·)⟩ = ⟨k(·, y), αT (x)k(·)⟩ = ⟨k(·, y), αT (x)(k(·, 0), k(·, 1))T ⟩ = ⟨k(·, y), α1 (x)k(·, 0) + α2 (x)k(·, 1)⟩ = ⟨k(·, y), α1 (x)k(·, 0)⟩ + ⟨k(·, y), α2 (x)k(·, 1)⟩ = α1 (x)k(0, y) + α2 (x)k(1, y) = k(y)T α(x) = k(y)T (k(x)T G−1 )T = k(y)T G−1 k(x).
T
−1
Similarly, let k(y) G
T
= β1 (y) β2 (y) = β (y):
⟨k(y)T G−1 k(·), k(·, x)⟩ = k(x)T G−1 k(y) = k(y)T G−1 k(x),
(10)
34
Vicky Geneva Haney, Kamel Lahouel, Victor Rielly, and Bruno M. Jedynak
where the last equality is a property of the inner product. Then finally, T
⟨k(y)T G−1 k(·), k(x)T G−1 k(·)⟩ = ⟨β (y)k(·), αT (x)k(·)⟩ = ⟨β1 (y)k(·, 0) + β2 (y)k(·, 1), α1 (x)k(·, 0) + α2 (x)k(·, 1)⟩ = ⟨β1 (y)k(·, 0), α1 (x)k(·, 0)⟩ + ⟨β1 (y)k(·, 0), α2 (x)k(·, 1)⟩ +⟨β2 (y)k(·, 0), α1 k(·, 1)⟩ + ⟨β2 (y)k(·, 0), α2 (x)k(·, 1)⟩ = β1 (y)α1 (x)k(0, 0) + β1 (y)α2 (x)k(0, 1) +β2 (y)α1 (x)k(1, 0) + β2 (y)α2 k(x)(1, 1) T k(0, 0) k(1, 0) = β1 β2 α1 α2 k(1, 0) k(1, 1) T
= β Gα = k(y)T G−1 G(k(x)T G−1 )T = k(y)T G−1 GG−1 k(x) = k(y)T G−1 k(x).
So Equation (10) becomes: k0 (x, y) = k(x, y) − k(y)T G−1 k(x) T
(11)
−1
= k(x, y) − k(x) G k(y). Proposition 1. The pd kernel of H0 is k0 Proof. We want to show that k0 ∈ H0 and the reproducing property holds. From Equation (11), we obtain the function
k0 (·, y) = k(·, y) − k(·)T G−1 k(y)
which is a linear combination of k(·, y), k(·, 0) and k(·, 1), all of which belong to H, our RKHS. Moreover: k0 (0, y) k(0, y) k(0, 0) k(0, 1) −1 k(0, y) = − G = 0, k0 (1, y) k(1, y) k(1, 0) k(1, 1) k(1, y) since GG−1 = I. Therefore, k0 (·, y) ∈ H0 . Next, we verify the reproducing property
Diffeomorphic Time Warping (DiffTW)
35
holds. Take any f ∈ H0 ⟨f, k0 (·, y)⟩H
= ⟨f, k(·, y) − k(·)T G−1 k(y)⟩H = ⟨f, k(·, y)⟩ − ⟨f, k(·)T G−1 k(y)⟩H = ⟨f, k(·, y)⟩ − ⟨f, k(y)T G−1 k(·)⟩H T
= f (y) − ⟨f, β k(·)⟩H = f (y) − ⟨f, β1 (y)k(·, 0) − β2 (y)k(·, 1)⟩H = f (y) − β1 (y)⟨f, k(·, 0)⟩H − β2 (y)⟨f, k(·, 1)⟩H = f (y) − β1 (y)f (0) − β2 (y)f (1) T f (0) = f (y) − β f (1) f (0) = f (y) − k(y)T G−1 f (1) = f (y),
since f ∈ H0 , then f (0) = f (1) = 0. Therefore, ⟨·, ·⟩ ∈ H is the inner product in H0 . We observe for {x1 , x2 , . . . , xn }, for xi ∈ X , then
k0 (xi , xj ) = k(xi , xj ) − k(xi )T G−1 k(xj ),
where the elements of the usual Gram Matrix, K is represented by [K]ij = k(xi , xj ). Let us represent k(xi )T G−1 k(xj ) in matrix form. Notate
T
k(x1 , 0) k(x1 , 1) k (x1 ) . . .. . . L= . . = . , T k (xn ) k(xn , 0) k(xn , 1)
36
Vicky Geneva Haney, Kamel Lahouel, Victor Rielly, and Bruno M. Jedynak
then
T
T
k (x1 ) α (x1 ) . −1 . .. G . k(x1 ) · · · k(xn ) = . k(x1 ) · · · k(xn ) T T k (xn ) α (xn ) . T T α (x1 )k(x1 ) · · · α (x1 )k(xn ) .. .. ... = . . T T α (xn )k(x1 ) · · · α (xn )k(xn ) Hence, K0 = K − LG−1 LT In the case of an explicit kernel, we can write k(x, y) = γ(x)T γ(y), where γ : X → Rp is a feature map. We compute:
T
−1
T
k0 (x, y) = γ(x) γ(y) − γ(x) γ(0), γ(x) γ(0) G γ(y) γ(0), γ(y) γ(1) T T −1 T T T T = γ(x) γ(y) − γ(x) γ(0), γ(x) γ(1) G γ(0) γ(y), γ(1) γ(y) T T T −1 = γ(x) γ(y) − γ(x) γ(0), γ(1) G γ(0), γ(1) γ(y) " # T
T
T
T
= γ(x)
T
−1
I − γ(0), γ(1) G
γ(0), γ(1)
Let us define Γ = γ(0) γ(1) ∈ Rp×2 . where
G = ΓT Γ γ(0) = γ(0) γ(1) γ(1) T T γ(0) γ(0) γ(0) γ(1) = γ(1)T γ(0) γ(1)T γ(1) k(0, 0) k(0, 1) = . k(1, 0) k(1, 1)
γ(y).
T
Diffeomorphic Time Warping (DiffTW)
37
Then we can write "
k0 (x, y) = γ(x)T I − γ(0), γ(1) G−1 γ(0), γ(1)
T # γ(y)
= γ(x)T [I − ΓG−1 ΓT ]γ(y). Let P = ΓG−1 ΓT . Then P satisfies the properties of a projection matrix onto the column space of Γ. This follows from the standard projection formula P = Γ(ΓT Γ)−1 ΓT when projecting on the column space of a matrix Γ. P is symmetric and idempotent: P
= ΓG−1 ΓT = (ΓG−1 ΓT )T = PT
P 2 = (ΓG−1 ΓT )2 = ΓG−1 ΓT ΓG−1 ΓT = ΓG−1 GG−1 ΓT = ΓG−1 ΓT = P. Hence P is an orthogonal projection matrix. Thus, the matrix
A = I − P = I − ΓG−1 ΓT ,
represents the projection onto the orthogonal complement of the space spanned by γ(0) and γ(1). Since P projects onto a 2-dimensional subspace, rank(P ) = 2. Then A = I − P is the complement of P , hence by the Rank-Nullity Theorem
rank(A) + nullity(A) = p,
and since nullity(A) = rank(P ) = 2, we get rank(A) = p − 2. Now we can represent the modified kernel as: k0 (x, y) = γ(x)T Aγ(y). By spectral decomposition, A = U ΛU T where U T U = I, U = (u1 , . . . , up−2 ), the matrix
38
Vicky Geneva Haney, Kamel Lahouel, Victor Rielly, and Bruno M. Jedynak
of eigenvectors of A and Λ the matrix of positive p − 2 eigenvalues of A. Then,
γ0 (x) = Λ1/2 U T γ(x),
is the feature vector of the explicit kernel k0 with dimension p − 2.