[ENH] Add MADRID multi-length discord anomaly detector

[JayeshSuryavanshi]

Reference Issues/PRs

Fixes #1710.

What does this implement/fix? Explain your changes.

This PR adds MADRID, a multi-length discord anomaly detector for series data, based on Lu, Srinivas, Nakamura, Imamura and Keogh, "Matrix Profile XXX: MADRID: A Hyper-Anytime and Parameter-Free Algorithm to Find Time Series Anomalies of All Lengths" (ICDM 2023). MADRID is the faster successor to the existing MERLIN detector: instead of committing to a single subsequence length, it runs the DAMP left-discord matrix-profile method across a candidate set of subsequence lengths and combines the length-normalised discord profiles into a single per-point anomaly score.

• New estimator aeon/anomaly_detection/series/distance_based/_madrid.py with MADRID(BaseSeriesAnomalyDetector), mirroring the MERLIN file structure and using numba @njit for the DAMP/MASS inner loops.
• _predict returns a real-valued per-point anomaly score (the column-wise maximum of the length-normalised multi-length discord table). Tags: anomaly_output_type="anomaly_scores", learning_type:unsupervised=True.
• Parameters: min_length, max_length, step_size, and train_test_split (accepts an int index, a float fraction, or None for a default warm-up). Input validation mirrors MERLIN, plus a constant-region warning.
• Registered in aeon/anomaly_detection/series/distance_based/__init__.py and in the API reference; tests added in tests/test_madrid.py.

Verification: passes check_estimator(MADRID) (20/20 checks), the docstring doctest, and ruff/pre-commit. The output was also validated against an unofficial reference implementation of the paper: same detected discord location, and the full multi-length discord table matches to within floating-point tolerance (~1e-13).

Does your contribution introduce a new dependency? If yes, which one?

No. It relies only on numba, which is already a core aeon dependency (as MERLIN is). The z-normalised MASS distance profile is computed with a direct sliding dot product, deliberately avoiding any FFT / rocket-fft dependency.

Any other comments?

Disclosure: this contribution was developed with the help of an AI coding assistant. I have reviewed and understand the implementation, verified it against the paper and a reference implementation, and will maintain it.

PR checklist

For all contributions

• I've added myself to the list of contributors (will use the @all-contributors bot after merge).
• The PR title starts with [ENH].

For new estimators and functions

• I've added the estimator to the online API documentation.
• (OPTIONAL) I've added myself as a __maintainer__ at the top of relevant files.

[aeon-actions-bot]

Thank you for contributing to aeon

I have added the following labels to this PR based on the title: [ enhancement ].
I have added the following labels to this PR based on the changes made: [ anomaly detection ]. Feel free to change these if they do not properly represent the PR.

The Checks tab will show the status of our automated tests. You can click on individual test runs in the tab or "Details" in the panel below to see more information if there is a failure.

If our pre-commit code quality check fails, please run pre-commit locally and push the fixes to your PR branch.

Don't hesitate to ask questions on the aeon Discord channel if you have any.

PR CI actions

These checkboxes will add labels to enable or disable CI functionality for this PR. This may not take effect immediately, and a new commit may be required to run the new configuration.

• Run pre-commit checks for all files
• Run mypy typecheck tests
• Run all pytest tests and configurations
• Run all notebook example tests
• Run numba-disabled codecov tests
• Disable numba cache loading
• Regenerate expected results for testing
• Push an empty commit to re-run CI checks

[JayeshSuryavanshi]

For context: this picks up from the earlier draft #1846 by @acquayefrank, which was closed as outdated with an invitation to continue. The implementation here is written from scratch to match aeon's current estimator format (BaseSeriesAnomalyDetector, tags, _get_test_params, tests, and a doctest). Happy to coordinate if @acquayefrank would still like to be involved.

[TonyBagnall]

hi, thanks for this, is there an external package we can test against for correctness and speed? Dont need exact numeric similarity, but some reassurance that its not sig diff to some reference implementation would be good

[JayeshSuryavanshi]

