How To Analyze Multi-dimensional Time Series with mstump

[Darveesh]

If my time series has more than 1 value, say a geolocation - i.e. lat and long - over time, can you provide an example or how to use STUMPY to analyze the the time series signal?

[seanlaw]

It depends on what you are trying to accomplish. The first thing you could try is to compute the matrix profile for each variable independently from the other using stumpy.stump. Chances are, if no motifs exists in a single dimension then they are not likely to exist in multiple dimensions.

Alternatively, one could compute what is called a "multi-dimensional matrix profile" using stumpy.mstump but we do not currently have a tutorial for this (it is still a work in progress and a lot more complex than you'd think and so support may be limited). We recommend going through this paper before venturing down this path as it can be rather nuanced.

What have you tried already? What are you trying to accomplish?

[seanlaw]

Additionally, if you haven't had a chance to look through then, you may find some useful ideas in our tutorials.

[Darveesh]

An example for what I am trying to accomplish is analyzing a bus route. A bus takes the same route daily so there is periodicity to the data (lat/lon of the bus tracked over time every day). I'd like to identify when the bus "deviated" from its route. This could be getting stuck (say in traffic) or perhaps even taking an alternate route (due to detour). The data I have is a time series data of lat/lon over a whole year. Hope that helps explain what I am trying to accomplish.

In this vein, maybe each variable could be treated independently? I would think a delay does affect both variables, same with taking a detour. So ideally perhaps this is multi-variate analysis and really suited for stump.mstump, but I can fake it with single variable analysis due to the dependency between the two variables.

Eventually though I suspect I need stumpy.mstump. Only because I imagine I will introduce other (related) signals, for example the causes for the diversion. So if you get to polishing up stumpy.mstump you got a customer here. :)

Thank you again for the guidance.

[seanlaw]

An example for what I am trying to accomplish is analyzing a bus route. A bus takes the same route daily so there is periodicity to the data (lat/lon of the bus tracked over time every day). I'd like to identify when the bus "deviated" from its route. This could be getting stuck (say in traffic) or perhaps even taking an alternate route (due to detour). The data I have is a time series data of lat/lon over a whole year. Hope that helps explain what I am trying to accomplish.
In this vein, maybe each variable could be treated independently? I would think a delay does affect both variables, same with taking a detour. So ideally perhaps this is multi-variate analysis and really suited for stump.mstump, but I can fake it with single variable analysis due to the dependency between the two variables.
Eventually though I suspect I need stumpy.mstump. Only because I imagine I will introduce other (related) signals, for example the causes for the diversion. So if you get to polishing up stumpy.mstump you got a customer here. :)
Thank you again for the guidance.

This is a pretty interesting example! So, technically, stumpy.mstump does compute the multi-dimensional matrix profile correctly (note that it isn't simply stacking individual 1-d matrix profiles on top of each other!!). The "work-in-progress" is referring to some of the tools that are needed for figuring out which dimensions are "important". That is a bit trickier.

In your case, you feed in a 2D time series (lat/lon) to stumpy.mstump and then compute the multi-dimensional matrix profile and assess the profile in the second dimension (i.e., it must consider both lat and lon). Then, if they take a detour then you might get a "high" value in your matrix profile. What's important to keep in mind is that STUMPY is comparing all "trips" with all other trips. So, if your bus had taken the same detour in the past then any upcoming detours that follow the same (past) detour path will "match" and therefore result in a small matrix profile value. So, if you want to detect deviations then you'd want to first start with a clean dataset where the bus NEVER took a detour and use that as your "normal" (one week worth of normal data should be enough). Then, you might append (or concatenate) new data to the normal data and assess from the matrix profile for the "new" part of the data has low values. If the newly appended data is a "never-seen-before" in the normal data then it will result in a higher matrix profile value. Then, each day, you remove yesterday's "new" data and replace it with today's "new" data (while never changing your "normal" data).

I hope that helps.

[seanlaw]

Btw, is this bus data publicly available? If so, I'd like to take a look at it to see if it could serve as an example in our tutorial.

[Darveesh]

Ok I'll take a look at mstump. Is there a tutorial for out I can use as guidance?

Ps: The bus data I have access to seems to come from a few sources and not a single public source. If I can share any of my data will let you know.

[seanlaw]

There is this tutorial that is incomplete and still evolving and so it might generate more questions than answers. I still recommend reading the original paper mentioned above.

[Darveesh]

I read through the paper and the in-progress tutorial, thank you. I am about to embark on feeding my data to mstump. I have isolated a single bus's (1512) location for about a ten day window (attached). My next question is around time spacing. The data I have at the moment is timestamped location data. For MP / stumpy to consume this data it needs to be at equal time intervals I suppose? If so, would I need to pre-fill/pad the data in such a way that all data points are exactly the same "distance" in time? Or can stumpy somehow process date/time stamps that are not evenly distributed?

bus 1512 time series.xlsx

[seanlaw]

I have isolated a single bus's (1512) location for about a ten day window (attached). My next question is around time spacing. The data I have at the moment is timestamped location data. For MP / stumpy to consume this data it needs to be at equal time intervals I suppose? If so, would I need to pre-fill/pad the data in such a way that all data points are exactly the same "distance" in time? Or can stumpy somehow process date/time stamps that are not evenly distributed?

Technically, STUMPY (or matrix profiles in general) do not actually take the "time" into consideration when computing a matrix profile. Instead, it processes and ordered sequence or series of data. At least, that's my mental model. You are correct in that STUMPY assumes that your data points should be (roughly) equally spaced and we do not provide any tools for that pre-processing work as there are many other specialized tools that can help with that (i.e., Pandas and NumPy are great). When you have missing data points, you may consider setting those values to np.nan or using other data imputation techniques like averaging or filling forward (with the neighboring known data points).

[Darveesh]

@seanlaw I saw there is a flag for discords. I am running stumpy.mstump as such (with my data in 1 minute increments):

m=30
mp, indices = stumpy.mstump(df[["lat", "lon"]], m, discords=True)

In the discord case, how do I access the index of the highest (absolute or 2 std above. mean) MP value? A bit confused on the all the dimensions. :)

