# Leveraging Data Seasonality and Matrix Profile for Anomaly Detection: Application to Climate Time Series

Guillaume Coulaud, Reza Akbarinia, Audrey Brouillet, Florent Masseglia

► To cite this version:

Guillaume Coulaud, Reza Akbarinia, Audrey Brouillet, Florent Masseglia. Leveraging Data Seasonality and Matrix Profile for Anomaly Detection: Application to Climate Time Series. RR-9572, Inria. 2025. hal-04906596

HAL Id: hal-04906596

https://inria.hal.science/hal-04906596v1

Submitted on 23 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Leveraging Data Seasonality and Matrix Profile for Anomaly Detection: Application to Climate Time Series

Guillaume Coulaud, Reza Akbarinia, Audrey Brouillet, Florent Masseglia

## RESEARCH REPORT

Nº 9572

January 2025

Project-Teams Iroko

ISSN 0249-6399

ISRN INRIA/RR-9572-FR+ENG
<figure><img src="image_3.png" /></figure>
Inria

# Leveraging Data Seasonality and Matrix Profile for Anomaly Detection: Application to Climate Time Series

Guillaume Coulaud*, Reza Akbarinia†, Audrey Brouillet‡,
Florent Masseglia†

Project-Teams Iroko

Research Report n° 9572 — January 2025 — 29 pages

**Abstract:** Seasonal time series analysis is fundamental in domains such as climate science, where detecting and understanding anomalies, patterns, and data changes are essential. The classical Matrix Profile approach does not consider the data's seasonality, failing to detect seasonal anomalies and patterns. This paper introduces the Interval Matrix Profile, a novel extension of the Matrix Profile specifically designed for analyzing periodic and seasonal time series data. The Interval Matrix Profile enables flexible interval-based comparisons across seasons, allowing the detection of anomalies that conventional approaches miss. We further propose the constrained k Nearest Neighbor Interval Matrix Profile, designed to identify anomalies that may appear across multiple periods, a common characteristic of abnormal climate events and extreme weather phenomena. Our approach leverages a scalable block-based algorithm that achieves significant performance gains through caching, vectorization, and parallelism. Additionally, we introduce a novel methodology to detect the first or last occurrence of a pattern, enabling the discovery of pattern emergence or disappearance within seasonal time series. The algorithms are demonstrated in case studies on temperature climate time series. They effectively capture seasonal anomalies and find pattern disappearance. Our results illustrate that the IMP consistently outperforms the classical Matrix Profile in the accuracy of seasonal anomaly detection and computational efficiency.

**Key-words:** anomaly detection, matrix profile, time series, climate

<footnote>* Université Montpellier, Inria, CNRS, LIRMM</footnote>

<footnote>† Inria, Université Montpellier, CNRS, LIRMM</footnote>

<footnote>‡ ESPACE-DEV, Université Montpellier, IRD, Université Guyane, Université Réunion, Université Antilles</footnote>
# Exploitation de la saisonnalité des données et du Matrix Profile pour la détection des anomalies : application aux séries temporelles climatiques

**Résumé** : L'analyse des séries temporelles saisonnières est fondamentale dans des domaines tels que la science du climat, où la détection et la compréhension des anomalies, des patterns et des évolutions sont essentielles. L'approche classique du Matrix Profile ne tient pas compte de la saisonnalité des données et ne permet pas de détecter les anomalies et les patterns saisonniers. Cet article présente l'Interval Matrix Profile, une nouvelle extension du Matrix Profile spécialement conçue pour l'analyse des données de séries temporelles périodiques et saisonnières. L'Interval Matrix Profile permet des comparaisons flexibles basées sur des intervalles répartis entre les saisons, ce qui permet de détecter des anomalies que les approches conventionnelles ne parviennent pas à déceler. Nous proposons, en outre, le constraint k Nearest Neigbors Interval Matrix Profile, conçu pour identifier les anomalies qui peuvent apparaître sur plusieurs périodes, une caractéristique commune des événements climatiques anormaux et des phénomènes météorologiques extrêmes. Notre approche s'appuie sur un algorithme efficace basé sur des blocs permettant d'obtenir des gains de performance significatifs grâce à la mise en cache, à la vectorisation et au parallélisme. Enfin, nous introduisons une nouvelle méthodologie pour détecter la première ou la dernière occurrence d'un modèle, ce qui permet de découvrir l'émergence ou la disparition d'un modèle dans les séries chronologiques saisonnières. Les algorithmes sont démontrés dans des études de cas portant sur des séries chronologiques de température et de climat. Ils capturent efficacement les anomalies saisonnières et détectent la disparition des modèles.

**Mots-clés** : détection d'anomalies, matrix profile, séries temporelles, climat
# Contents

<table><tr><td>1</td><td>Introduction</td><td>4</td></tr><tr><td>2</td><td>Background</td><td>6</td></tr><tr><td>3</td><td>Related Work</td><td>7</td></tr><tr><td>4</td><td>Seasonal and Interval Matrix Profile</td><td>8</td></tr><tr><td>4.1</td><td>Seasonal Matrix Profile</td><td>8</td></tr><tr><td>4.2</td><td>Interval Matrix Profile</td><td>9</td></tr><tr><td>4.3</td><td>Constrained kNN Interval Matrix Profile</td><td>15</td></tr><tr><td>4.4</td><td>Finding the Disappearance of a Pattern using the Interval Matrix Profile</td><td>16</td></tr><tr><td>5</td><td>Experimental Evaluation</td><td>17</td></tr><tr><td>5.1</td><td>Setup</td><td>18</td></tr><tr><td>5.2</td><td>Performance Evaluation</td><td>19</td></tr><tr><td>5.3</td><td>Precision Evaluation</td><td>20</td></tr><tr><td>5.4</td><td>Trivial Neighbors in the Classical kNN MP</td><td>21</td></tr><tr><td>5.5</td><td>Climate Case Study</td><td>21</td></tr><tr><td>6</td><td>Conclusion</td><td>24</td></tr></table>
# 1 Introduction

The detection of anomalies in time series has gained significant attention [1] across diverse fields, driven by its importance in identifying unusual patterns that may signal critical events or trends. From financial markets [2] to healthcare [3] and industrial monitoring [4], anomaly detection methods provide valuable insights for decision-making and predictive analysis. In the context of climate time series, anomaly detection plays a crucial role in uncovering extreme weather events, monitoring long-term environmental changes, and informing climate policies and adaptation strategies [5, 6]. Various techniques [7, 8, 9] have been developed to identify anomalies in time series. These techniques may be based on clustering, distances, densities, or statistics. This study focuses on the Matrix Profile (MP) [10], a well-established method for finding patterns and anomalies.

The MP is a distance-based approach for identifying the most similar or dissimilar subsequences within a time series. For a time series, the MP computes the closest matching distance for each subsequence in the time series, resulting in a series of values that represent the nearest-neighbor distance for each subsequence. Although it is a powerful tool, the MP faces significant challenges, both in terms of computational efficiency and comparison strategy. If this strategy is flawed, the resulting data will be of poor quality. This is especially true for climate-related data, which are inherently periodic (e.g., annual cycles), seasonal (e.g., summer or winter seasons), and cyclical (e.g., phenomena like El Niño). Furthermore, climate data exhibit significant variability due to the complexity of the underlying physical systems driving the different phenomena. These characteristics make analyzing and detecting anomalies in such data particularly challenging, as effective methodologies must account for interactions between periodicity, seasonality, cyclicity, and variability. For instance, comparing a winter subsequence to a summer subsequence solely based on similarity could lead to misinterpretation. Indeed, an abnormally warm winter week could be incorrectly classified as a normal subsequence when compared to a similar summer subsequence. Such a comparison not only fails methodologically but also risks obscuring critical anomalies, such as abnormal seasonal temperature extremes.

Consider the real-world example in our data as illustrated in Figure 1. It shows a scenario where a sudden drop in temperature is not detected by the classic MP approach. In this example, we compute the MP on the daily averaged Sea Surface Temperature (SST) from the Mediterranean Sea using ERA5 data [11, 12]. We visualize the months of June and October across all available years, highlighting in red a seven-day subsequence that can be considered abnormal from a climatological point of view, and in blue its closest match. The sudden one-degree Kelvin drop in temperature during June (in red) is not considered abnormal, since a similar subsequence (close match) is found in October (in blue). Interestingly, as the figure shows, such a temperature drop is typical for October but highly abnormal for June. This underscores the importance of accounting for seasonality when applying distance-based methods such as the MP for climate data anomaly detection.

Furthermore, the MP only focuses on identifying the closest matches. However, in climate, anomalies are expected to recur, so the closest match of an anomaly may also be an anomaly, and both of them would not be detected. To address these two limitations and develop a more robust tool for climatological analysis, we introduce an algorithmic enhancement of the MP: the Interval Matrix Profile, which incorporates seasonality into the analysis. This enhancement is combined with a reformulation of the kNN MP problem [13] to manage recurring anomalies effectively.

Our contributions are summarized as follows.

* We formulate the problem of the Interval Matrix Profile (IMP) for seasonal time series and propose the Block IMP (BIMP), an efficient algorithm based on blocking techniques
<figure><img src="image_4.png" /><figcaption>Figure 1: A daily averaged Sea Surface Temperature (SST) anomaly analysis from the Mediterranean Sea in June and October over the years. The seven-day subsequence (in red) and its closest match (in blue) is a temperature drop of 1°K, which is typical for October, but abnormal in June.</figcaption></figure>

to compute it in parallel. We also extend it to the left and right Interval Matrix Profile.

* We introduce the constrained kNN Interval Matrix Profile (ckNN IMP) to handle recurring discords.

* We propose a novel methodology, to discover emerging or disappearing patterns, that combines the left and right IMP and the ckNN IMP.

* We conduct experiments on synthetic data highlighting the BIMP algorithm scalability and sensibility to its parameters.

* We evaluate our methods on climate time series. The results illustrate that i) the IMP allows discovering relevant discords that would be ignored by the original MP, ii) the constrained kNN IMP is highly relevant for the use case of climate anomalies, and iii) it can detect the last occurrence of a seasonal pattern.

The paper is organized as follows. In Section 2, we present the background of the Matrix Profile. We present the related work in Section 3. The core of our contribution is then described in Section 4 with the Seasonal Matrix Profile, the Interval Matrix Profile, and the constrained kNN Interval Matrix Profile. A performance evaluation of our methods is presented in Section 5. We conclude in Section 6.
## 2 Background

A time series is a sequence of data points collected over time, denoted as $T = \{T_1, T_2, T_3, \dots, T_n\}$ where $T_i$ represents the data point at time $i$ and $n$ is the length (or size) of the series. A **subsequence** is a smaller part with continuous values of this series, represented as $T_{i,m} = \{T_i, T_{i+1}, \dots, T_{i+m-1}\}$, where $m$ is the length of the subsequence and $i$ is the starting index of the subsequence.

To compare two subsequences, $T_{i,m}$ and $T_{j,m}$, we measure their distance, which is denoted as $d(T_{i,m}, T_{j,m})$, where $d$ is a metric such as the Euclidean distance. When computing the distance between subsequences, we have to consider **self or trivial-match** which is a subsequence significantly overlapping with the one considered. The length of the overlap is user-defined but is typically taken to be either $m/2$ or 0. In such cases, the distance between the starting indices of the two subsequences must exceed an exclusion zone, $e$, of length $m/2$ or $m$.

In the context of time series analysis, the analysis of similarity between successive subsequences allows the extraction of two key concepts: motifs and discords. A **motif** is defined as a subsequence that has a very close match within the time series. Conversely, a **-discord** is a subsequence in a time series that has the largest distance to its nearest neighbor, meaning that it stands out as being the most dissimilar compared to any other subsequence in the series. Formally a discord is defined as

$$ \max_{i} \min_{j} \{d(T_{i,m}, T_{j,m}) \mid |i-j| \ge e\}, $$

where $|i-j| \ge e$ describes the non-trivial match for an exclusion zone $e$. We denote the top $k$ discords as the $k$ subsequences that have the most distant nearest neighbors in the time series.

The problem in constructing the Matrix Profile is to consider the set of all subsequences of size $m$, and for each subsequence, identify the nearest subsequence that is not a trivial match. Considering a time series $T$ of size $n$, the time series is divided into a set of subsequences using a sliding window of size $m$ with a shift of 1. This leads to a set of $n-m+1$ subsequences. The **Matrix Profile** (MP) of $T$ is the vector of the distances between each subsequence and its nearest subsequences. Formally, the elements of the MP are defined as follows for $i = 1, \dots, n-m+1$,

$$ \text{MP}_i = \min_{j} \{d(T_{i,m}, T_{j,m}) \mid |i-j| \ge e\}. $$

An inefficient way for computing the MP would be to build the distance matrix $D$, with $D_{i,j}$ the distance between the subsequence $i$ and the subsequence $j$, then retain the smallest value per row. The computational complexity of such an algorithm is $O(dn^2)$ where $d$ is the complexity of the distance computation between two subsequences.

Usually, the MP is computed using the z-normalized Euclidean distance. This distance allows for finding subsequences whose shapes are abnormal without taking into consideration the amplitude of the observation values, making it effective for identifying shape-based anomalies. However, this approach may not be suitable for cases where amplitude differences are essential to identify anomalies. For example, in climate applications such as temperature time series, two subsequences may have similar shapes but vastly different temperatures, for example, a 10 °C difference. In such situations, using z-normalization may obscure these differences, classifying the subsequences as similar despite the important amplitude variation. For the climate case studies, the standard Euclidean distance is more appropriate as it accounts for both amplitude and shape. This distinction is especially valuable when analyzing climate data, where amplitude-based anomalies, such as sudden spikes or drops in temperature, are often of significant interest.

Naively the Euclidean distance for the MP computation can be explicitly computed using
Equation 1, leading to an inefficient algorithm where each distance is computed individually.

$$D_{i,j}^2 = \sum_{k=0}^{m-1} (T_{i+k} - T_{j+k})^2 \qquad (1)$$

However, the distance between two subsequences can be incrementally computed with a recurrence as in the AAMP [14] using the previous coefficient of the diagonal as described in Equation 2. However, in order to calculate a row of the distance matrix using the recurrence, it is necessary to access the coefficients of the previous row.

$$D_{i,j}^2 = D_{i-1,j-1}^2 - (T_{i-1} - T_{j-1})^2 + (T_{i+m-1} - T_{j+m-1})^2 \qquad (2)$$

To compute all coefficients of the matrix, the explicit method (Eq.1) has a complexity of $O(mn^2)$, while the incremental method (Eq.2), which exploits the aforementioned recurrence, has a complexity of $O(n^2)$. The first method will be denoted as the brute force algorithm, while the second is the incremental algorithm (or AAMP). Whenever possible, the distance should be computed incrementally. This is what the STOMP algorithm [15] does for the z-normalized Euclidean distance. However, in finite arithmetic, such incremental calculations introduce numerical errors [16, 17] that must be considered. The brute force algorithm is a generic algorithm that works with any distances or metrics.

## 3 Related Work

The MP is a proven method to enhance anomaly detection and pattern discovery in time series. It has been applied in various fields, such as seismology [18], climate studies [19], and traffic data analysis [20]. The MP technique was first introduced by Yeh et al. with the STAMP algorithm [10]. This algorithm employs an efficient distance calculation technique, known as MASS, which utilizes the Fast Fourier Transform to calculate the z-normalized Euclidean distance from a subsequence to all the others. Zhu et al. introduce the STOMP algorithm [15], which takes advantage of incremental computations to reduce the complexity of the distance calculation to $O(1)$. This algorithm is well suited for GPU devices and can also be distributed with the SCAMP algorithm [17, 21]. The previous algorithms use the z-normalized Euclidean distance, but other distances are possible, such as the classical Euclidean distance with AAMP or Minkowski distance [14], Dynamic Time Warping [22], or even subsequence features [23].

The MP can also be tailored to detect only anomalies with DAMP [24], and MADRID [25] its generalization to multiple subsequence lengths. These algorithms efficiently compute the left Matrix Profile, by only looking on the left of each subsequence, with early abandoning for subsequences that cannot be discords.

The MP is sensitive to the *twin freaks* problem (or recurring discords). A twin freak is a discord not discovered as the abnormal subsequence appears at least twice. Different works have been conducted to address this issue by introducing the *k*Nearest Neighbors Matrix Profile (kNN MP) [13, 26]. The kNN MP enhances anomaly detection by examining the *k*th nearest neighbor, thereby improving the identification of recurring discords. In the same spirit [27] Huynh et al. introduce the K-distance which is the sum of the distances from a subsequence to its k-nearest neighbors.

The MP performs a global request, where a subsequence is compared to all others. Some research has been conducted to modify this property. Zhu et al. proposed incorporating domain knowledge with an annotation vector [28], while Dau et al. introduced the STUMP algorithm [29] to facilitate arbitrary range queries.
Similarly to the two previous works, we propose to reduce the query range. The objective of this reduction is to take into account the seasonality of the data. Seasonality is a key aspect of anomaly detection in fields such as climatology. Our aim is therefore to provide a more robust and efficient variant of the MP for seasonal data.

# 4 Seasonal and Interval Matrix Profile

We first introduce the Seasonal Matrix Profile, specifically designed to deal with seasonality in time series. Next, we develop the Interval Matrix Profile, a robust enhancement of the initial method, and we present an efficient block algorithm for computing it. Then, we propose a constrained kNN version for the Interval Matrix Profile to detect recurrent discords. Lastly, we present a methodology to discover seasonal patterns that disappear using the Interval Matrix Profile and its kNN formulation.

## 4.1 Seasonal Matrix Profile

Consider a periodic time series in which each period is divided into distinct seasons. For example, a year represents a period divided into seasons as the common four climate seasons in northern and southern mid-latitudes or broader categories like the monsoon and dry seasons. To illustrate this, Figure 2 presents a periodic function over four periods. Each period is split into two seasons: one where the values increase (orange) and the other where the values decrease (blue). In this scenario, to identify an anomaly within a season, subsequences must only be compared with other subsequences from the same season. For example, subsequences from orange seasons should only be compared with subsequences from other orange seasons.

<figure><img src="image_5.png" /><figcaption>Figure 2: Simple example of a periodic time series with two seasons</figcaption></figure>

An intuitive approach to take into account the seasonality in time series data is to restrict the search for similar subsequences to those occurring within the same season across different periods. For example, a subsequence observed during winter is only compared to subsequences from other winters. We call such an approach the **Seasonal Matrix Profile (SMP)**. This seasonal restriction is interesting in climate analysis, where significant events, for example, a drought in the fall, could be missed if compared to subsequences from unrelated seasons like summer. Formally, given a time series T, the elements of the SMP vector are defined as follows for $i = 1, \dots, n - m + 1$,

$$ \text{SMP}_i = \min_{j} \{d(T_{i,m}, T_{j,m}) \mid |i - j| \ge e \wedge \text{Season}(i) = \text{Season}(j)\}, $$

where e is the exclusion length, and Season(·) is a function that maps an index to its corresponding season.

We visually highlight the SMP comparisons in Figure 3. Figure 3a depicts the subsequence perspective of the SMP. In this approach, each subsequence is compared only with other subsequences from the same season to find the closest match. For example, a subsequence (query) from Season 2, shown in yellow, is compared solely to subsequences from Season 2, shown in blue. Figure 3b presents the perspective of the SMP distance matrix that highlights the distances to compute. The matrix is structured into square blocks, where each block represents the pairwise
distances between subsequences within the same season for two years. For example, an orange block describes the distances of subsequences belonging to the first season. In contrast, white blocks represent distances that are not considered as they fall between different seasons.

<figure><img src="image_6.png" /><figcaption>Figure 3: Schematic representation of the Seasonal Matrix Profile with two seasons per period</figcaption></figure>

**Limitation.** By computing the Matrix Profile within each season, the SMP allows the detection of seasonal anomalies, trends, and patterns that may be missed when using the classical Matrix Profile. However, a notable limitation of this method is the abrupt seasonal partitioning, which does not account for the gradual transitions between seasons, as seen in meteorological data. For example, a subsequence starting on June 19 would not be compared to a subsequence starting on June 22 because the first falls in spring and the second in summer. This is a significant limitation for applications in climate and meteorology. Furthermore, an inherent challenge of the SMP is the selection of appropriate seasons. The beginning and end of the seasons are not necessarily constant; there can be a shift without this being anomalous. Conversely, the SMP, with its rigidity, will detect anomalies if the partitioning does not include the shifted dates, which will result in the generation of false discoveries.

To address these issues, we propose a more refined approach in the next section.

## 4.2 Interval Matrix Profile

The **Interval Matrix Profile** (IMP) enhances the capabilities of the SMP by allowing the comparison of subsequences within predefined intervals of the time series data. The core concept is similar to the SMP as it aims to restrict the search space, but the IMP defines this space in terms of fixed-size intervals centered on the subsequence's position within the period. This use of intervals allows for smoother transitions between seasons, avoiding the abrupt arbitrary boundaries that occur with strict seasonal partitioning.

To illustrate the IMP, consider an interval of length $L$ and a subsequence starting on June 21, 2023. The search space consists of all intervals, each of length $L$, centered on June 21 across all years. In other words, if $L = 100$ the subsequence is compared to all subsequences from May 2 to August 10 for all years.

Formally, given a time series $T$ the elements of the IMP are defined as follows for $i = 1, \dots, n - m + 1$,

$$IMP_i = \min_j \{d(T_{i,m}, T_{j,m}) \mid |i - j| \ge e \land \forall p \in \mathcal{P}, |\Pi_p(i - P_i) - j| \le [L/2]\},$$
where *e* represents the exclusion size, *P* is the list of periods (e.g., years), *P*<sub>*i*</sub> is the starting index of the period containing *i*, and *i* − *P*<sub>*i*</sub> gives the relative position of *i* within its period. The function Π<sub>*p*</sub> projects a relative position within the interval in a specific period *p*. For example, when using years as periods, if *i* represents January 10, 2010, then *P*<sub>*i*</sub> corresponds to January 1, 2010, and the relative position *i* − *P*<sub>*i*</sub> would be 10. The projection Π<sub>2018</sub>(*i* − *P*<sub>*i*</sub>) would map *i* to January 10, 2018<sup>1</sup>.

It is important to note that this definition does not account for the periodicity of the first and last years. For example, the last days of each year are not compared to the first days of the first year. Conversely, the first days of each year are not compared to the last days of the last year. This comparison will later be covered by triangular blocks in the blocking approach.

**Visual Illustration** To visually illustrate the difference between seasonal and interval for-
mulations, Figure 4 contrasts the block division of a single year, represented by season-based
squares, with the new division using fixed-size moving intervals. In the seasonal approach, the
first square (in orange) represents the first season, with its top left corner corresponding to the
first day of the year, while the second square (in blue) represents the second season, with the
bottom right corner marking the last day of the year. However, in the IMP formulation, each
day of the year becomes the center of its interval. As a result, moving from one day to the
next corresponds to shifting the interval's center by one day, effectively moving down one row
in the matrix. This shift produces parallelogram-shaped blocks. These parallelogram shapes
enable more pertinent comparisons, particularly at the beginning and end of the season, where
the search space transitions more smoothly. Additionally, the first days of the year are compared
with the last days if they fall within the same interval range, and vice versa. This leads to two
triangular-shaped regions, which are essentially truncated parallelograms, providing continuity
in the comparison process across year boundaries.

<figure><img src="image_7.png" /><figcaption>Figure 4: Comparison between seasonal and interval blocks partitioning of a year</figcaption></figure>

Building on the visual distinction between the seasonal and interval-based formulations, we now shift focus to the specific structure of the search space and its impact on the distance matrix. Figure 5 illustrates how the search space is constructed for a given subsequence in the IMP. The search space for the subsequence (in yellow) includes the subsequences (in blue) that lie within intervals of length *L*, centered on the projected positions of the subsequences in the second and third periods obtained by projections Π<sub>2</sub> and Π<sub>3</sub>. Figure 5b further demonstrates the partitioning of the corresponding distance matrix. Notably, in a periodic time series, subsequences like January 1 may need to be compared with those from December. This periodicity means that the interval covering the first days of a year also covers the last days of the last year. The height of the blocks corresponds to the length of a period (e.g., a year with 365 or 366 days), while the width represents the length of the interval *L*. The matrix is organized into rows of blocks, where

<footnote n="1">¹If the relative position of *i* is 366, corresponding to the last day of a leap year, it would be projected to the 365th day for non-leap years.</footnote>
each row covers the IMP for a specific year<fnref n="2" />. This symbolic matrix partitioning is well-suited to
a block algorithm.

<figure><img src="image_8.png" /><figcaption>Figure 5: Schematic representation of the Interval Matrix Profile</figcaption></figure>

**Blocking techniques** The IMP offers a computational advantage by enabling the exploitation of the distance matrix structure through matrix blocking, thereby enhancing data processing efficiency. In High-Performance Computing, matrix blocking is a well-established optimization strategy that has been demonstrated to enhance memory utilization and cache performance [30] allowing for significant performance boosts.

This blocking principle has previously been used in the SCAMP algorithm [17], which repre-
sents the current state of the art for the efficient computation of Matrix Profiles on distributed
GPUs, by partitioning the matrix into triangles and parallelograms of arbitrary sizes. The IMP
builds on this foundation by aligning its computations with the parallelogram block structure,
thereby ensuring compatibility with a block algorithm. In contrast to SCAMP, which employs
arbitrary partitioning, the IMP problem leads to an already partitioned matrix, as previously
described in Figure 5b, where the block dimensions depend on the problem’s parameters, i.e.
the interval length and the period length. In the IMP context, computing the coefficients of the
blocks can be done using the loop tiling [31] approach which consists of dividing the loops into
blocks (or tiles). For example, instead of using a naive approach that iterates over each row and
each column, loop tiling transforms these loops in order to have a loop iterating over the blocks
and one iterating within the blocks. This method has been developed to maximize parallelism
[32] and to reduce cache miss.

Moreover, by organizing the computations within these blocks, the IMP facilitates the vec-
torization of the distance calculations of Equation 2, as the block structure ensures that the
computations are performed on contiguous data. This structured partitioning not only enhances
computational efficiency, but also provides a scalable framework for parallelization, making the
IMP well-suited for either CPU or GPU parallelization using a task-based paradigm [33]. This
innovative adaptation provides an efficient algorithm that takes advantage of modern architec-
ture.

**Algorithm.** To compute the IMP, we propose the **Block IMP** (BIMP) algorithm. This algorithm takes advantage of the nature of the distance matrix and the blocking technique using

<footnote n="2">²Let us note that perfect alignment between blocks is not always possible. For example, leap years introduce a one-day shift, resulting in block heights of 365 or 366 days, depending on the year.</footnote>
<figure><img src="image_9.png" /><figcaption>Figure 6: Illustration of the block computations</figcaption></figure>

the aforementioned blocks. The first row inside each block is computed explicitly (Eq.1) to initialize the recurrence (Eq.2). Although the initialization process may be computationally expensive for large values of *m*, the advantage is twofold. Firstly, it allows each block to be computed independently; secondly, it provides greater numerical stability [17]. Figure 6 describes the computation organizations with the blocks. Figure 6a highlights the computation within a block shaped like a parallelogram. The goal is to produce local IMPs per block. Such local IMPs are obtained by computing the minimum per row of the block's coefficients. These intrablock computations are further described in Algorithm 1 in the case of a non-truncated parallelogram whose coordinates are given by its top left corner *i*, *j*. The distances contained in a row of a parallelogram are stored in a vector *v* of length *L*, while the local IMP is the vector of length *H*, where *H* is the height of the block, which can be either 365 or 366 for daily observations. Lines 3-4 describe the explicit computation of the distance. The nested loop of lines 6-7 describes the iteration over the rows and columns of the block. The loop over the columns (line 7) incrementally computes the distances using those of the previous rows, this loop is vectorized. Finally, the minimum of each row is written in the local IMP.

For the case of truncated (or partial) parallelograms, there are two types of truncation, right or left. For the right-truncated parallelograms, the number of columns to consider is not constant; therefore, the loops on the columns (lines 2 and 6) must be adapted. Especially since the total number of columns of the IMP matrix is *n* − *m* + 1, the loops must verify *j* + *c* ≤ *n* − *m* + 1. For the left-truncated blocks, the number of columns is also not constant, but this impacts the computations, as the recurrence must be initialized accordingly. There are two cases of left-truncated parallelograms polygons and triangles as described in Figure 7. The polygon has *L*/2
<table><caption>**Algorithm 1:** Block calculation routine</caption><thead><tr><th></th><th></th><th></th></tr></thead><tbody><tr><td>1</td><td>Initialize a vector <i>v</i> of length <i>L</i>;</td><td></td></tr><tr><td></td><td><code>/* Explicit computation */</code></td><td></td></tr><tr><td>2</td><td>for <i>c</i> = 0,...,<i>L</i>-1 do</td><td></td></tr><tr><td>3</td><td><code>  v<sub>c</sub> = d(T<sub>i,m</sub>, T<sub>j+c,m</sub>);</code></td><td></td></tr><tr><td>4</td><td>local IMP<sub>0</sub> ← min <i>v</i>;</td><td></td></tr><tr><td></td><td><code>/* Incremental computation */</code></td><td></td></tr><tr><td>5</td><td>for <i>r</i> = 1, ..., <i>H</i>-1 do</td><td></td></tr><tr><td>6</td><td><code>  for <i>c</i> = 0,...,<i>L</i>-1 do</code></td><td></td></tr><tr><td>7</td><td><code>    v<sub>c</sub> = v<sub>c</sub> - (T<sub>i+r-1</sub> - T<sub>j+c-1</sub>)<sup>2</sup> + (T<sub>i+r+m-1</sub> - T<sub>j+c+m-1</sub>)<sup>2</sup></code></td><td></td></tr><tr><td>8</td><td>local min ← min<sub>c</sub> <i>v</i>;</td><td></td></tr><tr><td>9</td><td><code>  local IMP<sub>r</sub> ← min(local min, local IMP<sub>r</sub>);</code></td><td></td></tr></tbody></table>

columns in the first row ($r = 0$) then $L/2 + 1$ in the second row ($r = 1$) and $L$ columns for $r \ge L/2$. When the number of columns is smaller than $L$, the first element of the row must be computed explicitly, and the others are computed incrementally. When the number of columns is $L$, all distances are computed incrementally using the loops of the routine. For the triangle, the first row has one element, the second two elements, and the last one $L/2$ elements. For each row of the triangle, the first element is computed explicitly, and the others are computed incrementally.

<figure><img src="image_10.png" /><figcaption>Figure 7: Schematic of the left-truncated parallelograms</figcaption></figure>

Figure 6b illustrates the process of computing the IMP from the local IMPs across partitioned blocks. The matrix is divided into symbolic blocks along each row, with each block yielding a local IMP of length $H$. To compute the final IMP, we perform an in-place reduction by iteratively taking the element-wise minimum across local IMPs within the row. Specifically, once the local IMP for a block is computed, it is passed to the next block in the row, which updates the local IMP in place rather than creating a new vector. This process continues until the last block in the row has been processed, at which point the final IMP is obtained.

The algorithm can be parallelized by adopting a task-based paradigm as follows. Each row of blocks is assigned to a separate task, which computes the block coefficients and updates the corresponding IMP for that row of blocks. Algorithm 2 summarizes our method for computing BIMP in parallel with a task-based paradigm, it highlights the loops to iterate over the blocks (lines 2 and 5). Furthermore, this parallel algorithm can also be applied to compute the SMP using the appropriate partitioning and blocks described in Section 4.1.
<caption>**Algorithm 2: Parallel BIMP**</caption>

**Input:** A time series $T$, a subsequence length $m$, the interval length $L$, the blocks' coordinates
**Output:** The IMP vector
1. Initialize the IMP;
2. **for** each row of block **do**
3. &nbsp;&nbsp;&nbsp;&nbsp;**Insert Task**
4. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Initialize a local IMP;
5. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**for** each block in row **do**
6. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Perform the blocked calculations with Algorithm 1;
7. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Write the local IMP to IMP;

**Complexity.** Let $P$ be the number of periods and $H$ be the height of the block, i.e. the period length. In each row of blocks, there are up $P+2$ blocks, with $P$ parallelograms and one or two triangles. With exactly $P$ rows of blocks, the total number of blocks to compute is approximately $O(P^2)$. Each block has a width of $L$, and the average height of the block is $H = n/P$, where $n$ is the time series length. Within each block, only the coefficients of the first row are explicitly computed, with a complexity of $O(mL)$. For climate-related applications, where $m \ll n/P$, the complexity to compute the coefficients of a single block is $O(Ln/P)$ leading to a total complexity for all block computations of $O(nPL)$. For a single row of blocks, there are at most $P+2$ local IMPs of length $n/P$, resulting in a reduction complexity of $O(n)$. Therefore, the total complexity of the reduction across all rows is $O(nP)$. To sum up, the total computational complexity to construct the IMP is $O(nPL)$. The algorithm stores the time series and the IMP of length $n$, the local IMP of length $H$, and a vector of length $L$ to perform the incremental computations leading to a space complexity of $O(n)$ as $n \gg H+L$. In the parallel configuration, the two vectors of length $H$ and $L$ are stored for each thread, we denote $n_t$ the number of threads, then the space complexity for the parallel version becomes $O(n+n_t(H+L))$

**Extension to the left and right IMP** The aforementioned example illustrates the computation of a full IMP. However, the approach can be adapted to compute alternative profiles, such as the left or right IMP. The study of left or right comparisons of distributions and patterns is particularly studied in climate under the prism of the *emergence* to identify new possible climate regimes [34, 35, 36].

The **left IMP** is defined as a variant of the IMP that focuses on detecting anomalies in the present by comparing each subsequence only to past observations. Formally, for a given index $i$, the left IMP is computed by finding the minimum distance between the subsequence starting at $i$ and all previous subsequences i.e. we consider all indices $j < i$. This approach helps to detect the first appearance of unusual subsequences within a seasonal interval, capturing anomalies that emerge in the current period compared to past subsequences.

To compute the left IMP efficiently with the BIMP algorithm, we restrict the computation to only the sub-diagonal blocks of the distance matrix. This ensures that each subsequence is only compared with subsequences from earlier periods, preserving the temporal directionality required for the left IMP. Additionally, we remove diagonal blocks to exclude comparisons within the same year, which eliminates trivial matches in the same period.

In contrast, a **right IMP** detects anomalies by comparing each subsequence at time $t$ with data from future periods, i.e. $j > i$. It is computed on the over-diagonal blocks of the matrix.
A discord detected by the right IMP reveals the last occurrence of a subsequence, which may
suggest the disappearance of a previously recurring pattern.

For these computations, a task-based parallelization strategy is particularly effective. Since
the yth period is compared only to the y − 1 previous periods, the computational load between
different periods (rows of blocks) becomes unbalanced. Task-based parallelization helps distribute
these computations efficiently, as tasks are dynamically assigned to threads.

This extension enhances the utility of the IMP by providing flexible options for anomaly
detection, whether focusing on past observations (left IMP) or future observations (right IMP).

4.3 Constrained kNN Interval Matrix Profile

When analyzing periodic time series with seasonal fluctuations, it is reasonable to expect that
a specific subsequence will appear in a significant number of periods. In contrast, anomalous
subsequences are likely to appear in only a few periods, often as isolated occurrences. The
standard IMP method effectively identifies anomalies that appear in only a single period, but
it may miss anomalies that recur sporadically across a few periods e.g., in k periods. This last
case is likely to occur when considering time series over a long historical record, in this context,
an anomaly can be defined as a subsequence that appears in at most k periods. For instance,
in climate applications, an abnormal event, such as a flash flood, may occur intermittently over
several years e.g., in two years among 60. An effective anomaly detection approach should
therefore consider subsequences that appear in no more than k distinct periods.

