Add MASS V4 Challenger

[seanlaw]

While looking up the original MASS code, I noticed a new MASS V4. It claims to be:

1. Numerically stable (especially for longer time series) and without producing imaginary numbers
2. Exploits Direct Cosine Transform (DCT), which is suppose to make MASS faster

% This code is created by Sheng Zhong and Abdullah Mueen
% The overall time complexity of the code is O(n log n). The code is free to use for research purposes.
% The code does not produce imaginary numbers due to numerical errors 
% k should greater than or equals to floor((3m+1)/2)

function [dist] = MASS_V4(T, Q, k)
    n = length(T);
    m = length(Q);
    Q = zNorm(Q);
    dist = [];
    batch=get_batch_size(k, m);
    
    for j = 1 : batch-m+1 : n-m+1
        right = j + batch - 1;
        if right >= n
            right = n;
        end 
        dot_p = DCT_dot_product(T(j:right), Q);
        sigmaT = movstd(T(j:right),[m-1 0],1);
        d = sqrt(2*(m-(dot_p./sigmaT(m:end)))) ;        
        dist = [dist; d];
    end
    
end

function dct_batch_size = get_batch_size(target_batch_size, m)
    dct_batch_size = floor((2*target_batch_size-2)/3) -1;
    if dct_batch_size < m
        dct_batch_size = m;
    end
    pad_length = dct_batch_size + floor((dct_batch_size-m+1)/2) + floor((m+1)/2);
    while pad_length < target_batch_size
        dct_batch_size = dct_batch_size + 1;
        pad_length = dct_batch_size + floor((dct_batch_size-m+1)/2) + floor((m+1)/2);
    end
    if (pad_length > target_batch_size)
        dct_batch_size = dct_batch_size - 1;
    end
end


function dot_p = DCT_dot_product(x, y)
    n = length(x);
    m = length(y);
    [x_pad, y_pad, si] = DCT_padding(x,y);
    
    N = length(x_pad);
    
    Xc = dct(x_pad, 'Type', 2);
    Yc = dct(y_pad, 'Type', 2);
    
    dct_product = Xc .* Yc;
    dct_product(N + 1) = 0;
    dct_product(1) = dct_product(1) * sqrt(2);
    dot_p = dct(dct_product, 'Type', 1);
    dot_p(1) = dot_p(1) * 2;
    dot_p = sqrt(2*N) * dot_p(si: si+n-m);
end

function [x_pad, y_pad, start_index] = DCT_padding(x, y)
    n = length(x);
    m = length(y);
    p2=floor((n-m+1)/2);
    p1=p2+floor((m+1)/2);
    p4=n-m+p1-p2;
    x_pad = zeros(n+p1, 1);
    x_pad(1+p1:end) = x;
    y_pad = zeros(m+p2+p4, 1);
    y_pad(1+p2:1+p2+m-1) = y;
    start_index = p1-p2+1;
end

[NimaSarajpoor]

@seanlaw
Noticed that when I was looking at their webpage. Should have created an issue for it! Thanks for creating this! 👍

[NimaSarajpoor]

The MASS webpage also point to the authors' MASS paper published in 2024. See section 4.4 to read more about "MASS V4".

[seanlaw]

V4 executes 46% to 48% more FLOPs than V3 due to the extra padding with zeros and most importantly, the DCT always executes more FLOPs than FFT (Johnson and Frigo 2007; Shao and Johnson 2008)

More FLOPs generally means "faster" but, in this case, they mean "Total number of floating point operations to compute the sliding dot product" (so "more FLOPs" is actually "slower" - i.e., a smaller number is faster in the table below)!

Based on this, V3 is the fastest and then followed by V4. V1 and V2 are considered "FFT"/"DFT" based methods and slightly slower. One thing to note is that there is no guarantee that STUMPY has implemented V2 because we haven't directly performed zero-padding on Q so that it matches the length of T. We've simply left everything up to the scipy.signal.convolve function. I wonder if we can achieve the same/similar 2x speed up observed in MASS_V2 by:

1. Manually zero-padding Q to be the same length as T
2. Apply pocketfft rfft/irfft ourselves
3. I wonder if this will match the performance of FFTW?