Hi Tony, thanks for taking a look.

Yes, there is a reference implementation to check against: k-kotera/MADRID-python, an unofficial numba port of the authors' MATLAB code from the ICDM 2023 paper. It is what I validated this implementation against while writing it.

Correctness. I compared the full multi-length discord tables (not just the final score profile) on synthetic series with injected anomalies. On an 8192-point series with lengths 16..64 step 4 (13 lengths):

• discord table max abs difference: 9.0e-14
• best-so-far discord scores max abs difference: 4.4e-15
• best-so-far discord locations: identical for every length
• top discord location: identical (5028, inside the anomaly injected at 5000:5060)

The residual 1e-14 noise is consistent with fastmath reassociation in the numba kernels here; the two implementations are otherwise numerically identical. Other series and length grids I tried during development gave the same picture (max table diff 5.8e-14, same discord locations).

Speed. On the same 8192-point series, after JIT warm-up for both, this implementation ran in ~0.11 s vs ~21 s for the reference port on my machine (Apple Silicon). I have not dug into why the reference port is slower, but this one is not giving anything away on speed.

MADRID-python is not packaged (loose madrid.py/damp.py files, not on PyPI), so it is awkward as a CI dependency. If you would like a regression guard, I am happy to add a test that pins the discord table values from the reference on a small fixed series.

Reproduction script (needs a clone of MADRID-python on sys.path)

import sys
import time

import numpy as np

sys.path.insert(0, "/path/to/MADRID-python")  # clone of k-kotera/MADRID-python
from madrid import MADRID as madrid_ref

from aeon.anomaly_detection.series.distance_based._madrid import _madrid

rng = np.random.default_rng(42)
n, split = 8192, 2048
x = np.sin(np.linspace(0, n / 20 * np.pi, n)) + rng.normal(0, 0.05, n)
x[5000:5060] += 2.0
m_set = np.arange(16, 65, 4, dtype=np.int64)

# JIT warm-up for both implementations
warm = x[:1024].copy()
madrid_ref(warm, 256, 16, 32, 4, "manual")
_madrid(warm, 256, np.arange(16, 33, 4, dtype=np.int64))

t0 = time.perf_counter()
_, table_ref, bsf_ref, bsf_loc_ref, m_set_ref = madrid_ref(x, split, 16, 64, 4, "manual")
t_ref = time.perf_counter() - t0

t0 = time.perf_counter()
table, bsf, bsf_loc = _madrid(x, split, m_set)
t_aeon = time.perf_counter() - t0

assert np.array_equal(m_set, m_set_ref)
print("discord table max |diff|:", np.max(np.abs(table - table_ref)))
print("best-so-far scores max |diff|:", np.max(np.abs(bsf - bsf_ref)))
print("best-so-far locations equal:", np.array_equal(bsf_loc, bsf_loc_ref))
print("top discord location: ref =", bsf_loc_ref[np.argmax(bsf_ref)],
      "| aeon =", bsf_loc[np.argmax(bsf)])
print(f"runtime: reference {t_ref:.2f}s | aeon {t_aeon:.2f}s")

[TonyBagnall]

I think the core MADRID implementation is fine, but I m not sure how the result is converted into an aeon anomaly score.

BaseSeriesAnomalyDetector defines float output as a pointwise anomaly score, i.e. each value should represent the anomalousness of the corresponding time point. MADRID/DAMP naturally produces scores for subsequences starting at each position and for multiple candidate lengths. The current

return np.max(discord_table, axis=0)

therefore gives the maximum score of a subsequence starting at each point, rather than a pointwise anomaly score.

For comparison, other window-based aeon anomaly detectors such as STOMP map window scores back onto the points covered by those windows using reverse_windowing.

I think MADRID should do the same adaptation. For each candidate length, map its subsequence scores back onto the original time points, then combine the resulting pointwise profiles across lengths, probably using the maximum across lengths. The exact aggregation of overlapping windows (mean or max) needs a little thought, but I don't think taking max directly over the raw discord table has the right aeon semantics.

