5

This bug is due to Matlab being too smart for its own good.

I have something like

for k=1:N
    stats = subfun(E,k,stats);
end

where statsis a 1xNarray, N=5000 say, and subfun calculates stats(k)from E, and fills it into stats

function stats = subfun(E,k,stats)
    s = mean(E);
    stats(k) = s;
end

Of course, there is some overhead in passing a large array back and forth, only to fill in one of its elements. In my case, however, the overhead is negligable, and I prefer this code instead of

for k=1:N
    s = subfun(E,k);
    stats(k) = s;
end

My preference is because I actually have a lot more assignments than just stats. Also some of the assignments are actually a good deal more complicated.

As mentioned, the overhead is negligable. But, if I do something trivial, like this inconsequential if-statement

for k=1:N
    i = k;
    if i>=1
        stats = subfun(E,i,stats);
    end
end

the assignments that take place inside subfun then suddenly takes "forever" (it increases much faster than linearly with N). And it's the assignment, not the calculation that takes forever. In fact, it is even worse than the following nonsensical subfun

function stats = subfun(E,k,stats)
    s = calculation_on_E(E);
    clear stats
    stats(k) = s;
end

which requires re-allocation of stats every time.

Does anybody have the faintest idea why this happens?

5
  • What is E? Just some constant data? Commented Sep 26, 2013 at 18:05
  • E might well be constant, or random, yes. Its size does not change. See my edit.
    – Patrick
    Commented Sep 26, 2013 at 18:06
  • 2
    The JIT of recent versions of Matlab knows not to create a new array, but to do modifications in-place. One of the requirements is that you define a function as function x = modify(x, y) and not function x_out = modify(x_in, y). Your examples seem to follow this rule, but did you check your real code? Commented Sep 26, 2013 at 18:12
  • 2
    As @Bas mentions, MATLAB has some heuristics built into it to determine whether an array can be modified in-place. I'm guessing that by adding a branch in the loop, you've confused the JIT and it no longer thinks in-place modification is possible, leading to reallocation every iteration. Given that MATLAB has no notion of reference parameters to functions you really should rewrite the function so you're calling stats(k) = subfun(E,k); within the loop. That would be better practice instead of relying on the JIT and its pattern matching (or whatever it does).
    – Praetorian
    Commented Sep 26, 2013 at 18:17
  • @Praetorian I've pondered this, but the problem is that my stats is actually a struct containing multiple arrays. I can of course have an array of structures, but I find it more awkward to get the fields from such an object.
    – Patrick
    Commented Sep 26, 2013 at 18:21

2 Answers 2

5

This might be due to some obscure detail of Matlab's JIT. The JIT of recent versions of Matlab knows not to create a new array, but to do modifications in-place in some limited cases. One of the requirements is that the function is defined as

function x = modify_big_matrix(x, i, j)
x(i, j) = 123;

and not as

function x_out = modify_big_matrix(x_in, i, j)
x_out = x_in;
x_out(i, j) = 123;

Your examples seem to follow this rule, so, as Praetorian mentioned, your if statement might prevent the JIT from recognizing that it is an in-place operation.

If you really need to speed up your algorithm, it is possible to modify arrays in-place using your own mex-functions. I have successfully used this trick to gain a factor of 4 speedup on some medium sized arrays (order 100x100x100 IIRC). This is however not recommended, could segfault Matlab if you are not careful and might stop working in future versions.

2
  • You do mean function x = modify(x, y) and not simply x = modify(x, y), right?
    – Patrick
    Commented Sep 26, 2013 at 18:31
  • Yes, my real code also has the required format. Although I would assume that it is JIT that messes up, it is truly strange that the if-statement should do it.
    – Patrick
    Commented Sep 26, 2013 at 18:31
0

As discussed by others, the problem almost certainly lies with JIT and its relatively fragile ability to modify in place.

As mentioned, I really prefer the first form of the function call and assignments, although other workable solutions have been suggested. Without relying on JIT, the only way this can be efficient (as far as I can see) is some form of passing by reference.

Therefore I made a class Stats that inherits from handle, and which contains the data array for k=1:N. It is then passed by reference.

For future reference, this seems to work very well, with good performance, and I'm currently using it as my working solution.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.