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

Question

atharva aalok 2023-8-19

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2010327-skewness-calculator-using-one-pass-algorithm-code-speed-up-needed

编辑： Bruno Luong 2023-8-20

采纳的回答： Bruno Luong

在 MATLAB Online 中打开

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!

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

atharva aalok 2023-8-19

编辑：atharva aalok 2023-8-19

在 MATLAB Online 中打开

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 2023-8-19

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2010327-skewness-calculator-using-one-pass-algorithm-code-speed-up-needed#answer_1289862

编辑：Bruno Luong 2023-8-20

在 MATLAB Online 中打开

skewness_mex_mexw64.mat

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

Answer 2

Bruno Luong 2023-8-19

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2010327-skewness-calculator-using-one-pass-algorithm-code-speed-up-needed#answer_1289692

编辑：Bruno Luong 2023-8-19

在 MATLAB Online 中打开

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

3 个评论
显示 1更早的评论隐藏 1更早的评论

Bruno Luong 2023-8-19

编辑：Bruno Luong 2023-8-19

在 MATLAB Online 中打开

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 2023-8-19

编辑：Bruno Luong 2023-8-19

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

请先登录，再进行评论。

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

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

采纳的回答

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

更多回答（1 个）

3 个评论
显示 1更早的评论隐藏 1更早的评论

另请参阅

类别

标签

Community Treasure Hunt

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

1 个评论 显示 -1更早的评论隐藏 -1更早的评论

采纳的回答

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

更多回答（1 个）

3 个评论 显示 1更早的评论隐藏 1更早的评论

另请参阅

类别

标签

Community Treasure Hunt

WeChat

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

3 个评论
显示 1更早的评论隐藏 1更早的评论