[Darveesh]

To expand a bit on my confusion, I went back to motifs and I added the following debug lines to try to understand the lows and highs in the MP, but I think I am missing the idea of a multi-dimensional MP. In the tutorial, when each time signal is analyzed individually, the motif indices are close but not identical. Later when the multi-dimensional MP analysis is done, only the motif indices from the first time signal are used (and I presume assumed to be equal to the 2nd signal). On my dataset, when running the code below:

logger.info("Starting stumpy...")
m=30
mp, indices = stumpy.mstump(df[["lat", "lon"]], m, discords=False)
logger.info("Stumped!")

logger.debug(f"MP shape: {mp.shape}")
logger.debug(f"Motif Indecis lat: {np.argpartition(mp[:, 0], 2)[:2]}")
logger.debug(f"Motif Indecis lon: {np.argpartition(mp[:, 1], 2)[:2]}")
logger.debug(f"Max index lat: {np.argwhere(mp[:, 0] == mp[:, 0].max()).flatten()[0]}")
logger.debug(f"Max index lon: {np.argwhere(mp[:, 1] == mp[:, 1].max()).flatten()[0]}")

I see this output:

12/28/2020 04:23:43 PM Starting stumpy...
12/28/2020 04:23:58 PM Stumped!
12/28/2020 04:23:58 PM MP shape: (13573, 2)
12/28/2020 04:23:58 PM Motif Indecis lat: [8232 8186]
12/28/2020 04:23:58 PM Motif Indecis lon: [5326 4159]
12/28/2020 04:23:58 PM Max index lat: 12965
12/28/2020 04:23:58 PM Max index lon: 7571

So the np.argpartition(mp[:, 0], 2)[:2] approach - which I think is trying to find the index of the 2 lowest MP values - is not yielding the same index locations on signal 1 vs signal 2 (for me). Further, if I calculate the max index (for discord detection), I again get different locations. So if the indices are supposed to be the same how come they aren't? And if they can be different, then how is this different from doing 2 1-dimensional stumpies as in the tutorial?

Thanks.

[seanlaw]

@Darveesh There is a lot to unpack here and I'm afraid that trying to follow a tutorial that is still a "work-in-progress" may not be the best idea as the code and explanations are still continually being improved and refined. If there is no rush then I recommend that you wait until the tutorial is completed (I am working on it right now and hope to have it published by next week). As I had mentioned before (and as you correctly pointed out), the multi-dimensional matrix profile is indeed not a simple a stacking of two (or more) 1-dimensional matrix profiles. When you read the Matrix Profile VI paper in detail, "Algorithm I, lines 10-15" are the key. Essentially, what mstump is doing is the following:

