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
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 个评论
采纳的回答
Bruno Luong
2023-8-19
编辑：Bruno Luong
2023-8-20
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
更多回答（1 个）
Bruno Luong
2023-8-19
编辑：Bruno Luong
2023-8-19
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 个评论
atharva aalok
2023-8-19
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_Bruno = timeit(@()skewness_onepass_BLU2(timeseries),1) % 0.0392 second
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
类别
在 帮助中心 和 File Exchange 中查找有关 Write C Functions Callable from MATLAB (MEX Files) 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
