Add Top-K Nearest Neighbors to Matrix Profile (normalize=True)

[seanlaw]

Currently, functions like stump, stumped, and gpu_stump only return the top-1 nearest neighbor. Depending on the added complexity, we may want to consider adding a top-K nearest neighbor matrix profile.

This would add additional columns to the output. Namely, it would add 2 additional columns for every nearest neighbor.

Also, we should use np.searchsorted in place of heapq

See: #639 and #640 as they are closely related to this issue

[NimaSarajpoor]

@seanlaw
I am going through the stump right now. I am proposing the following design:

We can create a new function stump_topk, and then:

# in stump
def stump(...., k=1):
    if k>1:
        return stump_topk(...)

    # stump code

And, the output mp in mp = stump_topk(...) does not have four columns. Instead, it has 2k columns:
P_topk = mp[:, :k]
I_topk = mp[:, k:]

The reason is that I think returning the right and left NN may not be necessary in matrixprofile_topk, and this rerouting to another function can result in cleaner code.

What do you think? Or, do you prefer to just apply changes to the current function, and be worried about the design later?

[NimaSarajpoor]

@seanlaw
I think I should better submit a PR first (in which I just modify the current module without creating new function). Then, you can provide your comments there. I will update you.

[seanlaw]

Or, do you prefer to just apply changes to the current function, and be worried about the design later?

I think I would prefer trying to modify the existing functions by adding the relevant and needed data structures. If you add a stump_topk function then it will cause too much code branching and much of the code will be highly redundant. This means that if we change _stump then we'll likely need to change stump_topk which is undesirable. Right now, when we change the normalized code, we also need to change the non-normalized code but this was a specific and necessary design choice that has given us more flexibility. In this case, I think stump_topk would lead to more branching and more work in the long run.

To be clear, I don't know what the best solution is so we will have to be flexible as the right solution isn't clear to me either. I have a few points:

1. Before we do any coding, we need to decide on the output data structure. Whatever we decide should also minimize any change to the existing code logic where possible and affect users the least. Again, I don't know the right answer but I suspect that we should append the additional 2 * k - 2 columns (k-1 distances and k-1 indices) to the existing 4 columns (e.g., mp distances, mp indices, left indices, right indices) of the current matrix profile. This way, if somebody accidentally sets k>1 then the first 4 columns are still the same. The only difference would be the addition of two columns. Of course, it'll be a bit of bookkeeping to select for the top-k values from this array but maybe it's fine.
2. After we've decided on the output, we should be smarter and actually follow "Test Driven Development", whereby we write our unit test first along with the expected output. This is often referred to as red, green, refactor. It is important that we DON'T change any other tests and we only write one new test. When we run this new test, the test should fail because the output will not match our test expectations. Then, we modify our code (not our tests) to make our test(s) pass. So, the first commit should be ONLY the new unit test. Sure, you can add a stump_topk function to naive.py if you like.
3. Finally, as per above in 2., you modify the existing functions and we try to learn what will work and what will not work and really think through the design.

Any thoughts/questions? Criticism/feedback?

[NimaSarajpoor]

I would prefer trying to modify the existing functions by adding the relevant and needed data structures. If you add a stump_topk function then it will cause too much code branching and much of the code will be highly redundant. | To be clear, I don't know what the best solution is so we will have to be flexible as the right solution isn't clear to me either.

Right! In fact, I was thinking of creating new private functions as well. That may become a disaster :) I am flexible. So, I will modify the stump function as it seems to be a more rational solution at this moment.

1. This way, if somebody accidentally sets k>1 then the first 4 columns are still the same. The only difference would be the addition of two columns.

I overlooked this point! I will do as suggested and we can then see if it sounds alright.

2. This is often referred to as red, green, refactor. It is important that we DON'T change any other tests and we only write one new test. When we run this new test, the test should fail because the output will not match our test expectations.

Thanks for the link! Very interesting!! I think I got it.

Just to confirm:
stump calls the private functions _stump, which calls another private function _compute_diagonal. And, if I understand correctly, the unit test is only written for stump. So, I assume the new unit test should be written for stump only (and not for those two private functions). Correct?

[seanlaw]

stump calls the private functions _stump, which calls another private function _compute_diagonal. And, if I understand correctly, the unit test is only written for stump. So, I assume the new unit test should be written for stump only (and not for those two private functions). Correct?

For now, let's go with "yes" :) What might fail is 100% code coverage and then we'll have to make sure we handle that. I don't want to get too far ahead of ourselves and over-analyze everything. Given the possible complexities and newness, let's just "try something and see what happens". Maybe we'll break something and we need to either fix it or abandon this altogether. Either way, we learn something new. Sounds good?

At the end of the day, I believe that only ONE user has asked about VALMOD and maybe only 2 more have asked about top-k, which isn't a lot. If it's not a big problem to add top-k then cool because, computationally, it shouldn't be a huge additional expense! However, if it really makes our code unmaintainable then I don't know how much sense it would make to keep it given that there's not a ton of demand for it. When I develop, I want to make sure that I am focused on adding the things that are the most useful to the most users. Often times, it is adding documentation and tutorials that has the largest impact and, rarely, is it adding a new feature. It doesn't make much sense solving a problem for a party of one when we have limited time and resources. We need to be pragmatic and perfect the things that are actually worth perfecting and that affects the largest number of people.