1. Let's say that we have a three dimensional time series (just like in the tutorial example) and they all have the same length
2. For the ith subsequence (with length m) in each time series, we compute all pairwise distances between this subsequence and all other subsequences within the same time series. So, we compute all of the distances between T1[i : i + m] against all other subsequences in T1. Then we repeat this with T2[i : i + m]and, separately, with T3[i : i + m]. Note that T1 subsequences are only compared with T1 and never with T2 or T3. These pairwise distance vectors (for the ith subsequence) are what is referred to as a "distance profile" and let's refer to them as D1, D2, and D3, respectively. Technically, D1, D2, and D3 are the same length and, conceptually, we stack them one on top of each other and call this D123. So, D123 will contain three rows (i.e., one row for each dimensions). This will be relevant later. For example, let's say that D123 looks like this:

D123 = [[1.1, 0.2, 3.4, 0.4],
        [0.3, 1.8, 1.2, 3.0],
        [2.5, 0.1, 2.2, 1.6]]

Note that each time you slide i across the time series to another subsequence location, you will create a brand new D123!
3. Next, we perform a column-wise sort (see "Algorithm I, line 10") on D123 (i.e., the values within a column are sorted from smallest to largest). And, yes, this may seem silly at first since this sorting seems to cause us to lose track of which distance corresponds to which time series but know that it is okay and we'll recover this information when we deal with "subspaces". Nonetheless, D123 still contains three rows/dimensions except that they are sorted by columns now:

D123_sorted = [[0.3, 0.1, 1.2, 0.4],
               [1.1, 0.2, 2.2, 1.6],
               [2.5, 1.8, 3.4, 3.0]]

4. Now, we can compute the "matrix profile values" for the ith subsequence window (this is what eventually gets output by mstump). To fill in the first dimension of the matrix profile (i.e., the first column or mp[i, 0]) we are allowed to choose one time series (or dimension) to "use" out of the three. In this case, we choose the smallest value along the first row/dimension of D123_sorted and use that value and index position as our "matrix profile value" and "matrix profile index" for the first dimension. Since we are using D123_sorted where information was lost, what we're basically doing here is saying, "I don't care which time series out of the three, just tell me the value and location of the smallest distance between the ith subsequence and its nearest neighbor within a single dimension!" Notice that we don't actually need to know exactly "which" time series this "value and location" are derived from. In this case, the matrix profile value is 0.1 and the index location is 1.
5. Next, we can fill in the second dimension of the matrix profile (i.e., the second column or mp[i, 1]) but this time we get to choose two time series (or dimensions) to "use" out of the three. In this case, we choose the smallest average value along the first and second rows/dimensions of D123_sorted and use that average value and index position as our "matrix profile value" and "matrix profile index" for the second dimension. Here, what we're basically doing here is saying, "I don't care which two time series out of the three, just tell me the value and location of the average smallest distance between the ith subsequence and its nearest neighbor but it needs to come from two dimensions!" In this case, we take D123_sorted and average along the columns but only for the first two rows:

[0.7, 0.15, 1.7, 1.0]

Since the smallest average value is 0.15 then the matrix profile value for the second dimension is 0.15 and the index location is (coincidentally) 1.
6. Finally, we can fill in the third dimension of the matrix profile (i.e., the third column or mp[i, 2]) but this time we get to choose three time series (or dimensions) to "use" out of the three. In this case, we choose the smallest average value along the first, second, and third rows/dimensions of D123_sorted and use that average value and index position as our "matrix profile value" and "matrix profile index" for the third dimension. Here, what we're basically doing here is saying, "I don't care which three time series out of the three, just tell me the value and location of the average smallest distance between the ith subsequence and its nearest neighbor but it needs to come from all three dimensions!" In this case, we take D123_sorted and average along the columns but we include all three rows:

[1.3, 0.7, 2.26666667, 1.66666667]

Since the smallest average value is 0.7 then the matrix profile value for the second dimension is 0.7 and the index location is (coincidentally) 1. Note that as you move from dimension to dimension the location of the index (to its nearest neighbor) and move and need not stay the same.

Please read this over and feel free to ask any clarifying questions. It took me a long time to fully grasp all of the intricacies but I promise you that the logic described in the paper is sound.

[Darveesh]

Thank you for the detailed writeup @seanlaw. I think I am beginning the grasp the idea: At position i -

1. Compute the individual time series distance profiles.
2. Then pick your count of time series to consider (i.e. dimensions). Note: The returned MP has the full range of values from 1 to X time series count/"dimension".
3. Average the lowest X distance profiles and pick the smallest as your MP value for that X "dimension".

So now I get what you mean by the MP not being individual signal MPs stacked on top. The dimensions of the returned MP from stumpy.mstump correspond to how many signals the distance profiles are averaged/computed over. They are not the individual signal MPs.

So in my case, I have 2 times signals (lat and lon) and I want both dimensions to be considered. so mp[:,1] is what I am after. Beyond that, it's back to the basics of finding the minima and maxima of mp[:,1] like regular stumpy.stump.

