Matlab large array minus

Question

I am trying to optimize a Matlab code for a statistic calculation for a large array of data (1e6 values). I tried several methods, with loops or fun functions, with diff or basic math. For me this code runs in approx 8 seconds. Though, using an additional for loop runs fine 'run section' mode, but it is unstable while running the script.

Is there any way to improve the time for this code. I tried different methods for parfor, but i could not set it up in my loop. Is anyone here experienced enough with parfor to tell me how do I get out of broadcast variable issues?

%% Matlab Question
L=400000;
t_clk = rand(1, L);
t_clk = t_clk-0.5;
plot (t_clk)
%
disp(' ')
tic
N = 1000; %2000
M = length(t_clk)-N;
temp_Pkp = zeros(1, N);
temp_Pkn = zeros(1, N);
temp_Std = zeros(1, N);
myMat = zeros(1,M);
% Time to execute: 'run section' / 'run' / 'run and time'
%parfor xx = 1 :1  : N  %2.3 -> broadcast variable issues
for xx = 1 :1  : N  %2.3

    myMat =  bsxfun(@minus,t_clk(xx+1 : xx+M) , t_clk(1:M)); %% Time to execute: 8.1s / 8.2s / 7.76s
    %myMat =  t_clk(xx+1 : xx+M) - t_clk(1:M); %% Time to execute: 8.1s / 8s / 86s
    %myMat = zeros(1,M); % not used
    %for yy = 1:1:M %% Time to execute: 4.6s / 4.6s / 23.6s ('run and time' execution time is very high)
    %    myMat(yy) =t_clk(xx+yy) - t_clk(yy);
    %end
    temp_Mean= mean(myMat) ;
    temp_Pkp(xx) = max(myMat(:)) - temp_Mean  ; % max - min
    temp_Pkn(xx) = temp_Mean  - min(myMat(:))  ; % max - min
    temp_Std(xx) = sqrt(sum(sum((myMat-temp_Mean).^2))/M);
end
toc
plot(temp_Std)