[JaKasb]

k-NN for each distance-matrix row ?
k-NN for each distance-matrix column ?
k-NN for each motif ?

The most useful method is k-NN for each motif.
It is also the fastest method.
Furthermore it is also easy to implement.
But it won't compute a true k-NN neighborhood search, because there is no neighborhood information for time-points which were not classified as motifs.

1. compute regular matrix profile
2. perform peak detection on negative distance profile
   (scipy signal find_peaks, enforce minimum temporal distance between peaks -> exclusion-zone width)
   https://docs.scipy.org/doc/scipy-1.8.0/html-scipyorg/reference/generated/scipy.signal.find_peaks.html
3. perform MASS(time_series[x:x+w], time_series) for each x in peak_indices
4. for each MASS output: np.argsort -> computes the location and similarity-ordering of neighbors.
5. for each np.argsort(mass_result) : select the first k-entries
6. Despair, because now you have to deal with trivial-matches of order k instead of order 1.

Alternative Algorithm, entirely based on Black Magic:
https://github.com/xiaoshengli/LL

[seanlaw]

@JaKasb We are talking about distance matrix rows in this case. Given that all pairwise subsequence distances are computed, I suspect that this is more about whether there is a clean way to handle the bookkeeping that also doesn't reduce the maintainability of the code base. Any pointers or suggestions is greatly appreciated!

I have not come across this linear time paper. Thanks for sharing it!

[NimaSarajpoor]

(@seanlaw You beat me to it!)

@JaKasb
Thanks for sharing your knowledge and the paper!

k-NN for each distance-matrix row ? k-NN for each distance-matrix column ?

Are these two different? I mean...the distance matrix is symmetric. If there is a difference, could you please explain what that is?

The most useful method is k-NN for each motif.

I can see your point! As you noticed, top-k smallest values in each distance profile may not be useful for getting k-NN best match to each subsequence. Because, it does not exclude trivial matches. Therefore, after removing trivial matches, we may end up with a fewer number than k. The goal of this top-k smallest values in each distance profile is that we want to use them later in the VALMOD algorithm (see PR #586). And, I think the VALMOD algorithm does not work properly if we do not return top-k smallest values of each distance profile (NOT top-k smallest with minimum temporal distance).

---

I am trying to understand the thing you suggested

a true k-NN neighborhood search, because there is no neighborhood information for time-points which were not classified as motifs."

Could you please elaborate/clarify this sentence?

2. perform peak detection on negative distance profile

Did you mean "matrix profile" here?

5. for each np.argsort(mass_result) : select the first k-entries
6. Despair, because now you have to deal with trivial-matches of order k instead of order 1.

So, I assume you want to get top-k best matches (considering minimum temporal distance between these top-k matches) for each x in peak_indices. Is that correct?

And then, at the end, we have M motifs (# of peaks detected from negative of matrix profile, i.e. M smallest values considering the exclusion of trivial matches) and, for each, we find top-k matches. Correct? And, then there should be a way to sort these motifs considering their k-NN. (maybe average to their k neighbors?)

[NimaSarajpoor]

@seanlaw

Sounds good?

Yes! That makes sense!

We need to be pragmatic and perfect the things that are actually worth perfecting and that affects the largest number of people.

Yep! Let's go with MVP and see how things will go!

[NimaSarajpoor]

I am trying to write test_stump_topk.py which is basically the same as test_stump.py. I just need to replace naive.stump with naive.stump_topk (a new function that I added to naive.py).

However, I can also see that test_stump.py calls the function naive.stamp in some of its tests, and that means I need to create naive.stamp_topk as well. Is that correct?

---

SIDE NOTE:
naive.stamp calls naive.mass function, which is completely different than core.mass.
Why not using the name mass_PI, the name that is used in stumpy.stamp?

[NimaSarajpoor]

@seanlaw
I decided to submit a PR so you can see if I got your point correctly or not. For now, I skipped the use of naive.stamp in the tests.

[seanlaw]

I am trying to write test_stump_topk.py which is basically the same as test_stump.py. I just need to replace naive.stump with naive.stump_topk (a new function that I added to naive.py).

Actually, I don't think you need naive.stump_topk. I think you should modify naive.stump so that it can accept an additional parameter k and then update all of the functions that naive.stump depends on.

However, I can also see that test_stump.py calls the function naive.stamp in some of its tests, and that means I need to create naive.stamp_topk as well. Is that correct?

At the end of the day, we'll call stump.stump with k > 1 and so we'll want a naive (i.e., slow) way to return back the correct result.

Why not using the name mass_PI, the name that is used in stumpy.stamp?

Yes, I don't recall "why" (maybe due to laziness or haste) but, indeed, we should have called it mass_PI instead. In fact, we should've split the mass part out into naive.mass and then created naive.mass_PI that calls naive.mass. Would you mind doing this and then updating any tests that use naive.mass?