[MRG] Fix GAK returning NaN for time series longer than ~405 samples (#450)

[samim-reza]

Summary

gak and cdist_gak returned NaN for every pair of time series longer than
about 405 samples. This normalizes the Global Alignment Kernel in log space, so
the result stays accurate whatever the length of the inputs.

Fixes #450 (and the duplicates #510, #188 that were closed in its favour).

Problem

The GAK recursion sums over every alignment path, so the unnormalized kernel
grows like the central Delannoy numbers, roughly $(3 + 2\sqrt{2})^{sz}$. It
therefore leaves the range of a 64-bit float at sz ≈ 406:

import numpy as np
from tslearn.metrics import gak, cdist_gak

x = np.sin(2 * np.pi * np.linspace(0, 1, 500)).reshape(-1, 1)
gak(x, x, sigma=5.0)          # nan   -- and gak(x, x) is 1 by construction
cdist_gak(np.stack([x, x]))   # all nan

gak normalized by dividing those already-overflowed values, so
inf / inf turned every result into NaN. Because TimeSeriesSVC,
TimeSeriesSVR and KernelKMeans build their kernel matrix with
cdist_gak, they failed outright on long series:

from tslearn.generators import random_walk_blobs
from tslearn.svm import TimeSeriesSVC

X, y = random_walk_blobs(n_ts_per_blob=5, sz=500, d=1, n_blobs=2, random_state=0)
TimeSeriesSVC(kernel="gak").fit(X, y)
# ValueError: Input X contains NaN.

That is the symptom originally reported in #188 and #450, and #514 raised the
threshold from ~204 to ~405 without removing it.

Solution

The quantity users ask for is the normalized kernel, which is always in
(0, 1] and perfectly well defined at any length — only the intermediate
values overflow. So the normalization is now done in log space:

gak(x, y) = exp( log k(x, y) - ½·log k(x, x) - ½·log k(y, y) )

and cdist_gak normalizes its matrix the same way instead of multiplying by
diagonal matrices of 1/sqrt(k(x, x)).

To obtain log k(x, y), _log_unnormalized_gak runs the existing linear
recursion first and only re-runs a new log-space recursion (log-sum-exp over
the three predecessors) when the linear value has overflowed to inf or
underflowed to 0. Series that already worked therefore keep both their
current values and their current speed; the extra cost is paid only in the
range that used to return NaN.

The log-space recursion is added through the same __make_..._from_gram_matrix
factory as the existing one, so the numba and PyTorch backends stay in sync.

unnormalized_gak is unchanged and still returns inf past ~405 samples —
the exact value really is outside float64 range — which its docstring and the
user guide now state explicitly.

Testing

• tests/test_metrics.py::test_gak_long_time_series — new non-regression test:
  gak(x, x) == 1, cdist_gak finite with unit diagonal, and the
  one/two-dataset paths agree, at sz=500.
• tests/test_svm.py::test_gak_svm_long_time_series — new non-regression test
  covering the originally reported TimeSeriesSVC / TimeSeriesSVR failure.
• tests/test_metrics.py::test_gak_log_and_linear_recursions_agree — new test
  that the two recursions give the same answer where both are valid.

Both non-regression tests fail on main (NaN, and ValueError: Input X contains NaN) and pass here.

Full suite on Linux / Python 3.12 with all_features, before and after:

main:  1063 passed, 31 skipped, 13 xfailed, 34 xpassed
here:  1070 passed, 31 skipped, 13 xfailed, 34 xpassed

(the difference is the 3 new tests plus 4 doctests of tslearn/metrics/_gak.py;
no other test changed status). flake8 reports no new warnings on the touched
files. The PyTorch backend was exercised for gak, unnormalized_gak and both
cdist_gak paths, and its log-space recursion was checked against the numba one.

Timings for gak on two random walks (numpy backend, best of 5):

sz	main	this PR
300	6.0 ms	5.2 ms
600	34.7 ms	56.3 ms

sz=300 uses the linear recursion, as before. sz=600 pays for the log-space
pass, and previously returned NaN.

[Copilot]

Pull request overview

Fixes GAK returning NaN for long time series by adding log-space normalization and recursion.

Changes:

• Added stable log-space GAK computation.
• Updated normalization, documentation, and changelog.
• Added regression tests for metrics and SVM usage.

Reviewed changes

Copilot reviewed 5 out of 5 changed files in this pull request and generated no comments.

Show a summary per file

File	Description
tslearn/metrics/_gak.py	Implements stable log-space GAK computation.
tests/test_svm.py	Tests SVMs with long GAK inputs.
tests/test_metrics.py	Tests long-series and recursion behavior.
docs/user_guide/kernel.rst	Documents overflow behavior and normalization.
CHANGELOG.md	Records the GAK fix.

Suppressed comments (2)

tests/test_metrics.py:555

• The new fallback implementation selected when PyTorch is available is not exercised by the added regression coverage: this test imports and calls only the _njit_* functions, while the 500-sample GAK and SVM tests use NumPy arrays. Please add a long-series be=Backend("pytorch") case covering both the direct and cdist_gak paths so the newly added PyTorch recursion and normalization are verified in CI.

    from tslearn.metrics._gak import (
        _njit_gak_from_gram_matrix,
        _njit_log_gak_from_gram_matrix,
    )

tslearn/metrics/_gak.py:291

• The fallback only inspects the final linear result. An individual entry of exp(log_gram) can underflow to zero while other alignment paths keep value finite, so this branch accepts a result that has silently discarded all paths through that entry; with long series, the many discarded paths can contribute materially. Treat any zero gram entry as an unusable linear pass and run the log recursion in that case.

    value = gak_from_gram(backend.exp(log_gram))
    if 0.0 < value < backend.inf:
        return backend.log(value)

---

💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 95.20%. Comparing base (4658696) to head (ae5aa89).

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #715      +/-   ##
==========================================
+ Coverage   95.15%   95.20%   +0.04%     
==========================================
  Files          83       83              
  Lines        7822     7896      +74     
==========================================
+ Hits         7443     7517      +74     
  Misses        379      379              

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[rtavenar]

Hi!

Thanks for this PR, it is of great interest for our lib! I just have one naive question: why not directly perform the computation in log space, to save time? In other words, why do you need the call to gak_from_gram in the first place?

[samim-reza]

Hi @rtavenar, thanks for the quick look!

Not naive at all — it's a deliberate trade-off, but the numbers are the real answer.

The linear pass is there because the log-space recursion is much more expensive per cell: the inner loop trades two additions and a multiply for three exp and a log. Timing the numba kernels alone (square gram, best of 7):

sz	linear	log-space
200	0.08 ms	1.08 ms	14x
400	0.35 ms	4.52 ms	13x
800	1.57 ms	18.4 ms	12x

So it isn't a free swap: going straight to log space would slow down every user whose series are short enough to work today. End-to-end on gak (numpy backend, best of 5):

sz	this PR	always log-space
100	0.76 ms	1.49 ms
300	4.67 ms	11.4 ms
400	12.5 ms	24.4 ms
600	51.9 ms	52.7 ms
800	116 ms	95 ms

Below the threshold the linear pass is ~2x faster; above it you're right that the speculative pass is wasted, and it costs ~20%.

A way to get both, if you want it

The overflow point depends only on the lengths, so it can be decided up front instead of by trial. Every gram entry is u / (2 - u) with u = exp(-d / 2σ²) ∈ (0, 1], hence ≤ 1, so the unnormalized value is bounded by the number of alignment paths — the Delannoy number D(sz1-1, sz2-1).

That bound turns out to be tight exactly where users hit the bug:

log(DBL_MAX)      = 709.78
log D(404, 404)   = 708.59   -> sz = 405 cannot overflow
log D(405, 405)   = 710.35   -> sz = 406 can

which is precisely the 405/406 cutoff reported in #450 / #510 / #188. So _log_unnormalized_gak could pick the recursion from sz1, sz2 in O(1) and never run a speculative pass — keeping the linear speed below the threshold and dropping the ~20% overhead above it.

Happy to push that, or to drop the linear path entirely if you'd rather keep the code simple and accept the ~2x on short series. Your call — just say which and I'll update.

Two smaller notes while I'm here:

• I checked the concern that a gram entry underflowing to 0 could make the linear pass silently drop paths: over 400 randomised cases where 4–13% of gram entries underflowed and the guard accepted the linear result, it agreed with the log-space value to 4e-16. A path through a cell with gram < 5e-324 carries weight ≤ 5e-324, so it only matters when essentially every path is killed — and then the total underflows to 0, which the 0.0 < value guard already rejects.
• The PyTorch log-space recursion is currently only verified locally, not in CI (a 500-sample pure-Python DP would be far too slow there). I can add a cheap test that exercises it directly on a small gram matrix against the numba one — happy to include it in whichever change you prefer.

[charavelg]

The current implementation in this PR seems legit: no overcost for existing working code is good IMHO.