The tests should then also check the pointwise result, e.g. that an injected anomalous interval receives elevated scores across the interval, rather than only checking that argmax is close to its starting location.

@SebastianSchmidl any thoughts?

[JayeshSuryavanshi]

Thanks @TonyBagnall, that's a fair point about the semantics. Fixed in 55d618c: within each length the subsequence scores are mapped onto the points they cover via reverse_windowing, then combined across lengths with max, so a point scores as the most anomalous subsequence of any length covering it.

I used max rather than STOMP's mean within a length: a discord score is a property of the whole subsequence, and a mean would dilute a single-window discord with its normal neighbours. Happy to switch to mean for consistency if you'd prefer.

One caveat on the tests. I assert the detected discord's cover is elevated as a block and beats the background median, but not calibrated scores across the whole injected interval, because DAMP's pruning leaves unrefined estimates away from the top discord. MADRID only guarantees the top-discord region, which is the trade that makes it fast. If a full calibrated profile is the requirement, that points to STOMP rather than a change here.

Curious what @SebastianSchmidl thinks on mean vs max too.

[TonyBagnall]

thanks for this, it all looks structurally good but we are a bit pedantic about submitting new estimators. The main issue is how the scores are reconstructed to pointwise anomalies. The new reverse_windowing(..., reduction=np.max) fixes the previous issue that the raw output was indexed by subsequence start rather than time point, but DAMP/MADRID does not calculate an exact discord score at every subsequence position: pruned positions can inherit the previous value.

That means projecting the complete discord_table back onto the series may assign pointwise anomaly scores based partly on values that were never actually evaluated for those windows. The tests also acknowledge that only the top-discord region has a strong guarantee.

Could we either:

1. justify why the full approximate/pruned table is appropriate as a pointwise anomaly profile, or
2. construct the aeon output from the discord locations/scores that MADRID actually identifies?

also I think it would be good to be able to recover the more complete information returned by the reference implementation? Cant store these from predict, but could have something like

def _predict(self, X):
    M, bsf, locations, lengths = self._run_madrid(X)
    return self._to_pointwise_scores(M, lengths)

def predict_discords(self, X):
    M, bsf, locations, lengths = self._run_madrid(X)

[TonyBagnall]

do these parameters come from the reference implementation? I think they match MERLIN in aeon, but the original its data adaptive I think? We could add a None option? Not a blocker, just a thought

[TonyBagnall]

MERLIN uses 5 and 50, I'm not sure whats best

[TonyBagnall]

Thanks, these changes address my main concern. Building the aeon pointwise output only from the per-length discords actually identified by MADRID is much cleaner than projecting the complete pruned table, and exposing the underlying results separately with predict_discords is useful. The fit_is_empty and _mass fastmath changes also look good, especially given that you rechecked the numerical agreement after enabling fastmath.

A couple of remaining points:

1. Could predict_discords also return the overall best discord interval? The reference implementation returns best_discord_loc in addition to M, BSF_m, BSF_loc_m and m_set, and the new method currently exposes equivalents of the latter four only. Since it is just the location/length corresponding to argmax(scores), it should be essentially free to retain this part of the native MADRID output too.

2. the reference automatic mode derives maxL from the split and samples about 50 candidate lengths, whereas this implementation defaults to every integer length from 8 to 50. Is there a reason for deliberately dropping the automatic/reference length selection? An explicit aeon interface is fine, but perhaps max_length=None could retain MADRID's automatic behaviour, particularly given the algorithm's parameter-free motivation?

[JayeshSuryavanshi]

Both done in 7c3125d. predict_discords now returns best_interval, the (start, end) of the argmax discord, matching the reference's best_discord_loc. And max_length=None is now the default, reproducing the reference's automatic mode: max length from split // 20 and ~50 unique candidate lengths sampled via linspace. An explicit max_length keeps the previous every-integer sweep with step_size. No reason for dropping it originally beyond having mirrored MERLIN's interface; the parameter-free default is truer to the paper.