Hope I got that somewhat correct.

[seanlaw]

So now I get what you mean by the MP not being individual signal MPs stacked on top. The dimensions of the returned MP from stumpy.mstump correspond to how many signals the distance profiles are averaged/computed over. They are not the individual signal MPs.

That's right! In the paper, they say k out of d dimensions and so what you are calling X is what they are calling k. The tricky part is that, for each dimension, you pick the k smallest distances along each column of the multi-dimensional distance profile and you average those distances to give you the kth dimension of the ith matrix profile distance (and index).

So in my case, I have 2 times signals (lat and lon) and I want both dimensions to be considered. so mp[:,1] is what I am after. Beyond that, it's back to the basics of finding the minima and maxima of mp[:,1] like regular stumpy.stump.
Hope I got that somewhat correct.

Yesss, I'm so happy (and impressed) that you got it!! I'm also super relieved that my explanation made sense to you but please don't hesitate to point out where it may have been confusing and I can try to edit it to be more clear. I'm sure other users would benefit from this discussion.

[seanlaw]

In your case, since you actually care about 1) including both dimensions and 2) looking for discords (i.e., both lat and lon have to deviate and not just one) then I think you want to do something like:

import numpy as np
import stumpy

T = np.random.rand(30000).reshape((2, 15000))
mp, idx = stumpy.mstump(T, m=50, discords=True)

discord_idx = np.argsort(mp[1, :], kind="mergesort")[::-1][0]  # Note that this is different from np.argpartition above
nearest_neighbor_idx = idx[1, discord_idx]

It's getting late here so hopefully I didn't misspeak in any of the above.

[seanlaw]

It's actually for work. We have lots of time series data and are looking into ways to detect anomalies and correlations among the signals. Came across Stumpy in my research and here we are. :)

That's awesome! I rarely ever get to hear about the use cases or learn which companies are leveraging so I'd love to know what company this is (but no pressure).

One more question for today if its ok - what does the nearest neighbor indices[discord_idx, 1] signify? The second highest MP value (i.e. the next anomaly / largest distance profile)?

No problem! Essentially, the discord_idx tells you the location of where the discord is. Notice that this is a single index value and applies to k dimensions. Whereas indices[1, discord_idx] tells you where the closest match (or nearest neighbor) to the discord is located and, again, it is a singular index value and applies to k dimensions. While this second index may seem like something that is more useful for motifs, for discords, you can visually compare the discord and see precisely how much it differs from its closest match.

[Darveesh]

@seanlaw In our example above, when the index for the discord is determined, is it possible to know the subsequence index against which this subsequence was the most different? I.e. what was the reference subsequence that resulted in such a large distance profile. Maybe I am asking for farthest neighbor in the lingo above, not sure.

[seanlaw]

Just so we are on the same page, a matrix profile always stores the distance between each subsequence and its nearest neighbor (not the farthest neighbor or most dissimilar neighbor). A discord is choosing the subsequence that has a nearest neighbor that has a distance that is the farthest away compared with all other subsequences and their nearest neighbors. In the multidimensional case, we use averaged distances to the nearest neighbor instead.

In our example above, when the index for the discord is determined, is it possible to know the subsequence index against which this subsequence was the most different? I.e. what was the reference subsequence that resulted in such a large distance profile. Maybe I am asking for farthest neighbor in the lingo above, not sure.

We can probably agree that "discord_index" refers to the index location for a discord subsequence with length m within our time series, T. So, the actual discord is: T[discord_idx : discord_idx + m]. When you call mstump, you also get the indices that tell you where the nearest neighbor for each subsequence is and that is T[indices[discord_idx] : indices[discord_idx] + m]. Notice that we are still saying "nearest neighbor" here. Discords (and its matrix profile) are still measuring how badly dissimilar your nearest neighbor is. However, based on your question, it sounds like once you find a discord (which is based on having a "poor" nearest neighbor distance), you then want to find out something different. Basically, you want to know, given T[discord_idx : discord_idx + m], what is the subsequence that looks the most different? Is that right?

Note that we've glossed over a detail here in that we are dealing in multiple dimensions. So, if you have 2-dimensions, then it is possible that while the discord_idx is the same for each dimension (i.e., T[0, discord_idx : discord_idx + m] and T[1, discord_idx : discord_idx + m]), the complimentary subsequences that looks "most different" to T[0, discord_idx : discord_idx + m] and T[1, discord_idx : discord_idx + m] will more than likely exist in different places along their respective time series. You might still be able to average the complimentary distances across the multiple dimensions but, essentially, what you'll want to do is compute the distance profile for each of T[0, discord_idx : discord_idx + m] and T[1, discord_idx : discord_idx + m] using a private function called mstump._multi_mass:

from stumpy.mstump import _multi_mass
from stumpy import core

T, M_T, Σ_T = core.preprocess(T, m)
D = _multi_mass(T[:, discord_idx : discord_idx + m],
                T,
                m,
                M_T,
                Σ_T,
                M_T[:, discord_idx],
                Σ_T[:, discord_idx],
)

Here, D is your multidimensional distance profile for the multidimensional subsequence T[0, discord_idx : discord_idx + m] and T[1, discord_idx : discord_idx + m]. Perhaps, you can then traverse each dimension (row) of D to find out which index gives you the farthest distance:

sorted_D = np.argsort(D, axis=1)[::-1]  # Sort the distances independently in each row from largest to smallest and return the indices
indices_with_largest_distance = sorted_D[0] 

Hopefully, that's what you're looking for?

Fair warning: Considering that you are venturing into using a Python private function, this may be beyond the scope of support that I can provide and you may need to dig into the code on your own.

[Darveesh]

Thank you @seanlaw for the detailed writeup. I am still trying to digest your response... so please forgive any misunderstandings. Going back to my bus route example, I was investigating the bus' motion at the flagged discord_idx / "time". My m = 30, so given the 1 minute sampling, I understood that the 30-min window at the discord_idx position is the most different. But the question in my head was compared to what? I imagine the "reference" is the bus' motion along the same path at other times? Or may be that's the wrong way to be looking at this. Not sure. Anyhow, I manually looked for times when then bus travelled the same route and tried to deduce how these differed from the flagged window. Interestingly enough, I found an example in my dataset where the bus travels the same path a few hours later but the bus's position in that 30-min window seemed delayed. So it seems MP had flagged a delayed (or early?) bus. Great. My instinct was can I find this other "reference" 30-min window of bus motion automatically? I need to read your response a few more times to understand it better, but hopefully this explanation gives more color to the context of my question.

[seanlaw]

My m = 30, so given the 1 minute sampling, I understood that the 30-min window at the discord_idx position is the most different.

If I may, I would like to be more particular about the wording here. Technically, the window at the discord_idx position is the window that has a nearest neighbor that is farther away than any other subsequence. Notice, again, that I say "nearest neighbor" because, relative to T[discord_idx : discord_idx + m], there are other subsequences that can look "more different" than the "nearest neighbor". In fact, by definition, since we are saying "nearest neighbor" then it must be true that all other subsequences must necessarily be "more different".

I imagine the "reference" is the bus' motion along the same path at other times?

Sorry, I'm not certain what you are referring to here as "reference"

So it seems MP had flagged a delayed (or early?) bus.

I don't think MP can be that deliberate in flagging a "delayed bus". It simply found a window that had a poor "best match" (i.e., a nearest neighbor that was the most far away compared with the distance of any other subsequence and their nearest neighbor). So, MP flagged "something" but it has no idea as to what it flagged aside from having a nearest neighbor that is far away.

My instinct was can I find this other "reference" 30-min window of bus motion automatically?

Unfortunately, I'm still having a hard time grasping what you mean by "reference" here as I can interpret it in many different ways.

[seanlaw]

ICYMI: The MSTUMP tutorial has been published

[Darveesh]

Very nice! Seems i just jumped the gun by a couple weeks and caused you a lot of headache. :)

[seanlaw]

No worries! Your questions forced me to formulate my thoughts more clearly and I happened to be knee deep in mstump already. It was, as the say, “the perfect storm”. Happy new year! 🎉

[Darveesh]

Awesome. Happy new year to you as well! 🍾 🎊

[geoHeil]

One question @Darveesh : why don`t you rasterize the geolocations such as https://github.com/uber/h3 H3 then you should be able to get rid of one dimension.

[Darveesh]

@geoHeil Thank you for the pointer. I get that this will reduce the need for multivariate analysis in my example above tracking a bus' movements in time. However, my broader intention is to eventually obtain correlation with other time series from among a bunch of signals (perhaps to gain insight into causation for the bus' deviations). Say another time series that tracks accidents or congestion on the bus' route. So eventually I still think I need multi-variate analysis and this was a good exercise for me to understand MP and mstump. Now whether that ultimate goal is feasible or not with MP is another story I suppose. As the paper shows including non-correlated signals in the analysis will produce random/garbage results. So jury is out on that objective, but if there are other ways to find correlations in time signals I am open to suggestions. Am a noob in all this so any pointers are much appreciated.