Skewness calculator using one pass algorithm. Code speed up needed.

Hello,

I am using a one pass algorithm to calculate skewness. My code is as follows:

n = 10000000;
timeseries = randi(1000, 1, n);

tic
sk = skewness_onepass(timeseries);
toc

Elapsed time is 0.075032 seconds.

function skewness = skewness_onepass(x)
    
    N = length(x);
    M1 = 0;
    M2 = 0;
    M3 = 0;
    
    % Algorithm taken from here:
    % https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Higher-order_statistics
    for n = 1: N
        delta = x(n) - M1;
        M1 = M1 + delta / n;
        M3 = M3 + delta * delta * delta * (n - 1) * (n - 2) / (n * n) - 3 * delta * M2 / n;
        M2 = M2 + delta * delta * (n - 1) / n;
    end
    
    std_dev = sqrt(M2 / N);

    % Calculate skewness
    skewness = (M3 / N) / (std_dev * std_dev * std_dev);
    
end

My Concern:

Is there anything I can do to speed up this code? My application is really time critical and I want to evaluate skewness as fast as possible.

Reading online seems to suggest that extending MATLAB using C++ (I know bit of C but not C++) could help, but I have no clue about this. If this is indeed something that can help, please provide code/snippet on how to execute this.

I also have other functions like this and I am hoping that lessons learnt on this can be used on those as well.

Thanks!

[atharva aalok]

Comparing performance of difference suggested answers:

n = 10000000;

time_my = 0;
time_bruno = 0;

% My function
for i = 1: 100
    timeseries = randi(1000, 1, n);
    
    tic_my = tic;
    sk = skewness_onepass(timeseries);
    val = toc(tic_my);
    time_my = time_my + val;
    
    % Bruno's function
    tic_bruno = tic;
    sk = skewness_onepass_bruno(timeseries);
    val = toc(tic_bruno);
    time_bruno = time_bruno + val;
end

fprintf('My average time = %f\n', time_my / 100);

My average time = 0.064831

fprintf('Bruno average time = %f\n', time_bruno / 100);

Bruno average time = 0.066629

Answer 1

[Bruno Luong]

Implement skewness in C-mex

#include "mex.h"
#include "matrix.h"
#include "math.h"

/* 
compile command
mex -O -R2018a skewness_mex.c 
*/

#define X prhs[0]
void mexFunction(int nlhs, mxArray *plhs[],
                 int nrhs, const mxArray *prhs[]) {

    double *x, delta, a, c, M1, M2, M3, skewness;
    size_t n, N;

    x = mxGetDoubles(X);
    N = mxGetN(X);
    M1 = M2 = M3 = 0;
    for (n=0; n<N; n++)
    {
        delta = x[n] - M1;
        a = delta / (n+1);
        M1 += a;
        c = delta * a * n;
        M3 += a * (c * (n - 1) - 3 * M2);
        M2 += c;
    }
    skewness = M3 / (M2 * sqrt(M2/N));   

    plhs[0] = mxCreateDoubleScalar(skewness);


} /* skewness_mex.c */

For both MSVS and Intel OneAPI compilers, as I have anticipated it is slower than MATLAB equivalent code, see timings in the code comments.

It shows that MATLAB EE is getting very good.

Your friend in Stackoverflow can come here and see the result.

For those who are on Windows and want to try the mex, pleas rename the attached mat file with .mexw64 extension. I do it to trick Answer.

n = 10000000;
timeseries = randi(1000, 1, n);

t_org = timeit(@()skewness_onepass(timeseries),1)           % 0.0432 second
t_Bruno = timeit(@()skewness_onepass_BLU2(timeseries),1)    % 0.0399 second
t_mex = timeit(@()skewness_mex(timeseries),1)               % 0.0414 second

% skewness_onepass(timeseries)
% skewness_onepass_BLU2(timeseries)
% skewness_mex(timeseries)

function skewness = skewness_onepass(x)

N = length(x);
M1 = 0;
M2 = 0;
M3 = 0;

% Algorithm taken from here:
% https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Higher-order_statistics
for n = 1: N
    delta = x(n) - M1;
    M1 = M1 + delta / n;
    M3 = M3 + delta * delta * delta * (n - 1) * (n - 2) / (n * n) - 3 * delta * M2 / n;
    M2 = M2 + delta * delta * (n - 1) / n;
end

std_dev = sqrt(M2 / N);
% Calculate skewness
skewness = (M3 / N) / (std_dev * std_dev * std_dev);

end

function skewness = skewness_onepass_BLU2(x)   
    N = length(x);
    M1 = 0;
    M2 = 0;
    M3 = 0;
    for n = 1: N
        delta = x(n) - M1;
        a = delta / n;
        M1 = M1 + a;
        c = delta * a * (n - 1);
        M3 = M3 + a * (c * (n - 2) - 3 * M2);
        M2 = M2 + c;
    end
    skewness = M3 / (M2 * sqrt(M2/N));    
end

Answer 2

[Bruno Luong]

Simple factorize common quantities used in updating M1, M2, M3

function skewness = skewness_onepass_BLU2(x)   
    N = length(x);
    M1 = 0;
    M2 = 0;
    M3 = 0;
    for n = 1: N
        delta = x(n) - M1;
        a = delta / n;
        M1 = M1 + a;
        c = delta * a * (n - 1);
        M3 = M3 + a * (c * (n - 2) - 3 * M2);
        M2 = M2 + c;
    end
    skewness = M3 / (M2 * sqrt(M2/N));    
end

EDIT: some more simplification

[atharva aalok]

@Bruno Luong This seems to be slower than the one that I suggested. Please see comment on question.

[Bruno Luong]

On my laptop my edited code is slightly faster, not much though. Up to you to stay with your code. On the flops point of view mine is mess demanding.

n = 10000000;
timeseries = randi(1000, 1, n);

t_org = timeit(@()skewness_onepass(timeseries),1) % 0.0429 second

t_org = 0.0651

t_Bruno = timeit(@()skewness_onepass_BLU2(timeseries),1) % 0.0392 second

t_Bruno = 0.0666

function skewness = skewness_onepass(x)

N = length(x);
M1 = 0;
M2 = 0;
M3 = 0;

% Algorithm taken from here:
% https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Higher-order_statistics
for n = 1: N
    delta = x(n) - M1;
    M1 = M1 + delta / n;
    M3 = M3 + delta * delta * delta * (n - 1) * (n - 2) / (n * n) - 3 * delta * M2 / n;
    M2 = M2 + delta * delta * (n - 1) / n;
end

std_dev = sqrt(M2 / N);
% Calculate skewness
skewness = (M3 / N) / (std_dev * std_dev * std_dev);

end

function skewness = skewness_onepass_BLU2(x)   
    N = length(x);
    M1 = 0;
    M2 = 0;
    M3 = 0;
    for n = 1: N
        delta = x(n) - M1;
        a = delta / n;
        M1 = M1 + a;
        c = delta * a * (n - 1);
        M3 = M3 + a * (c * (n - 2) - 3 * M2);
        M2 = M2 + c;
    end
    skewness = M3 / (M2 * sqrt(M2/N));    
end

[Bruno Luong]

Here is the comparison of operations for each loop iteration.

original code

• +-: 8
• *: 9
• /: 4
• =: 4
• indexing: 1

My code

• +-: 7
• *: 5
• /: 1
• =: 6
• indexing: 1

So my code reduces

• 1 (+-) operation;
• 4 (*) operation;
• 3 (/) operation;

but increases

• 2 affectations