Skip to content

[ENH] Implementation for Dynamic Time Warping (DTW) in cython - #5

Open
fnhirwa wants to merge 6 commits into
sktime:mainfrom
fnhirwa:dtw
Open

fnhirwa wants to merge 6 commits into
sktime:mainfrom
fnhirwa:dtw

Conversation

@fnhirwa

@fnhirwa fnhirwa commented Jul 27, 2026 •

Copy link
Copy Markdown
Member

Reference Issues/PRs

Part of #4

What does this implement/fix? Explain your changes.

This PR implements ahead-of-time (AOT) compiled Cython dynamic time warping (DTW) kernels (dtw_distance and dtw_cost_matrix), providing a zero-warmup, high-performance C-extension alternative to sktime's Numba implementation.

Key Changes

  1. Cython DTW Kernels (dtw_distance & dtw_cost_matrix)

    • Ported the dependent multivariate squared-Euclidean local cost (DTW_D) recurrence to Cython memoryviews (double[:, ::1]) with nogil.
    • Memory Optimization: dtw_distance computes the final DTW cost via a rolling two-row buffer, reducing memory overhead to $O(m_2)$ instead of allocating a full $O(m_1 \times m_2)$ matrix when alignment path recovery isn't required.
  2. Precomputed Vectorized Masking (unsigned char[:, ::1])

    • Replaced in-kernel floating-point sentinel checks with a precomputed 2D uint8 mask (np.isfinite(bounding_matrix)).
    • Sidesteps cross-platform compiler optimizations and strict-aliasing issues under -O2/-O3 / -ffast-math (e.g., Clang on macOS), providing ultra-fast $O(1)$ cell validity lookups inside nogil loops.

Benchmark Results

Below is the end-to-end performance comparison between sktime's Numba implementation (sktime.dists_kernels._numba_distances.dtw_distance) and the new ahead-of-time compiled Cython extension.

Times in microseconds (µs). Benchmark measured across 500 calls with pre-allocated bounding matrices.

Series Shape (d, m) Numba dtw_distance (µs) Cython dtw_distance (µs) Speedup Cython cost_matrix (µs)
d=1, m=50 24,681.96 10.74 2,298.51x 9.84
d=1, m=200 30,798.19 148.33 207.63x 141.62
d=3, m=200 31,617.42 151.80 208.29x 159.52
d=5, m=500 125,121.70 930.15 134.52x 970.37

Key Performance

  • Sub-Millisecond Execution: Even for multi-channel series (d=5, m=500), distance calculation finishes in < 1 ms (930 µs).
  • Zero Dispatching Overhead: Bypasses pure-Python wrapper allocation and argument verification overhead (~25–30 ms saved per call).
  • No JIT Warmup: AOT compilation avoids cold-start latency on first invocation.

@sssilvar sssilvar left a comment •

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I found several issues:

  1. [P1] Validate channel counts before entering the unchecked Cython loop

In _dtw_cython.pyx, _local_cost loops over x.shape[0] and reads y[k, j] with bounds checking disabled. If x has more channels than y, this reads beyond the memory allocated for y. A reproduction with shapes (2, 4) and (1, 4) returned the incorrect value 4.0 instead of raising an error. Both public entry points should require _x.shape[0] == _y.shape[0].

  1. [P1] Validate the bounding matrix shape

A custom bounding matrix is accepted without ensuring its shape is (m1, m2). The Cython kernels then access bm_mask[i, j] with bounds checking disabled, so an undersized matrix causes out-of-bounds reads and undefined behavior. This also affects automatically generated Itakura bounds for unequal-length series: _itakura_parallelogram creates (y_size, x_size), while the kernels require (x_size, y_size).

  1. [P1] Do not combine the infinity-based recurrence with incompatible fast-math assumptions

The extension uses INFINITY to initialize and block dynamic-programming cells, but setup.py compiles it with -ffast-math. Clang reports multiple use of infinity via a macro is undefined behavior warnings while building this PR. Precomputing the bounding mask avoids finite checks in the loop, but it does not protect the infinity values used by the recurrence. The DTW extension should be compiled without incompatible finite-math assumptions, or use a finite sentinel with overflow-safe handling.

  1. [P2] Expand legacy equivalence coverage beyond equal-length scalar distances

The current test_cython_matches_numba is a useful start, but it only exercises equal-length inputs and dtw_distance. Please mirror the Rocket regression approach across unequal lengths in both orientations, univariate and multivariate inputs, all bounding modes, custom masks, and dtw_cost_matrix. Add explicit validation tests for mismatched channel counts and incorrectly shaped bounding matrices. At least one regular CI job should install the dev extra so these comparisons run instead of being skipped; wheel tests can continue to skip them.

@fnhirwa

fnhirwa commented Sep 11, 2026

Copy link
Copy Markdown
Member Author

Thanks for the review!
all four items are addressed: channel-count and bounding-matrix shape validation on both entry points and in the kernels, the DTW extension now compiles without -ffast-math, and the legacy equivalence tests are expanded across unequal lengths in both orientations, univariate/multivariate, all bounding modes, custom masks, and dtw_cost_matrix, with one CI job set to fail rather than skip them.
Ready for another look.

@fnhirwa
fnhirwa requested a review from sssilvar September 11, 2026 00:16
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants