How to properly use the "discords" parameter in mstump for anomaly detection?

[dbolotov]

Hello,

Thank you for the very useful package! I am trying to adapt the doc example for mstump to find discords instead of motifs in multidimensional data.

• Is it necessary to set discords=True in stumpy.mstump() when looking for discords?
• From the doc for the discords parameter: "When set to True, this reverses the distance matrix which results in a multi-dimensional matrix profile that favors larger matrix profile values (i.e., discords) rather than smaller values (i.e., motifs).". What exactly is meant by reversing the distance matrix - is it flipped somehow?
• If discords=True, can the resulting MP still be used to find motifs?

Including code example below (reproduced from the documentation with changes to find discords). I first set discords=False in stumpy.stump(), and calculated the 3 top discord indices using discords_idx = np.argsort(-mps, axis=1)[:,:3]. Then I plugged discords_idx into where motifs_idx is used when finding the best k with the elbow method and the subspace. Then repeated the whole thing with discords=True. Also including two sets of plots - one where discords=False, and one where discords=True. The green vertical lines represent the top 3 discords for each k, and green segments show the area corresponding to the top discord for k=1.

Thank you.

import pandas as pd
import numpy as np
import stumpy
import matplotlib.pyplot as plt
import os

plt.style.use('https://raw.githubusercontent.com/TDAmeritrade/stumpy/main/docs/stumpy.mplstyle')

# Loading and Visualizing the Toy Data
df = pd.read_csv("https://zenodo.org/record/4328047/files/toy.csv?download=1")

# Multi-dimensional Matrix Profiles
m = 30
discords_md = True
mps, indices = stumpy.mstump(T = df, m = m, discords = discords_md, normalize = True)


# find k-dimensional motif
motifs_idx = np.argsort(mps, axis=1)[:, :2]

# find k-dimensional discords
discords_idx = np.argsort(-mps, axis=1)[:,:3] # index of top 3 discords for each k

fig, axs = plt.subplots(mps.shape[0] * 2, sharex=True, gridspec_kw={'hspace': 0})

for k, dim_name in enumerate(df.columns):
    axs[k].set_ylabel(dim_name, fontsize='20')
    axs[k].plot(df[dim_name])
    axs[k].set_xlabel('Time', fontsize ='20')

    
    axs[k + mps.shape[0]].set_ylabel(dim_name.replace('T', 'P'), fontsize='20')
    axs[k + mps.shape[0]].plot(mps[k], c='orange')
    axs[k + mps.shape[0]].set_xlabel('Time', fontsize ='20')

    axs[k].axvline(x=motifs_idx[1, 0], linestyle="dashed", c='black')
    axs[k].axvline(x=motifs_idx[1, 1], linestyle="dashed", c='black')
    axs[k + mps.shape[0]].axvline(x=motifs_idx[1, 0], linestyle="dashed", c='black')
    axs[k + mps.shape[0]].axvline(x=motifs_idx[1, 1], linestyle="dashed", c='black')
    
    if dim_name != 'T3':
        axs[k].plot(range(motifs_idx[k, 0], motifs_idx[k, 0] + m), df[dim_name].iloc[motifs_idx[k, 0] : motifs_idx[k, 0] + m], c='red', linewidth=4)
        axs[k].plot(range(motifs_idx[k, 1], motifs_idx[k, 1] + m), df[dim_name].iloc[motifs_idx[k, 1] : motifs_idx[k, 1] + m], c='red', linewidth=4)
        axs[k + mps.shape[0]].plot(motifs_idx[k, 0], mps[k, motifs_idx[k, 0]] + 1, marker="v", markersize=10, color='red')
        axs[k + mps.shape[0]].plot(motifs_idx[k, 1], mps[k, motifs_idx[k, 1]] + 1, marker="v", markersize=10, color='red')
    else:
        axs[k + mps.shape[0]].plot(motifs_idx[k, 0], mps[k, motifs_idx[k, 0]] + 1, marker="v", markersize=10, color='black')
        axs[k + mps.shape[0]].plot(motifs_idx[k, 1], mps[k, motifs_idx[k, 1]] + 1, marker="v", markersize=10, color='black')

    # mark the top 3 discords for each k with vertical line on the MP plots
    axs[k + mps.shape[0]].axvline(x=discords_idx[k, 0], linestyle="dashed", c='green')
    axs[k + mps.shape[0]].axvline(x=discords_idx[k, 1], linestyle="dashed", c='green')
    axs[k + mps.shape[0]].axvline(x=discords_idx[k, 2], linestyle="dashed", c='green')

    # plot the k=1 discord over the original data
    axs[k].plot(range(discords_idx[1, 0], discords_idx[1, 0] + m), df[dim_name].iloc[discords_idx[1, 0] : discords_idx[1, 0] + m], c='green', linewidth=4)


plt.suptitle('multidimensional MP with discords={0}'.format(discords_md), size=20)
plt.show()

# Finding the right k
plt.plot(mps[range(mps.shape[0]), discords_idx[:, 0]], c='red', linewidth='4')
plt.xlabel('k (zero-based)', fontsize='20')
plt.ylabel('Matrix Profile Value', fontsize='20')
plt.xticks(range(mps.shape[0]))
plt.plot(1, 1.3, marker="v", markersize=10, color='red')
plt.suptitle('Choosing k with discords={0}'.format(discords_md), size=20)
plt.show()

# Matrix profile subspace
k = 1
S = stumpy.subspace(df, m, discords_idx[k][0], indices[k][discords_idx[k][0]], k)

print(f"For k = {k}, the {k + 1}-dimensional subspace includes subsequences from {df.columns[S].values}")

[seanlaw]

@dbolotov Thank you for your question and welcome to the STUMPY community!

Is it necessary to set discords=True in stumpy.mstump() when looking for discords?

Yes, it is necessary to set discords=True in stumpy.mstump() when you are looking for discords.

From the doc for the discords parameter: "When set to True, this reverses the distance matrix which results in a multi-dimensional matrix profile that favors larger matrix profile values (i.e., discords) rather than smaller values (i.e., motifs).". What exactly is meant by reversing the distance matrix - is it flipped somehow?

Yes, exactly! I admit that it may quite dense but I encourage you to (re-)read this section of the mstump tutorial (and feel free to ask questions).

So, for a given multi-dimensional subsequence, we compute a corresponding multi-dimensional distance profile, D. D contains the distance from this multi-dimensional subsequence to every other subsequence within the same dimension (i.e., we only compute distances between T1 and T1 or T2 and T2 but never between T1 and T2). However, before we compute the multi-dimensional matrix profile using D, we first sort D by column:

https://github.com/TDAmeritrade/stumpy/blob/ed3abb5d082565815965a9f5afbe35dd08ff774d/stumpy/mstump.py#L816-L819

In the case of motifs, D is sorted in ascending order (smallest distances at the top and largest distances at the bottom). In the case of discords, D is sorted in descending order (largest distances at the top and the smallest distances at the bottom). Then, with D sorted, we compute the multi-dimensional matrix profile for this specific multi-dimensional subsequence:

https://github.com/TDAmeritrade/stumpy/blob/ed3abb5d082565815965a9f5afbe35dd08ff774d/stumpy/mstump.py#L680-L687

Thus, responding to your original question, yes, the multi-dimensional distance profile is reversed when discord=True and this is done so that we can force mstump to prioritize/prefer a larger (average) distances (i.e., discords) over smaller distances (i.e. motifs).

If discords=True, can the resulting MP still be used to find motifs?

Great question and the answer is a NO! Unfortunately, in order to get motifs, you must set discords=False.

Including code example below (reproduced from the documentation with changes to find discords). I first set discords=False in stumpy.stump(), and calculated the 3 top discord indices using discords_idx = np.argsort(-mps, axis=1)[:,:3].

This sounds right and I believe that you've identified a possible error in the tutorial. We should indeed be sorting using -mps as you've correctly identified. I've corrected that now, thank you!

Then I plugged discords_idx into where motifs_idx is used when finding the best k with the elbow method and the subspace. Then repeated the whole thing with discords=True. Also including two sets of plots - one where discords=False, and one where discords=True. The green vertical lines represent the top 3 discords for each k, and green segments show the area corresponding to the top discord for k=1.

Given the points that I made above (and please follow up with questions if anything is unclear), you cannot use the multi-dimensional matrix profile computed with discords=True to find motifs. Once you've set discords=True then you must set discords=True everywhere. Otherwise, the results are meaningless. For completeness, to find motifs, you must (re)-compute the corresponding multi-dimensional matrix profile with discords=False everywhere. Hopefully that makes sense.

Additionally, in the case of motifs, we may use the "elbow" method to choose k. Of course, for discords, this is reversed/flipped and so it'd be more appropriate to use the "knee" method to choose k.

[dbolotov]

Thank you for the explanation, things are a lot more clear now.

Reading through this tutorial section, I wasn't sure about two small points:

In our example, T has exactly l=5 multi-dimensional subsequences (i.e., we can incrementally slide a 4×3-shaped rectangle across the length of T only 5 times before we reach the end of T) and, for the ith multi-dimensional subsequence, we can iterate over each of its dimensions independently and compute an aggregated “multi-dimensional distance profile” (i.e., three 1-dimensional distance profiles stacked one on top of each other). Essentially, the ith multi-dimensional distance profile has shape d×l (or 3×5) and gives you the pairwise distances between the ith multi-dimensional subsequence and all other possible multi-dimensional subsequences within T.

• should "three 1-dimensional distance profiles" be "four 1-dimensional distance profiles"?
• should "d×l (or 3×5)" be "d×l (or 4×5)"?

A couple a more questions about mstump (please let me know if these should be in a separate thread):

• Are there any plans for streaming version, like there is stumpi for the univariate case?
• Your comment here mentions turning off z-normalization in stumpy.stump when looking for amplitude-based anomalies. Does this apply to mstump, and would there be any reason to not turn off z-normalization for discords in either univariate or multivariate case?

Thanks again.

[seanlaw]

should "three 1-dimensional distance profiles" be "four 1-dimensional distance profiles"?
should "d×l (or 3×5)" be "d×l (or 4×5)"?

Yes, you are absolutely right on both points! Originally, when the example was first written, it was a 3-dimensional example (i.e., T1, T2, T3) and then we realized that, since m = 3, it was less ambiguous/confusing to make it 4-dimensional (i.e., T1, T2, T3, T4). Of course, in doing so, we failed update the two points accordingly. I've just updated the tutorial but please let me know if you come across anything else. Your questions have been just as valuable and I'm sure that other users will appreciate your keen eye to detail as well so thank you!

Are there any plans for streaming version, like there is stumpi for the univariate case?

I think the short answer is "no" but not because we aren't interested but because, if I recall correctly, it has only been brought up once on Twitter (was that you?). mstumpi would be interesting but, considering how complex mstump was to implement, mstumpi might be challenging as well. It's certainly doable but one would need to assess the value and consider things like:

1. As the number of dimensions increase, will mstumpi be too slow for streaming? Will it be a problem when egress=False?
2. Currently, mstump returns the global multi-dimensional matrix profile (i.e., a nearest neighbor can be to the "left" or to the "right" of a given subsequence). However, for streaming data, we don't get to see the future and so we are almost always interested in the left multi-dimensional matrix profile (i.e., left nearest neighbor). This is currently unavailable in mstump and we'd need to also add this feature to the subspace function as well. This is not insurmountable but there are a good number of moving parts.
3. The choice of k is still somewhat of a human-in-the-loop process even for batch mstump. How could/would we update/automate this for streaming data? The original paper mentions a different approach called MDL for choosing k but we could not get it to work.
4. Overall, STUMPY tries to take a lot of care in what it delivers to users so to avoid having users accidentally "shoot themselves in the foot". Ultimately, users come back to us for questions/support and so we really want to have at least thought through most of those questions first. Additionally, with limited development time and resources, we'd want to prioritize accordingly so as to create a sustainable open source environment. Help is always welcome!

These are just a few thoughts/considerations off the top of my head. Do you happen to have an application for mstumpi? Would you be interested in helping to contribute this feature (I can help support you)?

Your comment here mentions turning off z-normalization in stumpy.stump when looking for amplitude-based anomalies. Does this apply to mstump, and would there be any reason to not turn off z-normalization for discords in either univariate or multivariate case?

Great question. By definition, anomalies are just hard to find. From my (limited) experience, if you are looking for spikes in your time series then setting normalization=False tends to be "better" suited for this task. Fundamentally, when you z-normalize a subsequence, you are basically saying "I don't care that my time series is shifted up or down. I only care about comparing the overall shape of one subsequence to the shape of all of the other subsequences." So, if you care about comparing subsequence shapes (even if the shape is shifted upward, say, due to seasonal trends) then normalization is your friend. However, if you care about changes in amplitude then normalization is your enemy. This is one of those cases where you kind of need to test-and-learn as the need for normalization is data-specific and there is no universal guidance that will cover all data scenarios.

[dbolotov]

Thanks for the answers. Our application is anomaly detection in the IoT space, working with streaming multivariate sensor data, but not bound to any specific type of measurements. We were just researching currently available methods, and looking for something user-friendly and with few parameters. We started by looking at algorithms like robust random cut forest, and eventually found out about matrix profiles & STUMPY (it was not me asking about it on Twitter :)).

Thanks for the list of challenges with a streaming version. I will have some internal discussions with my team about how much time we can spend contributing to implementation of mstumpi, we are currently not sure.

[seanlaw]

Sounds good!