MSTUMP Tutorial Clarification

[lukail]

Hi,

I was looking through the tutorial for mstump and either my understanding is incorrect or I believe there may be an error. Specifically, under Multidimensional Matrix Profiles, it is mentioned that,

Consequently, the “k-dimensional motif” can be found by locating the two lowest values in the correspond -dimensional matrix profile, mps, (these two lowest values must be a tie).

Why are we choosing the two lowest values here? Is this just arbitrarily selecting the two most probable motifs? From the code below this, it seems like we are selecting along index 1 (row 2), meaning the second dimension right? And also, why must these two lowest values be a tie?

An additional question: When choosing k for when we're interested in discords instead, is it as simple as setting discords=True when running stumpy.mstump and plotting the maximum matrix profile value in each dimension against k, then taking the max curvature?

[seanlaw]

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

Why are we choosing the two lowest values here? And also, why must these two lowest values be a tie?

So, this goes back to the concept of a "global motif" (or global minima) in that the global minima must be a pair. In the case of a 1D time series, the two subsequences that are identified at the global minima MUST be the one-nearest neighbor of each other (i.e., for a global minima, if subsequence A is the nearest neighbor to subsequence B then subsequence B must also be the nearest neighbor to subsequence A). This is why the two lowest values must be a tie. For a global minima, no other subsequences should be closer in distance (otherwise it wouldn't be called a global motif). Thus, by choosing the two lowest values, we are essentially choosing the index location of the two subsequences represented by the global minima.

Hopefully that makes sense but please feel free to ask for any clarification.

You now have me thinking whether this "tying pair" also holds for more than one dimension?! Hmmmm. I will have to think about that and I may need to edit that part out.

From the code below this, it seems like we are selecting along index 1 (row 2), meaning the second dimension right?

Are you referring to this code motifs_idx = np.argsort(mps, axis=1)[:, :2]? If so, let's break this down:

The mps should show:

array([[3.01630583, 2.83720883, 2.68261467, ..., 2.96434895, 3.18519577,
        3.49061656],
       [3.74779845, 3.58822096, 3.49334978, ..., 4.50403256, 4.46646621,
        4.55796841],
       [4.89094327, 4.61814682, 4.55491406, ..., 5.80773484, 5.48623554,
        5.37911065]])

And then np.argsort(mps, axis=1) means to sort the values along each dimension (the axis=1 corresponds to sorting within each row independently from the other rows) from smallest to largest but, instead of returning the values, return the index location:

array([[352, 152, 351, ..., 172,  73, 392],
       [351, 151, 352, ...,  73, 391, 516],
       [365, 108, 107, ..., 518, 517,  74]])

And then the additional indexing after sorting [:, :2] simply grabs all rows but only the first two columns from each row:

array([[352, 152],
       [351, 151],
       [365, 108]])

So, in the plot proceeding this code, we plot the profiles for P1, P2, and P3 since, technically, we don't yet know how many dimensions are "useful". The red arrows (in P1 and P2) indicate what we'd expect to be consistent with the "ground truth" but if we blindly chose P3 (i.e., to only allow patterns from all possible dimensions) then we'd identify a pattern that wasn't really "there".

Let me know if that helps.

When choosing k for when we're interested in discords instead, is it as simple as setting discords=True when running stumpy.mstump and plotting the maximum matrix profile value in each dimension against k, then taking the max curvature?

Well, in addition to setting discords=True in stumpy.mstump there are probably a few additional steps:

1. Instead of using np.argsort(mps, axis=1)[:, :2], which finds the index location of the smallest distances you'd want to use np.argsort(-mps, axis=1)[:, :2] in order to find the index location of the largest distances. Or alternatively, np.argsort(-mps, axis=1)[:, :2] is equivalent to np.argsort(mps, axis=1)[:, -2:] (though the order is different so the first is preferred).
2. When computing the subspace for a give dimension, you'll also want to set discords=True as in S = stumpy.subspace(df, m, discords_idx[k][0], indices[k][discords_idx[k][0]], k, include=[0, 1], discords=True)

This section of the tutorial provides some guidance but let me know if that is helpful!

[lukail]

Thank you so much for the welcome and the detailed answer Sean! This clarified a lot about my understanding, especially the pair and why they should be tie. Now, in the case of discords, am I correct in believing that discords are not a pair? Also, thanks for mentioning that I need to set discords=True when calling subspace. Missed that one!

[seanlaw]

Thank you so much for the welcome and the detailed answer Sean! This clarified a lot about my understanding, especially the pair and why they should be tie.

It was a great conceptual question! Glad that I could help clarify it.

Now, in the case of discords, am I correct in believing that discords are not a pair? Also, thanks for mentioning that I need to set discords=True when calling subspace. Missed that one!

Correct, the top discord (global maxima) is not guaranteed to be a pair. In fact, all nearest neighbors that are not the global minima has no guarantee to be a pair. A "pair" only corresponds to the global motif (i.e., global minima) and this is an exclusive property

[jackcan2]

Well, in addition to setting discords=True in stumpy.mstump there are probably a few additional steps:

1. Instead of using np.argsort(mps, axis=1)[:, :2], which finds the index location of the smallest distances you'd want to use np.argsort(-mps, axis=1)[:, :2] in order to find the index location of the largest distances. Or alternatively, np.argsort(-mps, axis=1)[:, :2] is equivalent to np.argsort(mps, axis=1)[:, -2:] (though the order is different so the first is preferred).
2. When computing the subspace for a give dimension, you'll also want to set discords=True as in S = stumpy.subspace(df, m, discords_idx[k][0], indices[k][discords_idx[k][0]], k, include=[0, 1], discords=True)

Hey @seanlaw - does stumpy.mdl have the same interpretation with discords? I am working on a multivariate anomaly detection problem and would like to know which dimensions contribute most to the anomaly. Taking the following code snippet from the tutorial, and modifying to use mdl, would we still take the k where the number of bits is smallest in mdls?

mps, indices = stumpy.mstump(df, m, discords=True)
discords_idx = np.argmax(mps, axis=1)
nn_idx = indices[np.arange(len(discords_idx)), discords_idx)
k = 2
mdls, subspaces = stumpy.mdl(df, m, discords_idx[k], nn_idx[k], k, discords=True)

I ask because I understand MDL to select the k with the optimal compression for two similar subsequences (in the motif case). However, I would expect the nearest neighbor for a discord to be dissimilar. As a result, wouldn't the compression be bad in every case? If so, are we then looking to take the k with the "least bad" compression in this case? Will this answer the question of which dimensions contribute most to the anomaly? Hope that makes sense - thanks!

[jackcan2]

To clarify my question, here is an example with a 3D time series and a discord injected at t = 100:

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

# Create time series
rng = np.random.default_rng(44)
ts1 = np.concatenate([np.repeat(0.25, 100) + rng.normal(scale=0.05, size=100),
                      np.repeat(0.85, 15) + rng.normal(scale=0.05, size=15),
                      np.repeat(0.25, 100) + rng.normal(scale=0.05, size=100)])
ts2 = np.concatenate([np.repeat(0.85, 100) + rng.normal(scale=0.05, size=100),
                      np.repeat(0.5, 15) + rng.normal(scale=0.05, size=15),
                      np.repeat(0.85, 100) + rng.normal(scale=0.05, size=100)])
ts3 = np.repeat(0.5, 215) + rng.normal(scale=0.05, size=215)
df = pd.DataFrame({'T1': ts1, 'T2': ts2, 'T3': ts3})

# Plot
fig = plt.figure(figsize=(15, 7))
ax = fig.gca()
for col in df.columns:
    ax.plot(df[col], alpha=0.6, label=col)
plt.legend()
plt.show()

Here you can see the discord at t = 100, where T1 and T2 are the sole contributors. If I run mstump and use a modified version of your plotting code, I get the following:

# Run mstump
m = 20
max_k = 3
mps, indices = stumpy.mstump(df, m, discords=True, normalize=False)
discords_idx = np.argmax(mps, axis=1)
nn_idx = indices[np.arange(len(discords_idx)), discords_idx]
mdls, subspaces = stumpy.mdl(df, m, discords_idx, nn_idx, discords=True, normalize=False)

# Plot
fig, axs = plt.subplots(mps.shape[0] * 2, sharex=True, gridspec_kw={'hspace': 0}, figsize=(10, 7))
for k, dim_name in enumerate(df.columns):
    axs[k].set_ylabel(dim_name, fontsize='20')
    axs[k].plot(df[dim_name].values)
    axs[k].set_xlabel('Time', fontsize ='20')
    axs[k + mps.shape[0]].set_ylabel(f"P{k}", 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=discords_idx[max_k - 1], linestyle="dashed", c='black')
    axs[k + mps.shape[0]].axvline(x=discords_idx[max_k - 1], linestyle="dashed", c='black')

    if k < max_k:
        axs[k].plot(range(discords_idx[max_k - 1], discords_idx[max_k - 1] + m), 
                    df[dim_name].iloc[discords_idx[max_k - 1] : discords_idx[max_k - 1] + m], 
                    c='red', linewidth=4)
        axs[k + mps.shape[0]].plot(discords_idx[k], mps[k, discords_idx[k]] + 1, marker="v", markersize=10, color='red')
    else:
        axs[k + mps.shape[0]].plot(discords_idx[k], mps[k, discords_idx[k]] + 1, marker="v", markersize=10, color='black')
plt.show()

As expected, the discord is easily detected where it should be. However, when inspecting the results of mstump.mdls, we see this:

print(mdls)
print()
print(subspaces)

This seems to imply that the optimal value of k should be k = 0, but I know that k = 1 here. Is my logic sound here? This is awesome at finding the anomaly - now I just need to figure out if it's possible to identify the dimensions driving any given anomaly. Thanks again!!!

[seanlaw]

@jackcan2 These are good questions and I'm afraid I don't have the answer. Multi-dimensional time series are quite hard to interpret and even I'm unable to easily think through whether or not the discord bit makes sense. I'm wondering if @SaVoAMP has any thoughts on this as she has spent a good amount of time thinking about multi-dimensional space?

[SaVoAMP]

Taking the following code snippet from the tutorial, and modifying to use mdl, would we still take the k where the number of bits is smallest in mdls?

I think using the mdl approach to find the anomalies doesn't make sense here. At least not when trying to compare a discord with its nearest neighbor. Mdl states that the model in which the observed data can be compressed the most is likely to be the "true" model. This means - in the sense applied to the dimensionality selection problem - that the value for the number of relevant dimensions can be chosen, whose k-dimensional motif requires the smallest number of bits to store. More precisely the number of bits needed to store each of the k-dimensional motifs is calculated. This is accomplished by compressing the multidimensional subsequence pair that spans the motif subspace proposed by the k-dimensional matrix profile. Purely intuitively, one can imagine that to describe two multidimensional subsequences that are visually very different from each other, one would need a lot of memory to highlight all the differences whereas one would to easily and not need much memory if one wanted to describe two particularly similar subsequences. Since an anomaly subsequence should be very different from its nearest neighbor, I think that it wouldn't make any sense to compare these two subsequences with mdl, because mdl only tells us in how many dimensions one subsequence is very similar to the other. But we don't want to find those dimensions in which strongly differing subsequences are the most similar, just as little as the number of dimensions in which the differing subsequences are most similar.

I haven't looked closely at the MDL approach yet, unfortunately, but I wonder if you don't need a completely different, complementary mathematical approach here. However, it may also be possible to modify the MDL approach so that you get some sort of "Maximum Description Length" (instead of Minimum). I only know that for the MDL frame the storage space requirement for lossless compression of the time series is relevant. To do this, one time series is used as a hypothesis time series to encode another. The encoding is then done by storing only the positions and values of the differing time series elements (this is why it is also referred to as a differential encoding method). Here a reduced description length function is used, which is necessary for the storage of a discretized time series, which is compressed relative to another time series:

T_H is the hypothesis time series here, m the length of the subsequence, b the number of bits and γ() is a function which counts the number of non-zero entries of the input vector and therefore the number of non-matching values of the time series T and T_H. The memory required to store the positions of non-matching values is contained in the term log_2m as well as in b the memory for their concrete values. I'm not sure but I guess it should be possible to use the RDL and try to find a Maximum Description Length. Whether you then get your desired result is probably still questionable. One would have to test that perhaps on a few simple constructed examples 🤔

[seanlaw]

Thanks for sharing your thoughts @SaVoAMP! To be 100% clear though, the multi-dimensional matrix profile is still "correct" and the largest values (in each dimension) still corresponds to the index location of the discords. Additionally, the mdl and subspsace values are also correct by definition and, with the discords=True parameter, the subspace is actually sorted from largest to smallest (instead of smallest to largest as is the case for motifs) before the top-k dimensions with the largest distances are used for the actual mdl calculation. What is important is that the subspace derived from the pair of multi-dimensional subsequence is NOT based on the mdl and is simply picked from the top-k matrix profile dimensions with the largest values. So, if mdl (compression) isn't right for discords then you can still make use of the subspace and figure out the which k dimensions is important. In fact, given that the kth subspace corresponds to the k dimensions with the largest matrix profile distance (between nearest neighbors), I'm wondering if it is sufficient to simply take the k with the LARGEST mdl value rather than the smallest mdl value? Ultimately, I guess it depends on your goal.

[jackcan2]

@seanlaw @SaVoAMP Thanks a ton for your insightful responses! I'll keep chugging away over here and follow up if I have anything new to report.