Note that while MASS_V3 is the fastest, MASS_V4 offers much better stability (no imaginary numbers) and it is still faster than FFT-based MASS methods. So, MASS_V4 is a reasonable compromise (assuming that it works) but we should also test the MASS_V2 with scipy to see if it helps at all

[seanlaw]

I wonder if we can achieve the same/similar 2x speed up observed in MASS_V2 by

Never mind. After some tinkering, I realized that we've already done exactly this in pocketfft_r2c_c2r_sdp

Nonetheless, I think MASS_V4 is still worth pursuing. Even if a SciPy-based MASS_V4 is still slower than FFTW/pyFFTW, it might still be faster than scipy.signal.convolve or pocketfft because it would be using a faster dct method. What do you think @NimaSarajpoor?

[NimaSarajpoor]

Never mind. After some tinkering, I realized that we've already done exactly this in pocketfft_r2c_c2r_sdp

Right! And, our version should be better as we use next_fast_len(n) instead of n

I think MASS_V4 is still worth pursuing. Even if a SciPy-based MASS_V4 is still slower than FFTW/pyFFTW, it might still be faster than scipy.signal.convolve or pocketfft

Right! Will check this out in the upcoming days...

[NimaSarajpoor]

Was hard for me to remember how the function DCT_padding (provided in this comment) constructs the padded Q and T. The following should help:

1. T_padded: $[0_{a}, 0_{b}, T]$
2. Q_padded: $[0_{a}, Q, 0_{b}, 0_{c}]$

where:

• $0_{a}$: $\lfloor \frac{n-m+1}{2} \rfloor$ numbers of zeros
• $0_{b}$: $\lfloor \frac{m+1}{2} \rfloor$ numbers of zeros
• $0_{c}$: n-m numbers of zeros

[NimaSarajpoor]

Maybe we should focus on V3... and only for large arrays? The section 4.5 of the paper says:

The performance show in Fig. 4, matches the observations made in Table 2. Overall, V3 is the fastest in all three metrics: time complexity, stopwatch time and FLOPs count.

According to the Fig. 4 provided in the paper, we can see that, for large input sizes, MASS_V3 is faster than MASS_V2, and MASS_V2 is faster than MASS_V4. Therefore, from performance perspective: V3 > V2 > V4

Notes:

1. The existing pyfftw_sdp.py module in this repo is actually "V2+", since we use next_fast_len(len(T)), which should be faster. So: V3 > V2+ > V2 > V4.... Well, actually V2+ (pyfftw_sdp) is in python and the others are in MATLAB. However, we did compare V2+ (pyfftw_sdp) with V2 (matlab) in Fixed #26 Use less-aggressive planning flag as default #27, and showed that "pyfftw_sdp" generally has better performance (see the plots in this comment and this comment)

2. For arrays in range 2^1...2^20, I can confirm that, on MATLAB, V2 outperforms V4 (batch logic excluded) (See the bottom of this comment)

3. The MASS webpage says:

The implementation [of MASS_V4] is a bit slower than MASS_V3, however, there is no numerical error due to leakage.

However, according to the Fig. 4. of the paper, MASS_V3 seems to be 2x faster than MASS_V4! So, IMO, that is NOT just "a bit..."

---

fft-based sdp vs dct-based sdp (on MATLAB)

We started to check the performance of DiscreteCosineTransform(DCT)-based SDP in PR #30 (with scipy) and in PR #31 (with pyfftw). We haven't yet added the "batch-processing" logic proposed in MASS_V4. Let's compare MASS_V2 with MASS_V4 (batch logic excluded) on MATLAB!

% DCT-based SDP


function [x_pad, y_pad, start_index] = DCT_padding(x, y)
    n = length(x);
    m = length(y);
    p2=floor((n-m+1)/2);
    p1=p2+floor((m+1)/2);
    p4=n-m+p1-p2;
    x_pad = zeros(n+p1, 1);
    x_pad(1+p1:end) = x;
    y_pad = zeros(m+p2+p4, 1);
    y_pad(1+p2:1+p2+m-1) = y;
    start_index = p1-p2+1;
end