In such cases, the classical *k*NN MP could theoretically help to identify these rare occurrences
as it finds the *k*th nearest subsequence. However, this approach can mistakenly capture trivial
matches within the same year or trivially neighboring subsequences within the same period. For
example, if *k* = 3, the closest matches might be found within a small temporal window (e.g.,
January 5, 6, and 7 of the same year), which does not provide information on whether the
subsequence is genuinely rare across multiple periods. In such circumstances, it is impossible to
identify an anomaly as previously described.

To address this issue, we propose a **constrained kNN IMP (ckNN IMP)** which limits the *k* nearest neighbors to distinct intervals. This restriction is particularly useful for detecting abnormal climatic events that may recur in a few isolated periods.

**Algorithm.** To compute the *ck*NN IMP, we adapt the BIMP algorithm. In the BIMP algorithm, for each block we compute the minimum distance per row after finding the nearest neighbor in each block. In the case of *ck*NN IMP, instead of retaining the first minimum value, we keep the *k* minimum values using a maximum-heap structure of size *k*. This allows us to efficiently compute the *k*th minimum distance while ensuring that each nearest neighbor comes from a distinct interval.

Algorithm 3 illustrates the parallel *ck*NN BIMP, where max-heaps are used to manage the *k* nearest subsequences. Unlike the standard BIMP algorithm, which pipelines the minimum reduction across the blocks of a row, *ck*NN BIMP requires the minimum of each row of each block to determine the *k* nearest neighbor. To ensure this, the local IMP is reinitialized in line 6 to store the minimum distances for each row of the block. As described in lines 9–12, the minimum values for each row are stored in a max-heap of size *k*, which dynamically tracks the *k*th smallest distance. After all blocks in the row have been processed, the global value for the row is the root value of the corresponding max-heap. Finally, in lines 13–14, the *k*th minimum from the heap is written in the IMP vector, completing the computation for that row.
**Complexity.** Let $N$ denote the total number of blocks. For each block, finding the minimum among the $N$ blocks has a complexity of $O(N)$. When using a max-heap to track the $k$th minimum, the complexity to find the $k$th nearest neighbor among $N$ elements becomes $O(N \log k)$. Based on the complexity analysis of BIMP, the total reduction complexity with max-heaps is $O(nP \log k)$. Since $\log k \ll L$, the overall computational complexity remains $O(nLP)$. However, the space complexity increases to $O(n + H \log k)$, where $n$ is the size of the IMP vector and $H \log k$ represents the storage requirement for each heap of a row of blocks. For the parallel version, the space complexity becomes $O(n + n_t(H \log k + L))$, where $n_t$ is the number of threads.

<table><caption>**Algorithm 3: Parallel ckNN BIMP**</caption><thead><tr><th></th><th></th></tr></thead><tbody><tr><td>Input:</td><td>A time series T, subsequence length m, interval length L, nearest neighbor parameter k, and block coordinates</td></tr><tr><td>Output:</td><td>The ckNN IMP vector</td></tr><tr><td>1</td><td>Initialize the Matrix Profile;</td></tr><tr><td>2</td><td>for each row of block do</td></tr><tr><td>3</td><td> Insert Task</td></tr><tr><td>4</td><td> Initialize a vector v of max-heaps, each of size k;</td></tr><tr><td>5</td><td> for each block in row do</td></tr><tr><td>6</td><td> Initialize a local IMP;</td></tr><tr><td>7</td><td> Perform the blocked calculations with Algorithm 1;</td></tr><tr><td>8</td><td> for each local IMP<sub>i</sub> do</td></tr><tr><td>9</td><td> if IMP<sub>i</sub> &lt; v[i].top then</td></tr><tr><td>10</td><td> v[i].pop;</td></tr><tr><td>11</td><td> v[i].push(local IMP<sub>i</sub>);</td></tr><tr><td>12</td><td> for row index i do</td></tr><tr><td>13</td><td> ckNN IMP<sub>i</sub> &larr; v[i].top</td></tr></tbody></table>

**Advantages** The ckNN IMP approach enables the robust identification of abnormal but recurring anomalies across a few isolated periods. Unlike the classical kNN MP method, the ckNN IMP interval constraint ensures that all neighbors are meaningful and non-trivial. This constraint enables us to leverage a priori knowledge. For instance, a climatologist may define an anomaly as a subsequence that appears in at most 5 years, thereby setting $k$ to 6.

## 4.4 Finding the Disappearance of a Pattern using the Interval Matrix Profile

In many applications, particularly in climate analysis, it is valuable not only to detect anomalies but also to identify the disappearance of recurring patterns. For example, certain seasonal phenomena (e.g. annual floods [37]) may occur consistently over many years but then vanish abruptly. As previously mentioned, such focus could benefit *emergence* topics in climate-related sciences. To detect such disappearances, we propose a method that combines the left and right IMP with the ckNN IMP.

The core idea is as follows: if the right IMP at a given index $i$ has a high value, it indicates that the subsequence starting at $i$ has no close match in any subsequent periods, suggesting that it may no longer be present after this observation. Conversely, if the left IMP at the same
index *i* is low, it indicates that the subsequence was observed in previous periods. Especially, when computing a *ck*NN left IMP, if the distance value at index *i* is low it indicates that the subsequence has been observed in *k* periods. It allows to indicate that the subsequence was a pattern when looking at anterior periods but an anomaly when looking at later periods. If *k* is large enough, this pattern disappearance can be identified as a significant change or anomaly in the time series.

This process can be reversed, with the left IMP and the right *ck*NN IMP also being computing
This enables the identification of the initial occurrence of a subsequence. In other words, we can
also detect the emergence of a pattern.

We define the process of finding the disappearance of a pattern as follows.

1. Compute the right IMP of the time series. This profile reflects the minimum distance of each subsequence to any similar subsequence in future periods.

2. Compute the *ck*NN left IMP for each subsequence, with the constraint that neighbors are drawn from distinct previous periods up to *k* occurrences. This profile indicates the level of recurrence of each subsequence in *k* past periods.

3. Find an index *i* such that the right IMP at *i* is high (indicating no similar subsequence in future periods) and the *k*NN left IMP at *i* is low (indicating the subsequence was frequently observed in prior periods), then the subsequence at *i* is likely a pattern that has disappeared.

**Interpretation and Applications** By combining the left and right IMPs with the *ck*NN IMP, this approach enables robust detection of pattern disappearances across large time series (or pattern appearances). Such disappearances can reveal meaningful changes in time series behavior, particularly in climate or environmental data, where the sudden disappearance of a previously common pattern may signal shifts in environmental conditions or anthropogenic impacts. For example, if a seasonal weather pattern or temperature fluctuation consistently appears over several decades but abruptly vanishes, this may indicate a notable climate anomaly or transition.

This method of detecting the disappearance of patterns provides an additional layer of knowl-
edge extraction beyond classical anomaly detection, helping researchers and practitioners identify
significant shifts in time series patterns that could have critical implications for long-term analysis
and forecasting.

5 Experimental Evaluation

In this section, we evaluate the performance of our proposed algorithms. First, we outline
the experimental setup. Next, we analyze the computational efficiency of our algorithms for
calculating the IMP and the *ck*NN IMP with respect to the different parameters of the algorithms.
Then, we study the numerical precision of the incremental algorithms. We also demonstrate
the limitations of the classical *k*NN MP technique through a specific example, highlighting its
susceptibility to trivial matches. Finally, we assess the practical utility of our methods using
climate time series, showcasing their effectiveness in identifying meaningful anomalies. Our
analysis is limited to the IMP, the generic and robust iteration of the SMP.
<figure><img src="image_11.png" /><figcaption>Figure 8: Execution time (in seconds) of BIMP, modified AAMP, and Brute Force algorithms as a function of (a) subsequence length, (b) interval length, (c) number of periods, and (d) period length. The default parameters are *m* = 7, *L* = 90, 100 periods, and a period length of 365</figcaption></figure>

## 5.1 Setup

Our different algorithms are implemented<fnref n="3" /> in C++20 and binded in Python using Pybind11 [38], the task parallelization is made with OpenMP. The experiments are conducted on an Intel i9-13900H @2.6GHz laptop with turbo boost disabled. The parallel experiments are performed on a 2x18-core Cascade Lake Intel Xeon Skylake Gold 6240@2.6 GHz (PlaFRIM testbed).

In our experiments, we use both synthetic and real-world data. For runtime experiments, we use two synthetic time series with uniform distributions, consisting of two types of time series: one with daily observations and the other with hourly observations. For climate case studies, we investigate two climate time series time series obtained from *Copernicus Climate Pulse*<fnref n="4" />[11, 12]. The first one represents the daily average of the global air surface temperature (TAS) of the Earth from 1940 to 2022. The second one represents the daily average sea surface temperature (SST) of the earth without the poles (*60°S – 60°N*) from 1979 to October 2024.

We compare the performance of our algorithms with two others: 1) Brute force algorithm that computes explicitly (Eq.1) all coefficients belonging to an interval; 2) Modified AAMP that incrementally computes all the coefficients of a row, then keeps the minimum distance belonging to an interval. Both algorithms are described in Section 2. For the Brute force implementation, the distance calculation is vectorized, and for the Modified AAMP the whole row calculation is vectorized.

<footnote n="3">³The code is available at https://github.com/GuillaumeCld/interval_matrix_profile</footnote>

<footnote n="4">⁴https://pulse.climate.copernicus.eu</footnote>
## 5.2 Performance Evaluation

We show in Figure 8 the impact of the algorithm parameters on the execution time on the synthetic time series with 100 years of daily observation. Figure 8a shows the execution time of BIMP, Brute force, and Modified AAMP by varying the length of the subsequence *m*. As expected, BIMP and modified AAMP are not impacted by the parameter *m*.

Figure 8b illustrates the execution time by varying the interval length *L*. We use the daily observations. All the algorithms' computation times increase linearly. Once more, the brute force scales poorly. The overhead of modified AMMP compared to BIMP remains constant. However, as *L* gets closer to a year's length (365), the closer the BIMP performance is to a classic AAMP as the search space covers almost all the time series.

Figure 8c highlights the impact of the number of periods (i.e., the number of years) on execution time. As seen, the runtime of all algorithms increases quadratically to the number of years. However, since modified AAMMP has a significant overhead compared to BIMP, it scales poorly. In this experiment, as *m* is small (*m* = 7), the brute force algorithm has a similar runtime to modified AAMP, but it would be worse for larger *m* values.

Figure 8d depicts the execution time as a function of the period length with *m* = 7, *L* = 90, and 100 periods. For the IMP algorithms, as *L* is fixed, the time grows linearly, while it grows quadratically for the AAMP. When the period length is 900, the IMP reduces that amount of calculations to 10% of the original MP.

In all four experiments, the BIMP algorithm consistently outperforms the modified AAMP in terms of speed. This performance gap becomes especially significant as the number of periods or the period length increases.

We also investigate the impact of the number of neighbors *k* on the *ckNN* BIMP algorithm compared to the BIMP algorithm. The experimentation is conducted on 100 years of synthetic hourly data with *m* = 168 and *L* = 2160. Table 1 shows that for *k* = 1, *ckNN* BIMP has a small overhead induced by the reduction step and the use of heaps. However, the value of *k* has little impact on the execution time.

<table><caption>Table 1: Execution times of the *ckNN* BIMP algorithm and BIMP for different values of *k* for 100 years of hourly data with *m* = 168 and *L* = 2160</caption><thead><tr><th></th><th>BIMP</th><th colspan="4">ckNN BIMP</th></tr></thead><tbody><tr><th>k</th><td>/</td><td>1</td><td>5</td><td>10</td><td>20</td></tr><tr><th>Time (min)</th><td>2.99</td><td>3.01</td><td>3.02</td><td>3.03</td><td>3.18</td></tr></tbody></table>

In Figure 9, we study the parallelization of the BIMP algorithm from 1 to 30 cores, using the hourly time series for 200 years (*n* = 1,752,000), *m* = 168 and *L* = 2160. We plot side by side the execution time and the speedup. The speedup is the ratio of the sequential time to the parallel time. The speedup is almost linear, illustrating that the algorithm is embarrassingly parallel. However, it is important to note that this relationship is not strictly linear. Instead, it can be described as a step function, as the level of parallelism depends on the number of years. Indeed, when the number of threads divides the number of years, the parallelism level is maximal as the computations are optimally distributed across the cores. Otherwise, the computations are not evenly distributed.
<figure><img src="image_12.png" /><figcaption>Figure 9: BIMP execution time (left axis) and corresponding speedup (right axis) as a function of the number of threads</figcaption></figure>

## 5.3 Precision Evaluation

As the algorithms use incremental calculations, it is important to evaluate their precision. To evaluate the numerical precision, we measure the absolute and relative error. Let $\hat{x}$ be the approximation of a true vector $x$, then the absolute error is $\|x - \hat{x}\|_2$ and the relative error is $\|x - \hat{x}\|_2/\|x\|_2$. For this evaluation $x$ is the IMP computed using the brute force method in double precision, and $\hat{x}$ is the IMP computed with both modified AAMP and BIMP in either single or double precision. Table 2 shows the absolute and relative error computed on a time series of daily temperature observations with 74 years, $m = 7$ and $L = 90$. The BIMP method is systematically an order of magnitude lower than the modified AAMP, indicating the benefits of re-initializing the incremental computations at the start of each block. When computing the IMP with a single precision, BIMP allows for the reduction of numerical errors. However, the absolute error is still high ($5.1 \times 10^{-3}$), the error comes from the incremental distance computation described in Equation 2. The recurrence allows to incrementally compute the element of a diagonal, where the diagonal length is the block height $H$, the second diagonal element is obtained after 1 iteration, and the last diagonal element is obtained after $H - 1$ iterations. The computation of the last element is equivalent to a recursive summation of length $H - 1$. This type of summation is known to lead to numerical errors [16] and the error is conditioned by the values that are manipulated i.e. the distance values and the time series values. Nevertheless, in double precision, this error is negligible for this application, but this increases the cost of calculations and double the use of memory, which can be a significant disadvantage.

<table><caption>Table 2: Absolute and relative error comparison between modified AAMP and BIMP for daily temperature observations with 74 years, $m = 7$ and $L = 90$</caption><thead><tr><th rowspan="2">Precision</th><th colspan="2">Modified AAMP</th><th colspan="2">BIMP</th></tr><tr><th>Single</th><th>Double</th><th>Single</th><th>Double</th></tr></thead><tbody><tr><td>Absolute Error</td><td>3.7 10<sup>-2</sup></td><td>6.8 10<sup>-11</sup></td><td>5.1 10<sup>-3</sup></td><td>9.4 10<sup>-13</sup></td></tr><tr><td>Relative Error</td><td>1.2 10<sup>-4</sup></td><td>2.3 10<sup>-13</sup></td><td>1.7 10<sup>-5</sup></td><td>3.2 10<sup>-14</sup></td></tr></tbody></table>
## 5.4 Trivial Neighbors in the Classical kNN MP

To illustrate the limitation of the classical kNN MP we conduct the following experiment: we compute, on the air temperature (TAS) time series, the kNN IMP without the constraint for $k \in [1, 10]$ with $m = 14$ and $L = 61$. For $k = 5$ and $k = 10$, we analyze whether the $k$th neighbor was a trivial neighbor of one of the preceding neighbors, defined by the exclusion zone $e = 7$. The results are summarized in Table 3, showing the number and percentage of trivial versus non-trivial neighbors for the 5th and 10th nearest neighbors. The results highlight that a significant proportion of the neighbors found using the classical kNN IMP is trivial, with trivial neighbors accounting for 38% of the 5th neighbors and 49% of the 10th neighbors. With the ckNN IMP, it is guaranteed to have non-trivial neighbors.

<table><caption>Table 3: Number (and percentage) of trivial and non-trivial neighbors among the 5th and 10th neighbors</caption><thead><tr><td>Neighbor</td><td>Non-trivial</td><td>Trivial</td></tr></thead><tbody><tr><td>5</td><td>18,925 (62%)</td><td>11,378 (38%)</td></tr><tr><td>10</td><td>15,364 (51%)</td><td>14,939 (49%)</td></tr></tbody></table>

## 5.5 Climate Case Study

We now evaluate the utility of our proposed methods on real-world climate data, demonstrating their effectiveness in detecting seasonal anomalies and pattern changes.

For the first time series (TAS) we compute both the MP and the IMP with $m = 14$ and $L = 61$. Figure 10a shows the MP (in blue) and the IMP (in red) values on November 6, 2024, where a discord is detected using the IMP. The distance value of this discord is higher than the top-1 discord detected with MP, which is indicated by the dashed blue line. In Figure 10b, we see the discord (in black) and the nearest neighbors found by the MP (in blue) and by the IMP (in red). The discord corresponds to a 14-day subsequence in November where the temperature rises, contrary to the expected decrease for that time of year. The nearest neighbor identified by the MP occurs in March and April, a period when temperatures naturally increase. However, the interval restriction imposed by IMP limits the search to subsequences between October and early December. As a result, no similar subsequences are found, thus the subsequence is a discord.

For the second time series (SST), we compare the MP, the IMP, and the ckNN IMP with $m = 7$, $L = 61$ and $k \in \{3, 5, 10\}$. In Figure 11a, we display the daily SST for different years, with a specific focus on the abnormal years 2023 (in red) and 2024 (in green) which highlight the need to provide constraints to the MP. Indeed, the month of August 2023, is similar to August 2024 which is also abnormal, but it is also similar to March in some other years. Additionally, the two last years stand out against the past years, making them ideal candidates for anomaly detection. Figure 11b shows the distance values calculated for 2023 and 2024 using the MP, the IMP, and the ckNN IMP, allowing us to assess how each method detects anomalies. The dashed black line represents the top-1 discord value identified by the MP indicating that the MP does not find any anomalies in the last two years. In particular, the subsequences from August have neighbors in March and April, making it difficult to identify truly anomalous behavior across seasons. By applying the IMP, we can avoid such misleading comparisons between distant parts of the year. While the IMP identifies some discords, particularly in January 2024, most subsequences have close neighbors in 2023 or 2024. Especially, August 2023 is close to August 2024 and is not detected as abnormal. This is problematic because all subsequences in these two months
<figure><img src="image_13.png" /><figcaption>(a) Zoomed comparison of the Matrix Profile and Interval Matrix Profile for the TAS</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(b) Zoomed representation of the discord and the closest match for the MP and IMP</figcaption></figure>

<caption>Figure 10: Illustration of a top 1 discord discovered with IMP</caption>

exhibit abnormal behavior and the IMP does not identify these anomalies. In contrast, the *ckNN* IMP method introduces a constraint that the anomalous subsequence must occur in *k* different years. This makes the *ckNN* IMP highly sensitive to detect widespread anomalies, such as those observed in August of both 2023 and 2024. The results show that the *ckNN* IMP captures the level of abnormality more effectively, as reflected in the distance values. For example, the top-1 discord identified by the MP has a distance value of 0.19, whereas the top-1 discord detected by the c3NN IMP has a distance value of 0.79, i.e. more than four times higher, illustrating a
<figure><img src="image_15.png" /><figcaption>(b) Zoomed comparison of the MP, IMP and ckNN IMP for $k = 3, 5, 10$ of the distance values for 2023 and 2024</figcaption></figure>

<caption>Figure 11: Application to climate pulse daily Sea Surface Temperature of the Earth [60°S–60°N]</caption>

stronger and clearer detection of anomalies. The ckNN method also reveals the specific timing of changes in distances (e.g., increasing from July to September 2023), highlighting periods with stronger anomalies relative to the interval dates. This information could enable climate scientists to focus their analyses on these critical time periods, assessing potential causes, impacts, and connections to related extreme events.

We now explore the methodology to find the disappearance of a pattern in the SST time series. We compute the right IMP and the ckNN left IMP with $k = 10$ and remove the first 10 years. This filtering is mandatory, as it is impossible to compute the 10th neighbor for the first 10 years. We also remove the last 10 years to compare the anomalies discovered by the right IMP
to at least 10 other years. In figure 12, we display the c10NN left IMP (in blue), the right IMP (in red), and the right discord index by a black dot. The discord verifies the condition to find the last occurrence of a pattern, i.e. a high right IMP value and a low ckNN left IMP value. The relatively high distances in the ckNN left IMP combined with the low distances in the right IMP up to 2023 indicate an unprecedented anomaly when compared with historical patterns. This approach offers a valuable indicator and method for identifying the emergence and abruptness of changes and anomalies in climate-related data [35]. In Figure 13, we show this last occurrence in black, we display the 10 nearest constrained neighbors for the left IMP in blue and the right IMP in red. The neighbors on the left all have similar shapes and amplitude, indicating that a similar subsequence is present for at least 10 years. However, the right neighbors all have different amplitudes indicating that after the last occurrence in 2011, the subsequence is never seen again. Additionally, the ckNN left IMP values have an important variability, which is consistent with climate-related literature about increasing temperature variability with global warming [39].

<figure><img src="image_16.png" /><figcaption>Figure 12: Zoomed comparison between the right IMP and c10NN left IMP on the daily sea surface temperature time series</figcaption></figure>

# 6 Conclusion

In this work, we investigated extensions to the classical Matrix Profile for periodic and seasonal time series, with a focus on climate applications. We introduced the Seasonal Matrix Profile, a naive approach to take into account seasonality. Then, we proposed the Interval Matrix Profile to allow more flexible comparisons of subsequences within similar date ranges, addressing the limitations of abrupt seasonal boundaries in the Seasonal Matrix Profile.

We also developed the Block IMP algorithm for efficient and scalable computation and extended it to a constrained kNN IMP algorithm to detect abnormal events occurring across different periods. Our experiments illustrated that the IMP and the ckNN IMP outperform the MP in identifying climate anomalies, offering more robust and interpretable results. In particular, the IMP detects anomalies missed by the MP when they are abnormal subsequences within their
<figure><img src="image_17.png" /><figcaption>Figure 13: Illustration of a pattern (in blue) on 10 different years and its last occurrence (in black). Red subsequences are the 10 right nearest neighbors of the black subsequence, blue subsequences are the 10 left nearest neighbors</figcaption></figure>

seasons, while the *ckNN* IMP identifies abnormal patterns present in only a few periods. Finally, we demonstrated how the aforementioned methods can be effectively employed to identify the disappearance of a seasonal pattern.

Future work will explore the application of the IMP and *ckNN*-IMP methodologies to other domains with periodic and seasonal data, such as healthcare, retail, and manufacturing. Additionally, we plan to investigate the integration of other distances such as Dynamic Time Warping within the IMP technique to further enhance its ability to capture seasonal anomalies and patterns. Finally, other search restriction approaches can be explored, like dynamic interval length or improved seasonal restriction.

## Acknowledgments

This work is supported by the Région Occitanie. Parallel experiments presented in this paper were carried out using the PlaFrim experimental testbed, supported by Inria, CNRS (LABRI and IMB), Université de Bordeaux, Bordeaux INP, and Conseil Régional d'Aquitaine.

## References

[1] Sebastian Schmidl, Phillip Wenig, and Thorsten Papenbrock. Anomaly detection in time series: A comprehensive evaluation. *Proceedings of the VLDB Endowment*, 15(9):1779–1797, May 2022. doi: 10.14778/3538598.3538602.

[2] Wenjie Yang, Ruofan Wang, and Bofan Wang. Detection of Anomaly Stock Price Based
on Time Series Deep Learning Models. In *2020 Management Science Informatization and Economic Innovation Development Conference (MSIEID)*, pages 110–114, December 2020. doi: 10.1109/MSIEID52046.2020.00029.

[3] Xue Yang, Xuejun Qi, and Xiaobo Zhou. Deep Learning Technologies for Time Series Anomaly Detection in Healthcare: A Review. *IEEE Access*, 11:117788–117799, 2023. doi: 10.1109/ACCESS.2023.3325896.

[4] Peng Yan, Ahmed Abdulkadir, Paul-Philipp Luley, Matthias Rosenthal, Gerrit A. Schatte, Benjamin F. Grewe, and Thilo Stadelmann. A Comprehensive Survey of Deep Transfer Learning for Anomaly Detection in Industrial Time Series: Methods, Applications, and Directions. *IEEE Access*, 12:3768–3789, 2024. doi: 10.1109/ACCESS.2023.3349132.

[5] Yang Chen, Wilfran Moufouma-Okia, Valérie Masson-Delmotte, Panmao Zhai, and Anna Pirani. Recent Progress and Emerging Topics on Weather and Climate Extremes Since the Fifth Assessment Report of the Intergovernmental Panel on Climate Change. *Annual Review of Environment and Resources*, 43(1):35–59, 2018. doi: 10.1146/annurev-environ-102017-030052.

[6] Sancho Salcedo-Sanz, Jorge Pérez-Aracil, Guido Ascenso, Javier Del Ser, David Casillas-Pérez, Christopher Kadow, Dušan Fister, David Barriopedro, Ricardo García-Herrera, Matteo Giuliani, and Andrea Castelletti. Analysis, characterization, prediction, and attribution of extreme atmospheric events with machine learning and deep learning techniques: A review. *Theoretical and Applied Climatology*, 155(1):1–44, January 2024. doi: 10.1007/s00704-023-04571-5.

[7] Ane Blázquez-García, Angel Conde, Usue Mori, and Jose A. Lozano. A Review on Outlier/Anomaly Detection in Time Series Data. *ACM Computing Surveys*, 54(3):56:1–56:33, April 2021. doi: 10.1145/3444690.

[8] Omar Alghushairy, Raed Alsini, Terence Soule, and Xiaogang Ma. A Review of Local Outlier Factor Algorithms for Outlier Detection in Big Data Streams. *Big Data and Cognitive Computing*, 5(1):1, March 2021. doi: 10.3390/bdcc5010001.

[9] Hamzeh Alimohammadi and Shengnan Nancy Chen. Performance evaluation of outlier detection techniques in production timeseries: A systematic review and meta-analysis. *Expert Systems with Applications*, 191:116371, April 2022. doi: 10.1016/j.eswa.2021.116371.

[10] Chin-Chia Michael Yeh, Yan Zhu, Liudmila Ulanova, Nurjahan Begum, Yifei Ding, Hoang Anh Dau, Diego Furtado Silva, Abdullah Mueen, and Eamonn Keogh. Matrix Profile I: All Pairs Similarity Joins for Time Series: A Unifying View That Includes Motifs, Discords and Shapelets. In *2016 IEEE 16th International Conference on Data Mining (ICDM)*, pages 1317–1322, December 2016. doi: 10.1109/ICDM.2016.0179.

[11] Hans Hersbach, Bill Bell, Paul Berrisford, Gionata Biavati, András Horányi, Joaquín Muñoz Sabater, Julien Nicolas, Carole Peubey, Raluca Radu, Iryna Rozum, Dinand Schepers, Adrian Simmons, Cornel Soci, Dick Dee, and Jean-Noël Thépaut. ERA5 hourly data on single levels from 1940 to present. Copernicus Climate Change Service (C3S) Climate Data Store (CDS), 2023. doi: 10.24381/cds.adbb2d47 (Accessed on 10-06-2024).