It would help if you described what the code is supposed to *do*. "Here is code, please make it faster" is not a well formed question. — Wolfie, Sep 10 '18 at 14:03
`parfor` is not a magic wand. Read [this question](https://stackoverflow.com/q/32146555/5211833) and answers thereon on its intricacies. Start by vectorising first, and once a code is fully optimised, check whether parallellisation can help you furter. — Adriaan, Sep 10 '18 at 14:13
Hello @Wolfie. I wrote in the beginning that it does some statistical calculation on a large set of data. In this case, vector t_clk is a zero crossing vector of a clock signal. By calculating myMat in the way I do, I can compute basic jitter values on this signal: RMS, max and min. — uhnucross, Sep 10 '18 at 20:29
"Some statistical calculation" is not a problem statement. Give us either the mathematical definition of what you're trying to do, or the desired output for a given input (using `rand` like this makes things ambiguous and non reproducible), or both! Describing things with vague words like "basic jitter" doesn't really help, it might be your whole problem can be done in 2 lines, but we'd have to guess whether a given suggestion is correct based on not a lot. Please read [ask] and provide a [mcve]. This would really help us to help you! — Wolfie, Sep 10 '18 at 21:16
Clock Jitter is defined as period deviation from an ideal clock. An ideal clock does not exist in real life, so I will have to consider my first clock period as reference (t_clk(1)). If I want to know the period jitter of my clock, I will have to make a difference between all clock edges and #1 then use hist function in matlab and get the statistics on it. Though, this is not enough. Sometimes a deviation after N periods is needed. Therefore an accumulated jitter must be calculated. This means that I decimate my clock N times and do the same calculations as before. — uhnucross, Sep 10 '18 at 21:41
In my case I would like to know the statistical data (min, max, rms/std) from N = 1 to 1000(or 2000) for L=400000 clock periods acquired.It would be enough to shift N by N (xx+1:N:xx+M, etc) but this would mean to loose some important information. Therefore my shift is always 1. This is increasing calculation time. Calculation is performed multiple times, so single execution time matters. — uhnucross, Sep 10 '18 at 21:46
I am using a rand of 400000 values in my tests to optimise the code. In reality I have 400000 period values of a clock, format longEng — uhnucross, Sep 10 '18 at 21:59
Please [edit] your question, comments aren't for adding new information critical to any answer. — Wolfie, Sep 11 '18 at 06:48
(1) `bsxfun` does nothing there, except cost time. (2) When timing code, put it inside a function, and time the function call. Functions are optimized better by the JIT, and therefore typically run faster. Ignore "run and time" stuff, what matters is how quickly the function runs. Use the profiler to find hotspots. (3) Optimize by collecting `mean(myMat)`, `max` and `min` in arrays, then computing the other values outside your loop in a vectorized manner. (4) Use `std`, don't roll your own. — Cris Luengo, Sep 11 '18 at 16:01
Hello @CrisLuengo. Thanks for your answer. 1. bsxfun minus is actually a bit faster in this case that using '-', especially for large N(=2000) 2. I will try to do so 3. I tried this once and it wasn't faster. This will generate myMat(1:1000,1:40000) and the time to calculate statistics on this exploded 4. Calculating std in this way is a bit faster than Matlab's std() — uhnucross, Sep 12 '18 at 11:12
Nevermind (3), I mis-read your code. But on (1): `bsxfun` in this case calls the function handle exactly once, because there are no singleton dimensions to expand. Try comparing `timeit(@()bsxfun(@minus,a,b))` vs `timeit(@()a-b)` with `a=rand(1,L)` and `b=rand(1,L)`. I see 3.3120e-04 in the first case, and 2.2348e-04 in the second. The difference is the overhead of the call to `bsxfun`. So `bsxfun` has a little bit of overhead, and it makes your code harder to read. You should not use it. And always use `timeit` for timing your code. It reduces variability by executing the function many times. — Cris Luengo, Sep 12 '18 at 14:26
Regarding (4): yes, `std` will re-compute the mean, which you already have. But you’re not computing the standard deviation, you need to divide by `M-1` for that. — Cris Luengo, Sep 12 '18 at 14:39

score 1 · Answer 1 · answered Sep 11 '18 at 15:14

I did some fiddling and it appears you're indexing a variable repeatedly and it's costing you a lot. This variable is t_clk and you're indexing by 1:M every time. If you create a temporary variable with this indexing you can drastically speed up the runtime.

%% Matlab Question
L=400000;
t_clk = rand(1, L);
t_clk = t_clk-0.5;
plot (t_clk)
%
disp(' ')
tic
N = 1000; %2000
M = length(t_clk)-N;
temp_Pkp = zeros(1, N);
temp_Pkn = zeros(1, N);
temp_Std = zeros(1, N);
myMat = zeros(1,M);
% Time to execute: 'run section' / 'run' / 'run and time'
%parfor xx = 1 :1  : N  %2.3 -> broadcast variable issues
t_clk_t = t_clk(1:M);
idx = 1:M;
for xx = 1 :1  : N  %2.3

    myMat =  bsxfun(@minus,t_clk(xx+1 : xx+M) , t_clk(1:M)); %% Time to execute: 3.043

    myMat =  bsxfun(@minus,t_clk(xx+1 : xx+M) , t_clk_t); % Time to execute 1.740

%     myMat = t_clk(xx+1 : xx+M) - t_clk(1:M);
    %myMat =  t_clk(xx+1 : xx+M) - t_clk(1:M); %% Time to execute: 8.1s / 8s / 86s
    %myMat = zeros(1,M); % not used
    %for yy = 1:1:M %% Time to execute: 4.6s / 4.6s / 23.6s ('run and time' execution time is very high)
    %    myMat(yy) =t_clk(xx+yy) - t_clk(yy);
    %end
    temp_Mean= mean(myMat) ;
    temp_Pkp(xx) = max(myMat(:)) - temp_Mean  ; % max - min
    temp_Pkn(xx) = temp_Mean  - min(myMat(:))  ; % max - min
    temp_Std(xx) = sqrt(sum(sum((myMat-temp_Mean).^2))/M);
end
toc
plot(temp_Std)

Hello @Durkee, Thanks for answer. I was working on this a bit yesterday and I noticed something similar. I had additional timing issues in my main code because I was working with structures. Plain variables + index skipping helps quite a lot. I have another idea that I will have to try: I noticed that dividing L/4.../10 time is is decreasing more than this division factor. I will have to use reshape + some additional stuff to optimize the matrix based on used indexes. Will try to work on this this week and will put here a code once is ready. — uhnucross, Sep 12 '18 at 11:17

Matlab large array minus

1 Answers1