function dot_p = dct_sdp(y, x)
    n = length(x);
    m = length(y);
    [x_pad, y_pad, si] = DCT_padding(x,y);
    
    N = length(x_pad);
    
    Xc = dct(x_pad, 'Type', 2);
    Yc = dct(y_pad, 'Type', 2);
    
    dct_product = Xc .* Yc;
    dct_product(N + 1) = 0;
    dct_product(1) = dct_product(1) * sqrt(2);
    dot_p = dct(dct_product, 'Type', 1);
    dot_p(1) = dot_p(1) * 2;
    dot_p = sqrt(2*N) * dot_p(si: si+n-m);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT-based SDP


function [z] = fft_sdp(y, x)
    %x is the data, y is the query
    m = length(y);
    n = length(x);

    y = y(end:-1:1);%Reverse the query
    y(m+1:n) = 0; %aappend zeros

    X = fft(x);
    Y = fft(y);
    Z = X.*Y;
    z = ifft(Z);
    z = z(m:n);

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% timing

function timing_arr = timing_fft_sdp(p_min, p_max, timeout)
    
    for q = p_min:p_max
        fprintf('--------------------\n');
        fprintf('q --> %d\n', q);
        
        m = 2^q;
        Q = rand(m, 1);
        T = rand(2^p_max, 1);
        
        fprintf('p --> ');
        for p = q:p_max
            T_p = T(1:2^p);
            fft_sdp(Q, T_p); % warm up
            
            timeout_tmp = timeout;
            if p >= 20
                timeout_tmp = 5.0;
            end

            total_time = 0;
            n_iter_count = 0;
            while total_time < timeout_tmp
                start_time = tic; % Start timer
                fft_sdp(Q, T_p);
                elapsed = toc(start_time); % Stop timer
                total_time = total_time + elapsed;
                n_iter_count = n_iter_count + 1;
            end
            
            avg_time =  total_time / n_iter_count;
            fprintf('%d, ', p);
            timing_arr(q, p) = avg_time;
        end
        fprintf('avg time for the last case: %d, ', avg_time);
        fprintf('\n');
    end
end



function timing_arr = timing_dct_sdp(p_min, p_max, timeout)
    
    for q = p_min:p_max
        fprintf('--------------------\n');
        fprintf('q --> %d\n', q);
        
        m = 2^q;
        Q = rand(m, 1);
        T = rand(2^p_max, 1);
        
        fprintf('p --> ');
        for p = q:p_max
            T_p = T(1:2^p);
            dct_sdp(Q, T_p); % warm up
            
            timeout_tmp = timeout;
            if p >= 20
                timeout_tmp = 2.0;
            end

            total_time = 0;
            n_iter_count = 0;
            while total_time < timeout_tmp
                start_time = tic; % Start timer
                dct_sdp(Q, T_p);
                elapsed = toc(start_time); % Stop timer
                total_time = total_time + elapsed;
                n_iter_count = n_iter_count + 1;
            end
            
            % Store the average time
            avg_time =  total_time / n_iter_count;
            fprintf('%d, ', p);
            timing_arr(q, p) = avg_time;
        end
        fprintf('avg time for the last case: %d, ', avg_time);
        fprintf('\n');
    end
end

%%%%%%%%%%%%%%%%

matlabVersion = version;
disp(['MATLAB Version: ', matlabVersion]);

num_of_thread = 'automatic'; 
maxNumCompThreads(num_of_thread);
n_threads = maxNumCompThreads;
disp(strcat("MAX N threads: ",num2str(n_threads)));
%%%%%%%%%%%%%%%%%%%%%%

p_min = 2;
p_max = 20;
timeout = 0.1;

disp('====================================');
disp('Timing fft_sdp');
timing_values_fft_sdp = timing_fft_sdp(p_min, p_max, timeout);

disp('====================================');
disp('Timing dct_sdp');
timing_values_dct_sdp = timing_dct_sdp(p_min, p_max, timeout);

save('matlab_sdp_p2_p20_with_fft_multithread.mat', 'timing_values_fft_sdp');
save('matlab_sdp_p2_p20_with_dct_multithread.mat', 'timing_values_dct_sdp');

And to plot based on the saved data:

import numpy as np
import matplotlib.pyplot as plt
import scipy


p_min = 2
p_max = 20

fft_sdp = scipy.io.loadmat('./fft_vs_dct/matlab_sdp_p2_p20_with_fft_multithread.mat')
fft_sdp = fft_sdp['timing_values_fft_sdp']


dct_sdp = scipy.io.loadmat('./fft_vs_dct/matlab_sdp_p2_p20_with_dct_multithread.mat')
dct_sdp = dct_sdp['timing_values_dct_sdp']

plt.figure(figsize=(10, 10))
plt.title(f'Performance Ratio of DCT-based SDP relative to FFT-based SDP')
ax = plt.gca()

for p1 in range(p_min, p_max + 1):
    for p2 in range(p1, p_max + 1):
        if p1 > p2:
            continue

        r = fft_sdp[p1-1, p2-1] / dct_sdp[p1-1, p2-1]

        plt.scatter(p1, p2, color='r')
        ax.annotate(
            f"{r:.2f}",
            (p1, p2),
            textcoords="offset points",
            xytext=(1, 1),
            ha="left",
            va="bottom"
        )

plt.xlabel('log2 of len(Q)')
plt.ylabel('log2 of len(T)')
plt.xticks(ticks=np.arange(p_min, p_max + 1), labels=np.arange(p_min, p_max + 1))
plt.yticks(ticks=np.arange(p_min, p_max + 1), labels=np.arange(p_min, p_max + 1))
plt.grid()
plt.show()

And here is the result: A value less than one means DCT-based SDP is slower:

[NimaSarajpoor]

UPDATE
After reading more about oaconvolve, I realized that the overlapping part is not happening at the input in overlap-add! The input is chunked into non-overlapping pieces, and the outputs will be overlapped and added. So, you may ignore the following. "MASS_V3" does overlap at input but OVERLAP-ADD does overlap at output.

---

As mentioned in my previous comment, The section 4.5 of the paper claims MASS_V3 is the fastest (particularly for large arrays).

After reading overlap-add method, it looks like both MASS_V3 and OVERLAP-ADD are basically following the same concept... Instead of computing (the dot product) for full array T, it breaks it down to smaller chunks and consider an overlap of m-1 to make sure all subsequences are covered. See Fig. 3 in this document

[seanlaw]

And here is the result: A value less than one means DCT-based SDP is slower:

So MASS_V4 is 8x slower than MASS_V2? Are you using the same code from the MASS website for both V2 and V4?

I guess this is consistent with V3 > V2 > V4. I misinterpreted Fig 4 in the paper and thought that V4 was faster than V2 so thanks for correcting my understanding!

Maybe we should focus on V3

Yes, I agree! I wonder if pocketfft + V3 could make it faster than FFTW + V2? I had always avoided V3` because it sounded like it was sensitive to or depended on the system it was running on

[NimaSarajpoor]

So MASS_V4 is 8x slower than MASS_V2?

Correct... for len(T)==2^20. Although the performance ratio depends on (len(Q), len(T)) , I noticed V4 is slower than V2 in all of those cases.

Are you using the same code from the MASS website for both V2 and V4?

correct. I took their SDP piece only.

I guess this is consistent with V3 > V2 > V4.

It is definitely consistent with V2 > V4 part. I didn't check V3 because of EXACTLY what you said:

I had always avoided V3 because it sounded like it was sensitive to or depended on the system it was running on

V3 requires an additional system-dependent parameter, the size of the chunk (See Fig.3 in the paper)

---

I wonder if pocketfft + V3 could make it faster than FFTW + V2?

it is worth considering. Not sure how it will look like for shorter arrays but we can explore and see. I am also interested in exploring scipy.signal._signaltools. _calc_oa_lens to find the optimal chunk size, and then write our own oaconvolve that uses pocketfft. IMO, scipy's oaconvolve has the following drawbacks:

• it calculates the chunk size by next_fast_len(opt_size, real=False). Can we use real=True?

• When it fallbacks, it fallbacks to fftconvolve. Can we fallback to pocketfft-based sdp?

• If it does not fallback, it performs irfftn(rfftn * rfftn) on chunks together (IIRC), not on 1D arrays. Can we use r2cn and c2rn instead?