Develop single-dimensional matrix profile for multi-dimensional time series

[NimaSarajpoor]

Motivation

I have noticed that there are some interest in using matrix profile in multi-dimensional time series. We can think of it as:

1. Multi-Dimensional Matrix Profile of Multi-Dimensional time series
2. Single-Dimensional Matrix Profile of Multi-Dimensional time series

I believe the first one is tackled in Matrix Profile VI, and its tutorial is under progress in PR #557 .

In the paper Matrix Profile VI, we can see the following note:

There are a few things to notice:

• So, it seems generalizing single-dimensional matrix profile works good for multi-dimensional data, particularly if the dimension is low (like 2 or 3) and if there is no irrelevant data there.
• The paper claims that if we generalize the concept to single-dimensional matrix profile, we can find the same motifs that were discovered before when matrix profile used for each time series. Well, in my opinion, it is hard to believe that the statement is true in all data sets. I mean, maybe having single-dimensional matrix profile gives us new information. (note that we should NOT simply add their matrix profiles to get one-dimensional matrix profile)

So, I think it is worth it to have support for single-dimensional matrix profile for multi-dimensional time series data.

---

Challenge-(I)
(this is based on what I read about generalizing single-dimensional matrix profile to multi-dimensional time series data. I do not remember the source though. I think it was from Eamonn Keogh.)

One of the main challenges is how to combine m distances across m-dimensional time series data.

Example:
Let's say we have two-dimensional data: T1 (first dimension) and T2 (second dimension). And, let's focus on subsequences at index i and j. Therefore:

Si_1 = T1[i: i+m]; Si_2=T2[i: i+m]
Sj_1 = T1[j: j+m]; Sj_2=T2[j: j+m]

d1 = dist(Si_1, Sj_1)
d2 = dist(Si_2, Sj_2)

D = f(d1, d2) # D: total distance between i and j

But, what is that function f for combining the two distances d1 and d2 (and gives one single value)? Two common options are:

• f(d, d') = d + d'
• f(d, d') = $(d^{p} + d '^{p})^{1/p}$ (p-norm)

This two are equal only when p=1. Otherwise, they are different. Both seems reasonable approach. However, we can go with the second approach and provide a module for that in stumpy (see section "implementation" below)

Challenge-(II)
Another challenge that I remember I read from the same source was how to avoid the domination of one dimension? well, this is probably an issue in non-normalized version. However, like many machine learning problems, it is usually up to the user on how to normalize time series. Maybe they normalize T1 and T2 by their maximum. Or, they may standardize the WHOLE data T1 and the WHOLE data T2. So, I believe this shouldn't be our concern.

---

Implementation
I think we can easily do $(d^{p} + d '^{p})^{1/p}$ for p-norm non-normalized matrix profile. The challenge might be in using Pearson correlation in normalized version. However, there is a nice solution for that!

$D^2 = d ^ 2 + d' ^ 2 = 2m (1-\rho) + 2m(1-\rho^{\prime}) = 2m(2 - (\rho + \rho^{\prime})) = 2[2m(1 - avg(\rho, \rho^{\prime})]$
Note that the factor 2 can be eliminated because if we scale all pairwise distances by the same number, it does not change the result.

Therefore:
$D = \sqrt{2m ( 1 - P)}$, where P is average of pearsons. So, we can use all those rolling/running variance stuff and simply just take average of pearsons!

Cool!

[seanlaw]

@NimaSarajpoor There's a lot to unpack here but here are a couple of points off of the top of my head:

1. You may already be aware of this (so ignore where appropriate) but the multi-dimensional matrix profile already combines dimensions by averaging the distances across a subset of dimensions where the distances are smallest. The details can be found here. However, not that multi-dimensional matrix profiles are not simply summing up distances for the k smallest distances and, instead, it is an average.
2. If your proposal is to simply compute the 1-D matrix profiles for each dimension independently of the others then, IMHO, there's nothing that STUMPY needs to "support" as the computation of 1-D matrix profiles is already optimized and any "combining" can already be done by the user in a post-processing step (outside of STUMPY). An example of this post-processing step might be a nice tutorial?

Maybe you can elaborate on what is currently missing as I may not be understanding your point clearly?

[seanlaw]

Adding @SaVoAMP @mihailescum to this conversation as they have thought about this a lot and may have some comments to contribute!

[NimaSarajpoor]

There's a lot to unpack here

Agree. I just wanted to share my idea before I forget about it.

1. You may already be aware of this (so ignore where appropriate) but the multi-dimensional matrix profile already combines dimensions by averaging the distances across a subset of dimensions where the distances are smallest. The details can be found here. However, not that multi-dimensional matrix profiles are not simply summing up distances for the k smallest distances and, instead, it is an average.

Thanks for sharing the link! I took a look and there are some similarities for sure. As provided in here:

ith_matrix_profile = np.full(d, np.inf)
ith_indices = np.full(d, -1, dtype=np.int64)

for k in range(1, d + 1):
    smallest_k = np.partition(ith_distance_profile, k, axis=0)[:k]  # retrieves the smallest k values in each column
    averaged_smallest_k = smallest_k.mean(axis=0) # line(*)
    min_val = averaged_smallest_k.min() 
    if min_val < ith_matrix_profile[k - 1]:
        ith_matrix_profile[k - 1] = min_val
        ith_indices[k - 1] = averaged_smallest_k.argmin()

IF I UNDERSTAND CORRECTLY:
For instance, say d=2 (two dimensional data T=[T1, T2]), then, in the last iteration where k=d, the averaged_smallest_k = smallest_k.mean(axis=0) is distance profile considering all dimensions (here, both two dimensions together). It gets average of each column in 2D distance profile.

What I am proposing is to simply keep min_val (in min_val = averaged_smallest_k.min()) of the last iteration k=d. So, skip that if-block, and simply return min_val as the value for each index in 1D matrix profile. (So, no for-loop or partition.) something like this:

1D_dist_profile = np.mean(ith_2D_distance_profile, axis=0)
# exclude trivial
idx = np.argmin(1D_dist_profile)

Instead of .mean(), we can combine the values of 2D_distance_profiles (across axis=0) by using the formula provided in my previous comment. We can take advantage of relationship between total_D and Avg_of_Pearsons(across dimensions). In fact, we can enhance current implementation of stump and stumped, and gpu_stump for this! No new module is needed.

Maybe I should find some 2D or 3D data to test it out and see if it gives me new/interesting insight about the data! I believe it gives me something new! However, I would like to test it in some real-world data @seanlaw @SaVoAMP @mihailescum do you have any suggestion for 2D or 3D data?

---

2. If your proposal is to simply compute the 1-D matrix profiles for each dimension independently of the others then,

That is not my proposal. I should have been more clear. Let's say D1 and D2 are two distance profiles of S1 and S2 at an index i. Then, min( combination of D1 and D2) is different than min(D1) + min(D2). The latter is simply the same as adding matrix profile of each dimension in post processing step. My idea is to calculate combined distance profile (similar to line(*) in code snippet above) and then return the minimum of that distance profile as the matrix profile value.

Note that we can use this for query matching as well! So, if I have 2D query Q=[Q1, Q2], I can find S=[S1, S2] close to Q considering both dimension. Again, this is not adding matrix profile of 1 and 2, but calculating 1D matrix profile by combining distance profile throughout the process. Is this useful? probably. The idea is that I now have a tool to explore and see if I can get new insight from data.

---

Is this reasonable when we have huge number of dimensions? Probably not. Because, as mentioned in the Eamon's paper, there might be some noise in some dimension of data.

---

Is it better than multi-dimensional matrix profile? I do not want to use the term "better". It may have its own advantage. I need to use it to see what kind of insight I can get from data that is different than multi-dimensional matrix profile. Also, I believe it is easier to understand. So, if I have 2 or 3 dimensional data, using this 1D matrix profile might be better in terms of interpretability. And, I think I can do matching in 2D! And, it can be simply(hopefully) implemented in the existing modules.

@seanlaw
Please feel free to close this issue or move it to discussion if necessary.

[SaVoAMP]

Maybe I should find some 2D or 3D data to test it out and see if it gives me new information! I believe it gives me something new! However, I would like to test it in some real-world data @seanlaw @SaVoAMP @mihailescum do you have any suggestion for 2D or 3D data?

I was working with a three-dimensional boxing data set (consisting of acceleration data of 8 different boxers) that is also labeled. I have found here that relatively similar results emerge when examining the data in one, two, or all three dimensions for a punch motif. I could obtain nice results, especially when concatenating all the punches of a boxer of the same type (for example, only frontal punches with the left hand), so that I only had to look for a single motif.

However, even with a 30-dimensional data set for analyzing human motion I could obtain relatively good results with different choices of k. Maybe this data sets will help to test your considerations.

[NimaSarajpoor]

@SaVoAMP
Thanks for sharing the data! Much appreciated! I will take a look and see if I can implement my idea and get some preliminary result. I will seek your help in the nearby future if I don't understand the data or the target (i.e, the problem we try to solve by using that data).

[seanlaw]

@NimaSarajpoor I think this is probably better suited for the Discussion section as it isn't quite an issue with the existing code and, instead, is a atypical and yet-to-be-confirmed case.

[NimaSarajpoor]

@seanlaw
Right! It is better to be put in discussion.

[NimaSarajpoor]

Toy Data
Let's say we have 3D time series data. In below, we show the data for each dimension.

>>> T1
array([10.1       ,  0.2       ,  0.3       ,  0.54488318,  0.4236548 ,
        0.64589411, 10.13      ,  0.24      ,  0.36      ,  0.38344152,
        0.79172504,  0.52889492,  0.56804456,  0.92559664,  0.07103606,
        0.0871293 ,  0.0202184 , 10.12      ,  0.18      ,  0.28      ])

>>> T2
array([ 0.9       , -1.19      ,  0.7       ,  0.78052918,  0.11827443,
        0.63992102,  0.93      , -1.17      ,  0.73      ,  0.41466194,
        0.26455561,  0.77423369,  0.45615033,  0.56843395,  0.0187898 ,
        0.6176355 ,  0.61209572,  0.86      , -1.175     ,  0.73      ])

>>> T3
array([0.3       , 0.8       , 0.1       , 0.06022547, 0.66676672,
       0.67063787, 0.1       , 0.6       , 0.2       , 0.36371077,
       0.57019677, 0.43860151, 0.98837384, 0.10204481, 0.20887676,
       0.16130952, 0.65310833, 0.26      , 0.83      , 0.25      ])

Let's plot these three time series:

NOTE: we are focusing on non-normalized Euclidean Distance

T_3D = np.transpose(np.c_[T1, T2, T3])
m = 3

---

Can we see any motifs? Let's use "multi-dimensional matrix profile" and "elbow method" to find the correct subspace that contain the motifs:

It seems we should better not to consider all 3 dimensions together. Right?

---

Let's calculate one-dimensional matrix profile for T_3D

P = np.full(l, np.inf, np.float64)
I = np.full(l, -1, np.int64)

excl_zone = math.ceil(m/4)
for i in range(l):
    D_square = np.zeros(l, dtype=np.float64)
    for j in range(3):
        D_square += np.square(core.mass_absolute(T_3D[j, i:i+m], T_3D[j]))
    
    D = np.sqrt(D_square)
    core.apply_exclusion_zone(D, i, excl_zone, np.inf)
    
    I[i] = np.argmin(D)
    P[i] = D[I[i]]
    
    if ~np.isfinite(P[i]):
        I[i] = -1

Let's take a look at matrix profile:

And, let's take a look at the discovered motif pairs (shown with thicker lines) :

Now, is it useful? Well, maybe! it depends on the data and application, and the domain expert to see what it means. What I am saying is that we cannot ignore the beauty of having 1-dimensional matrix profile for multi-dimensional time series.

Advantages:

• Can be easily implemented by modifying the current modules
• Easy interpretation
• Might give new insight
• As a user, I might be interested to see if there is any motifs across ALL dimensions.

---

Can we calculate 1-dimensional matrix profile for each time series data separately and get the same insight?
Let's take a look:

mp1 = stumpy.stump(T1, m, normalize=False)
mp2 = stumpy.stump(T2, m, normalize=False)
mp3 = stumpy.stump(T3, m, normalize=False)

# according to mp1: motif at index 0, and its NN at index 17
# according to mp2: motif at index 0, and its NN at index 6
# according to mp3: motif at index 6, and its NN at index 15

side-note: Are these 1-dimensional matrix profiles useful? Again, that depends on the application. Maybe a user wants to analyze each dimension separately.

[NimaSarajpoor]

Can you try the MDL method instead of using the elbow method to see if it makes a difference?

For some reason, one of the output that I am getting from one of the functions is empty. I am still working on it. Will update you if I get to some result/ conclusion

[NimaSarajpoor]

Can you try the MDL method instead of using the elbow method to see if it makes a difference?

For some reason, one of the output that I am getting from one of the functions is empty. I am still working on it. Will update you if I get a result

[NimaSarajpoor]

I've been thinking about this a little bit more and I think we could define a "transformation function" that accepts a multi-dimensional distance profile and returns a transformed multi-dimensional matrix profile.

This is cool :) and provides a good flexibility for the user!

[seanlaw]

This is cool :) and provides a good flexibility for the user!

Yeah, I am only interested in coding up what the original papers did so as to remain as faithful as possible. However, this doesn't mean that we can't develop/design things in a modular/replaceable way. Having said that, we should not encourage anything that doesn't work well (i.e., a new custom apply function that people may or may not use). Otherwise, we'll be stuck supporting it forever. So, before officially adding the feature, we should find a REAL dataset that shows that your "apply function" would be useful. By "REAL", I mean a non-contrived dataset that has a larger window and isn't made up data.

[NimaSarajpoor]

Yep..that makes sense! I need to go through the multi dimensional matrix profile tutorial (again), and then I will apply that as well as my own func to see what I can get for a real data. I will let you know if I find something interesting.