Updated derivation with improved numerical stability

[seanlaw]

Pull Request Checklist

Below is a simple checklist but please do not hesitate to ask for assistance!

• Read our Contributing Guide
• Referenced a Github issue (or create one if one doesn't already exist)
• Read and reviewed all of the comments in the Github issue that you've referenced (along with cross-referenced issues/pull requests) to ensure that the issue still requires a pull request
• Checked that the issue has not already been assigned to anybody else or is already being addressed in another pull request
• Left a meaningful comment on the original Github issue to discuss the detailed approach for your contribution and received confirmation from the maintainers before proceeding with this pull request
• Forked, cloned, and checked out the newest version of the code
• Created a new branch
• Made necessary code changes
• Installed black (i.e., python -m pip install black or conda install -c conda-forge black)
• Installed flake8 (i.e., python -m pip install flake8 or conda install -c conda-forge flake8)
• Installed pytest-cov (i.e., python -m pip install pytest-cov or conda install -c conda-forge pytest-cov)
• Ran black --exclude=".*\.ipynb" --extend-exclude=".venv" --diff ./ in the root stumpy directory
• Ran flake8 --extend-exclude=.venv ./ in the root stumpy directory
• Ran ./setup.sh dev && ./test.sh in the root stumpy directory and ensured that all tests are passing locally
• Check this box if AI code assistance was used to generate 15%+ of the code in this pull request

Please do not commit any code to avoid/circumvent a failing test and, instead, engage in a discussion (below) to determine the best course of action.

Only request a review after the checklist above is fully completed!

[gitnotebooks]

Found 1 changed notebook. Review the changes at https://app.gitnotebooks.com/stumpy-dev/stumpy/pull/1180

[seanlaw]

@NimaSarajpoor Would you mind taking a look when you have some time? Thanks in advance!

You may need to download the file in order to get the LaTex to render properly

[NimaSarajpoor]

@seanlaw
It was great to see all the steps involved in deriving equation (91). Thanks for putting them together so clearly! I have shared a few comments for your consideration.

[NimaSarajpoor]

It should be worth it to better demonstrate/explain the benefit of achieving this new equation. We can show that the new equation (91) has LARGER UPDATE and therefore should be more stable as it has less catastrophic cancelation.

Let's start with the equation (40), i.e.

$$cov(Q_{i,m}, T_{i,m}) = cov(Q_{i-1,m}, T_{i-1,m}) + update$$

Let's multiply both sides by m:

$$ m * cov(Q_{i,m}, T_{i,m}) = m * cov(Q_{i-1,m}, T_{i-1,m}) + m * update $$

Recall that m * cov is basically the summation part of "centered sum-of-products". So, this means the update term in the equation (91), i.e. $df^{T}.dg^{Q} + df^{Q}.dg^{T}$, is equivalent to m * update (NOT update). This quickly shows that the new equation has LARGER UPDATE.

[NimaSarajpoor]

Let's see what else we can get from the relationship between equation (40) and (91). Let's start from equation (40) again:

$$cov(Q_{i,m}, T_{i,m}) = cov(Q_{i-1,m}, T_{i-1,m}) + \frac{m-1}{m^{2}}U$$

where

$$ U = [(q_{i+m-1} - \mu_{Q_{i,m-1}})(t_{i+m-1} - M_{T_{i,m-1}}) - (q_{i-1} - \mu_{Q_{i,m-1}}) (t_{i-1} - M_{T_{i,m-1}})]$$

---

Let's multiply both sides of equation (40) by $m$:

$$m * cov(Q_{i,m}, T_{i,m}) = m * cov(Q_{i-1,m}, T_{i-1,m}) + \frac{m-1}{m}U$$

Again, recall that m * cov is basically the summation part of "centered sum-of-products". So, if I compare the equation above with the equation (91), I can see:

$$ \frac{m-1}{m}U = df^{T}.dg^{Q} + df^{Q}.dg^{T} $$

Note that $\frac{m-1}{m} < 1$. If $U$ is small, then $\frac{m-1}{m}U$ is smaller, meaning the value on the right hand side of the equation above is smaller than $U$. Therefore, if we say that $U$ involves the subtraction of two nearby numbers and that can lead to catastrophic cancellation, then, in that scenario, we will have "catastrophic cancellation" on the right hand side too as that is even smaller!

Can we simply use the update term $\frac{m-1}{m}U$ when we use the summation part of "centered sum-of-products"?

[seanlaw]

Maybe it's an issue on my end when things are being rendered in my Jupyter notebook but I'm only seeing a single equation number for a block of related/multiline equations:

If you are seeing the same thing, maybe we can refer to the lines as (14.1) - (14.5) (to mean lines 1-5 in equation block 14)? Otherwise, I am counting equations one-by-one each time :(

Or maybe you have a different workflow?

[seanlaw]

Can we simply use the update term $\frac{m−1}{m}U$ when we use the summation part of "centered sum-of-products"?

I've tried reading your comments a few times but I am unable to get your point.

Therefore, if we say that $U$ involves the subtraction of two nearby numbers and that can lead to catastrophic cancellation, then, in that scenario, we will have "catastrophic cancellation" on the right hand side too as that is even smaller!

I feel like I'm only reading this as a statement of fact but I am getting the sense that you have a much stronger point/argument that I am overlooking :(

You started with:

It should be worth it to better demonstrate/explain the benefit of achieving this new equation.

Maybe I've been staring at the equations too long but I'm not able to see "how" this is demonstrating the benefit any further (i.e., what is the special insight that "proves" the benefit clearly). I know you so I can feel that there is something (intuition?) there. I'm just not smart enough to infer your intuition but I would like to hear more.

[NimaSarajpoor]

Or maybe you have a different workflow?

I didn't download the notebook. I just opened the file on Github from your branch. It renders the file differently and I didn't know that it affects the equations numbers! Apologies for the inconvenience.

Please allow me to download the updated notebook and read it again to see the added information. I will use the equation numbers from there.

[seanlaw]

Hmm, I see your point and do not yet have a good answer...

[NimaSarajpoor]

I have been thinking about conducting an experiment that can demonstrate empirically which option works better when updating $\overline{QT}$.

If we can compute the actual values directly with higher precision, we can see which option works better.

What do you think?

[seanlaw]

What do you think?

I am certainly not opposed to it. Perhaps, the statement about "catastrophic cancellation" was not correct and that language can be softened. From a PR perspective, at what point is the derivation "correct" and what we need to do is adjust the language (maybe as a separate PR)?

[NimaSarajpoor]

Perhaps, the statement about "catastrophic cancellation" was not correct

Probably...If "a - b" has catastrophic cancellation, then the issue might be addressed if we change it to subtraction-free operations (with positive-only elements). However, this is not achieved here.

From a PR perspective, at what point is the derivation "correct"

I think the moment we obtained the original correction term, i.e. $[(q_{i+m-1} - \mu_{Q_{i,m-1}})(t_{i+m-1} - M_{T_{i,m-1}}) - (q_{i-1} - \mu_{Q_{i,m-1}}) (t_{i-1} - M_{T_{i,m-1}})]$ . We can then show that:

$$\overline{QT_{i}} = \overline{QT_{i-1}} + \frac{m-1}{m} \times correction$$

And that's it (unless we can show the reformulated version is better)

---

But, something is still bothering me. If the equation above is better (i.e., numerically more stable), then what about the following?

$$m * \overline{QT_{i}} = m * \overline{QT_{i-1}} + (m-1) \times correction$$

So, we compute $m * \overline{QT_{i}}$ and keep updating it by using $(m-1) \times correction$ (?!?)

[NimaSarajpoor]

[WIP]
Ok... I did an experiment to check which of the following equations results in less numerical error for $\overline{QT}$

Recursion Relation Using Original Formula

$$ \overline{QT_{i}} = \overline{QT_{i-1}} + \frac{m-1}{m} [(q_{i+m-1} - \mu_{Q_{i,m-1}})(t_{i+m-1} - M_{T_{i,m-1}}) - (q_{i-1} - \mu_{Q_{i,m-1}}) (t_{i-1} - M_{T_{i,m-1}})]$$

Recursion Relation Using New Formula

$$ \overline{QT_{i}} = \overline{QT_{i-1}} + df^{T}.dg^{Q} + df^{Q}.dg^{T} $$

Ground Truth

$$ \overline{QT_{i}} = (T[i:i+m] - M_{T_{i,m}}).(Q[i:i+m] - \mu_{Q_{i,m}})$$

Experiment

inputs
n = 2 ** 20
np.random.seed(0)
T = np.random.rand(n) * scale
Q = np.random.rand(n) * scale

# m: 10, 100, 1000
# scale: 0.001, 1, 1000

We compute the $\overline{QT}$ for i in range(n-m+1). We then compute the error for original and new formula, and compare their errors.

Code

# imports
import matplotlib.pyplot as plt
import numpy as np
from stumpy import core

# Ground Truth
def _mcov_precise(Q, T, M_T, μ_Q, i):
    mcov = np.dot(
            (T[i : i + m] - M_T[i]),
            (Q[i : i + m] - μ_Q[i]),
        )

    return mcov   


def mcov_precise(Q, T, m):
    if len(Q) != len(T):
        msg = 'The length of two arrays should be the same'
        raise ValueError(msg)
    M_T, _ = core.compute_mean_std(T, m)
    μ_Q, _ = core.compute_mean_std(Q, m)

    l = len(T) - m + 1
    out = np.empty(l, dtype=np.float64)
    for i in range(l):
        out[i] = _mcov_precise(Q, T, M_T, μ_Q, i)

    return out
    

## Original Formula
def _mcov_sliding_original(mcov, constant, cov_a, cov_b, cov_c, cov_d, i):
    mcov_next = mcov + constant * (
        cov_a[i] * cov_b[i]
        - cov_c[i] * cov_d[i]
    )

    return mcov_next


def mcov_sliding_original(Q, T, m):
    if len(Q) != len(T):
        msg = 'The length of two arrays should be the same'
        raise ValueError(msg)

    (
        Q, # T_A
        μ_Q,
        σ_Q_inverse,
        μ_Q_m_1,
        Q_subseq_isfinite,
        Q_subseq_isconstant,
    ) = core.preprocess_diagonal(Q, m, T_subseq_isconstant=None)
    
    (
        T, # T_B
        M_T,
        Σ_T_inverse,
        M_T_m_1,
        T_subseq_isfinite,
        T_subseq_isconstant,
    ) = core.preprocess_diagonal(T, m, T_subseq_isconstant=None)
    
    
    # compute for elements for updating cov or mcov
    cov_a = T[m - 1 :] - M_T_m_1[:-1]
    cov_b = Q[m - 1 :] - μ_Q_m_1[:-1]
    
    cov_c = np.empty(M_T_m_1.shape[0], dtype=np.float64)
    cov_c[1:] = T[: M_T_m_1.shape[0] - 1]
    cov_c[0] = T[-1]
    cov_c[:] = cov_c - M_T_m_1
    
    cov_d = np.empty(μ_Q_m_1.shape[0], dtype=np.float64)
    cov_d[1:] = Q[: μ_Q_m_1.shape[0] - 1]
    cov_d[0] = Q[-1]
    cov_d[:] = cov_d - μ_Q_m_1


    constant = (m-1) / m
    l = len(T) - m + 1
    out = np.empty(l, dtype=np.float64)
    out[0] = np.dot(T[:m] - M_T[0], Q[:m] - μ_Q[0])
    for i in range(1, l):
        out[i] = _mcov_sliding_original(
            out[i-1], 
            constant, 
            cov_a, 
            cov_b, 
            cov_c, 
            cov_d, 
            i
        )

    return out


# New Formula
def _mcov_sliding_new(mcov, Q, T, M_T, μ_Q, i):
    mcov_next = mcov + 1/2 * (T[i+m-1] - T[i-1]) * ((Q[i+m-1]-μ_Q[i]) + (Q[i-1] - μ_Q[i-1])) + 1/2 * (Q[i+m-1] - Q[i-1]) * ((T[i+m-1]-M_T[i]) + (T[i-1] - M_T[i-1]))

    return mcov_next


def mcov_sliding_new(Q, T, m):
    
    if len(Q) != len(T):
        msg = 'The length of two arrays should be the same'
        raise ValueError(msg)
    M_T, _ = core.compute_mean_std(T, m)
    μ_Q, _ = core.compute_mean_std(Q, m)

    l = len(T) - m + 1
    out = np.empty(l, dtype=np.float64)
    out[0] = np.dot(T[:m] - M_T[0], Q[:m] - μ_Q[0])
    for i in range(1, l):
        out[i] = _mcov_sliding_new(
            out[i-1], 
            Q, 
            T, 
            M_T, 
            μ_Q, 
            i
        )

    return out

Results