[12] Bill Bell, Hans Hersbach, Adrian Simmons, Paul Berrisford, Per Dahlgren, András Horányi, Joaquín Muñoz-Sabater, Julien Nicolas, Raluca Radu, Dinand Schepers, Cornel Soci, Sebastien Villaume, Jean-Raymond Bidlot, Leo Haimberger, Jack Woollen, Carlo Buontempo,
and Jean-Noël Thépaut. The ERA5 global reanalysis: Preliminary extension to 1950. *Quarterly Journal of the Royal Meteorological Society*, 147(741):4186-4227, October 2021. doi: 10.1002/qj.4174.

[13] Tanmoy Mondal, Reza Akbarinia, and Florent Masseglia. kNN matrix profile for knowledge discovery from time series. *Data Mining and Knowledge Discovery*, 37(3):1055-1089, May 2023. doi: 10.1007/s10618-022-00883-8.

[14] Reza Akbarinia and Bertrand Cloez. Efficient Matrix Profile Computation Using Different Distance Functions, January 2019.

[15] Yan Zhu, Zachary Zimmerman, Nader Shakibay Senobari, Chin-Chia Michael Yeh, Gareth Funning, Abdullah Mueen, Philip Brisk, and Eamonn Keogh. Matrix Profile II: Exploiting a Novel Algorithm and GPUs to Break the One Hundred Million Barrier for Time Series Motifs and Joins. In *2016 IEEE 16th International Conference on Data Mining (ICDM)*, pages 739-748, December 2016. doi: 10.1109/ICDM.2016.0085.

[16] Nicholas J. Higham. *Accuracy and Stability of Numerical Algorithms*. Society for Industrial and Applied Mathematics, Philadelphia, 2nd edition, 2002.

[17] Zachary Zimmerman, Kaveh Kamgar, Nader Shakibay Senobari, Brian Crites, Gareth Funning, Philip Brisk, and Eamonn Keogh. Matrix Profile XIV: Scaling Time Series Motif Discovery with GPUs to Break a Quintillion Pairwise Comparisons a Day and Beyond. In *Proceedings of the ACM Symposium on Cloud Computing*, SoCC '19, pages 74-86, New York, NY, USA, November 2019. Association for Computing Machinery. doi: 10.1145/3357223.3362721.

[18] Nader Shabikay Senobari, Peter M. Shearer, Gareth J. Funning, Zachary Zimmerman, Yan Zhu, Philip Brisk, and Eamonn Keogh. The Matrix Profile in Seismology: Template Matching of Everything With Everything. *Journal of Geophysical Research: Solid Earth*, 129(2): e2023JB027122, 2024. doi: 10.1029/2023JB027122.

[19] Hussein El Khansa, Carmen Gervet, and Audrey Brouillet. Prominent Discord Discovery with Matrix Profile: Application to Climate Data Insights. In *Computer Science & Technology Trends*, pages 65-79. Academy and Industry Research Collaboration Center (AIRCC), May 2022. doi: 10.5121/csit.2022.120806.

[20] Colin Hehir and Alan F. Smeaton. Calculating the matrix profile from noisy data. *PLOS ONE*, 18(6), June 2023. doi: 10.1371/journal.pone.0286763.

[21] Jose Carlos Romero, Angeles Navarro, Antonio Vilches, Andrés Rodríguez, Francisco Corbera, and Rafael Asenjo. Efficient heterogeneous matrix profile on a CPU + High Performance FPGA with integrated HBM. *Future Generation Computer Systems*, 125:10-23, December 2021. doi: 10.1016/j.future.2021.06.025.

[22] Sara Alaee, Kaveh Kamgar, and Eamonn Keogh. Matrix Profile XXII: Exact Discovery of Time Series Motifs Under DTW. In *2020 IEEE International Conference on Data Mining (ICDM)*, pages 900-905, November 2020. doi: 10.1109/ICDM50108.2020.00099.

[23] Sadaf Tafazoli, Yue Lu, Renjie Wu, Thirumalai Vinjamoor Akhil Srinivas, Hannah Dela Cruz, Ryan Mercer, and Eamonn Keogh. C22MP: The marriage of catch22 and the matrix profile creates a fast, efficient and interpretable anomaly detector. *Knowledge and Information Systems*, 66:4789-4823, May 2024. doi: 10.1007/s10115-024-02107-5.
[24] Yue Lu, Renjie Wu, Abdullah Mueen, Maria A. Zuluaga, and Eamonn Keogh. DAMP: Accurate time series anomaly detection on trillions of datapoints and ultra-fast arriving data streams. *Data Mining and Knowledge Discovery*, 37(2):627-669, March 2023. doi: 10.1007/s10618-022-00911-7.

[25] Yue Lu, Thirumalai Vinjamoor Akhil Srinivas, Takaaki Nakamura, Makoto Imamura, and Eamonn Keogh. Matrix Profile XXX: MADRID: A Hyper-Anytime and Parameter-Free Algorithm to Find Time Series Anomalies of all Lengths. In *2023 IEEE International Conference on Data Mining (ICDM)*, pages 1199-1204, December 2023. doi: 10.1109/ICDM58522.2023.00148.

[26] Tian Huang, Yongxin Zhu, Yafei Wu, and Weiwei Shi. J-Distance Discord: An Improved Time Series Discord Definition and Discovery Method. In *2015 IEEE International Conference on Data Mining Workshop (ICDMW)*, pages 303-310, November 2015. doi: 10.1109/ICDMW.2015.120.

[27] Thuy Thi Thu Huynh, Anh Tuan Duong, and Chau Thi Ngoc Vo. A new discord definition and an efficient time series discord detection method using GPUs. In *Proceedings of the 2021 3rd International Conference on Software Engineering and Development*, ICSED '21, pages 63-70, New York, NY, USA, February 2022. Association for Computing Machinery. doi: 10.1145/3507473.3507483.

[28] Yan Zhu, Chin-Chia Michael Yeh, Zachary Zimmerman, and Eamonn Keogh. Matrix Profile XVII: Indexing the Matrix Profile to Allow Arbitrary Range Queries. In *2020 IEEE 36th International Conference on Data Engineering (ICDE)*, pages 1846-1849, April 2020. doi: 10.1109/ICDE48307.2020.00185.

[29] Hoang Anh Dau and Eamonn Keogh. Matrix Profile V: A Generic Technique to Incorporate Domain Knowledge into Motif Discovery. In *Proceedings of the 23rd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining*, KDD '17, pages 125-134, New York, NY, USA, August 2017. Association for Computing Machinery. doi: 10.1145/3097983.3097993.

[30] Monica D. Lam, Edward E. Rothberg, and Michael E. Wolf. The cache performance and optimizations of blocked algorithms. *ACM SIGOPS Operating Systems Review*, 25(Special Issue):63-74, April 1991. doi: 10.1145/106974.106981.

[31] Jingling Xue. On Tiling as a Loop Transformation. *Parallel Processing Letters*, 07(04): 409-424, December 1997. doi: 10.1142/S0129626497000401.

[32] Jingling Xue. *Loop Tiling for Parallelism*, volume 575 of Springer Science & Business Media. Springer Science & Business Media, 2000.

[33] Cédric Augonnet. *Scheduling Tasks over Multicore Machines Enhanced with Accelerators: A Runtime System's Perspective*. These de doctorat, Bordeaux 1, December 2011.

[34] Andrew D. King, Markus G. Donat, Erich M. Fischer, Ed Hawkins, Lisa V. Alexander, David J. Karoly, Andrea J. Dittus, Sophie C. Lewis, and Sarah E. Perkins. The timing of anthropogenic emergence in simulated climate extremes. *Environmental Research Letters*, 10(9):094015, September 2015. doi: 10.1088/1748-9326/10/9/094015.

[35] Marco Gaetani, Serge Janicot, Mathieu Vrac, Adjoua Moise Famien, and Benjamin Sultan. Robust assessment of the time of emergence of precipitation change in West Africa. *Scientific Reports*, 10(1):7670, May 2020. doi: 10.1038/s41598-020-63782-2.
[36] Eric Pohl, Christophe Grenier, Mathieu Vrac, and Masa Kageyama. Emerging climate signals in the Lena River catchment: A non-parametric statistical approach. *Hydrology and Earth System Sciences*, 24(5):2817-2839, May 2020. doi: 10.5194/hess-24-2817-2020.

[37] J. A. Marengo and J. C. Espinoza. Extreme seasonal droughts and floods in Amazonia: Causes, trends and impacts. *International Journal of Climatology*, 36(3):1033-1050, 2016. doi: 10.1002/joc.4420.

[38] Wenzel Jakob, Jason Rhinelander, and Dean Moldovan. pybind11 – seamless operability between c++11 and python, 2017. https://github.com/pybind/pybind11.

[39] M. James Salinger. Climate Variability and Change: Past, Present and Future – An Overview. *Climatic Change*, 70(1):9-29, May 2005. doi: 10.1007/s10584-005-5936-x.
Inria

RESEARCH CENTRE
Centre Inria d'Université Côte d'Azur

2004 route des Lucioles - BP 93
06902 Sophia Antipolis Cedex

Publisher
Inria
Domaine de Voluceau - Rocquencourt
BP 105 - 78153 Le Chesnay Cedex
inria.fr

ISSN 0249-6399