u-normalized version of the Matrix Profile

[bjlobo]

I know that when normalize=True, stumpy computes the z-nomalized distance between subsequences where z-normalized means that each subsequence has it's mean subtracted before dividing by it's standard deviation before the distance is computed.

Is it possible in stumpy to specify a version where just the mean is subtracted (what I am terming u-normalization because you are only subtracting the mean) before computing the distance?

[seanlaw]

@bjlobo Thank you for your question and welcome to the STUMPY community. Unfortunately, what you are asking for is not possible in STUMPY. While it might seem trivial on the surface, STUMPY employs many complex mathematical manipulations and computational tricks behind the scenes to ensure that the matrix profile is computed as efficiently as possible. This means that there is only a certain class of distance functions that we can accommodate that can allow us to traverse the distance matrix in a diagonal fashion (rather than row-wise or column-wise) and leverage a hidden relationship between the diagonal elements, which is where the time savings comes from.

[seanlaw]

@NimaSarajpoor I haven't thought about this thoroughly but I wonder if an ugly hack would be to simply set σ_Q_inverse[:] = 1 and Σ_T_inverse[:] = 1 in the stumpy.stump function?

[bjlobo]

@seanlaw Understood. I did not think it was trivial, but thought it might be possible.

My motivation for the question is that I don't care about the vertical locations of two time series subsequences matching (i.e., where on the y-axis the time series subsequence is located) which is taken care of by dividing by the means, but I do care about the (unscaled) shapes matching (and division by the standard deviations scales the time series subsequences).

The figure below shows two time series subsequences which I would NOT want to match (and in this case the subsequences do NOT match under no nomalization (leftmost plot) but do match under z-normalization (rightmost plot), and hopefully would NOT match under u-normalization (middle plot)).

The figure below shows two time series subsequences which I WOULD want to match (and in this case the subsequences do NOT match under no nomalization but do match under z-normalization, and hopefully WOULD match under u-normalization).

[NimaSarajpoor]

@seanlaw
I thought about the hack but I think it does not work. The reason that I am saying this is because z-norm distance has linear relationship with $\sqrt{1 - \rho}$, where $\rho = \frac{cov}{m\sigma\Sigma}$. Therefore, setting $\sigma$ and $\Sigma$ to one means that we are setting the denom of the latter equation to m and the $\rho$ may exceed the range [-1, 1]. And, that means the distance may become undefined since the value under square-root sign can become negative.
[Additional note: Setting $\sigma$ and $\Sigma$ to one means that we are treating the subsequence T and the subsequence 2T the same since we are considering std value of one for both time series. However, IIUC, the u-normalized distance wants to make a difference between those two cases]

---

I think we should start with this STUMPY document which provides step-by-step proof on the relationship between z-norm distance and $\rho$.

When I start from the top and go through the equations to derive a formula for u-normalization, I will get the u-normalized distance provided below:

$d = \sqrt{2m(\sigma\Sigma - \rho')}$, where $\rho' = \frac{cov}{m}$

Now, I think stumpy.stump has necessary components of the formula shown above. So, with some changes, we should be able to get the distance above. Note that, in contrast to z-norm where $\rho$ values are compared first, the u-normalize distance needs to be computed and compared on the fly.

[NimaSarajpoor]

As a follow up to my previous comment:

Let's expand the parentheses by multiplying m to each term inside that parentheses:
$d = \sqrt{2(m\sigma\Sigma - cov)}$

And, now we can say:
$d = \sqrt{2r}$, where $r = m\sigma\Sigma - cov$.

So, in stumpy.stump, instead of

stumpy/stumpy/stump.py

Lines 205 to 206 in 6663f5f

	pearson = cov * Σ_T_inverse[uint64_j] * σ_Q_inverse[uint64_i]
	pearson = min(1.0, pearson)

,

we do:

r = m * Σ_T_inverse[uint64_j] * σ_Q_inverse[uint64_i] - cov 

In contrast to ρ, we are interested in lower value of r . we should be careful in updating the best-so-far r. For example, in stumpy.stump, we have the following line:

stumpy/stumpy/stump.py

Line 212 in 6663f5f

if pearson > ρ[thread_idx, uint64_i, 0]:

To update r in similar manner, we should use < (less than).

How about constant subsequences? I think one needs to think about that too and see if necessary changes need to be applied there or not for u-normalized distance.

[seanlaw]

For all of these reasons, I would NOT recommend doing this in STUMPY as it adds a tonne of complexity. @bjlobo If you are lucky enough to have a time series that is, say, less than 10K data points in length, then you might be able to get away with a poor-man's version of what you want:

from stumpy import config, core
import numpy as np
from scipy.spatial.distance import cdist
from numba import njit

T = np.random.rand(1000)
m = 50

@njit(fastmath=True)
def custom(u, v):
    return np.linalg.norm((u-u.mean()) - (v-v.mean()))

distance_matrix = cdist(core.rolling_window(T, m), core.rolling_window(T, m), metric=custom)

excl_zone = int(np.ceil(m / config.STUMPY_EXCL_ZONE_DENOM))
mask = (
    np.tri(distance_matrix.shape[0],distance_matrix.shape[1], k=excl_zone, dtype=bool) & 
    ~np.tri(distance_matrix.shape[0],distance_matrix.shape[1], k=-excl_zone, dtype=bool)
)
distance_matrix[mask] = np.inf

np.min(distance_matrix, axis=1)

[bjlobo]

@seanlaw Thanks for the code. Fortunately my problem does have less than 10k data points for each time series so I am hopeful that this solution will work.

[vlad-user]

@bjlobo, @seanlaw I'm looking to apply a similar normalization approach for approximately 3,000 data points. However, I'm in search of an implementation for streaming data, akin to how stumpy.stumpi is structured. Would it be challenging for me to adapt the existing implementation of stumpy.stumpi to accommodate this requirement?

[seanlaw]

@vlad-user I encourage you to try it and, if you get stuck, then please start a new discussion, reference this one, show your code, and we will do our best to help. With only 3000 data points, you are better off with:

1. Compute the (custom) matrix profile using your custom distance function
2. Allow a new data point to arrive
3. Compute the distance profile (this is simply the distance between the newest subsequence and the rest of the time series) for the new incoming data point using your custom distance function
4. Update the matrix profile