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

7 visualizzazioni (ultimi 30 giorni)
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 Commento
atharva aalok
atharva aalok il 19 Ago 2023
Modificato: atharva aalok il 19 Ago 2023
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

Accedi per commentare.

Risposta accettata

Bruno Luong
Bruno Luong il 19 Ago 2023
Modificato: Bruno Luong il 20 Ago 2023
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

Più risposte (1)

Bruno Luong
Bruno Luong il 19 Ago 2023
Modificato: Bruno Luong il 19 Ago 2023
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 Commenti
Bruno Luong
Bruno Luong il 19 Ago 2023
Modificato: Bruno Luong il 19 Ago 2023
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
Bruno Luong il 19 Ago 2023
Modificato: Bruno Luong il 19 Ago 2023
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

Accedi per commentare.

Categorie

Scopri di più su Programming in Help Center e File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by