How to make recursion faster in Matlab (Matlab seems to forget steps already done)?

I have defined a recursive function successfully in Matlab.

function result = f(a,i,j)
    subl = a(i:j);
    n=length(subl);
    first=subl(1);
    last=subl(end);
    if (i==j)
        result = a(i);
    else 
        result = max(first*(2*n-1)+f(a,i+1,j), last*(2*n-1)+f(a,i,j-1));
end

However, if I run the function for long enough array a. Say

a = [5 0 2 3 0 2 2 8 1 1 0 1 2 3 4 1 0 1 8 0 1 1 2 3 4]

The recursion gets really slow. This is because Matlab tends to calculate all the steps in the recursion again and again without remembering the between the steps. In Mathematica this can be easily changed by changing the recursive function a bit but I haven't found a way to do this in Matlab. Indeed, in Mathematica the same functions runs extremely fast.

Is there a way to do this easily?

[Fangjun Jiang]

"In Mathematica this can be easily changed by changing the recursive function a bit". How do you change it? I don't see anything obviously wrong. I can simplify it a bit.

[Harto Saarinen]

Well on mathematica it runs about 0.1 seconds with those initial values. And it can compute fast an array of length 200. So something must be different here.

For short arrays that both can calculate it gives exactly the same values.

[Fangjun Jiang]

+1, I find this question interesting!

[Harto Saarinen]

@Dana

There are only N(N+1) needed calculations of the function f to calculate f(a,1,N) where N is the lenght of the vector a.

For example if the array is of length 10 and we first calculate f(i,i) for all i which is trivial. Then continue using the recursion in the function it should not take 2^10 calls of the function but rather 10(10+1) calls.

[Harto Saarinen]

And in Mathematica if you define the function as

f[a_, i_, j_] := f[a, i, j] = ....

instead of just

f[a_, i_, j_] := ...

it will remember the steps.

[Stephen23]

Perhaps this might help (no idea how it works with recursive functions):

https://www.mathworks.com/help/matlab/ref/memoize.html

[Fangjun Jiang]

@Stephen Cobeldick, indeed, it can be improved.

>> tic;b=f1(a,1,20);toc
Elapsed time is 0.043371 seconds.
>> fh=@f1
fh =
  function_handle with value:
    @f1
>> memoizedFcn = memoize(fh)
memoizedFcn = 
  MemoizedFunction with properties:

     Function: @f1
      Enabled: 1
    CacheSize: 10
>> tic;b=memoizedFcn(a,1,20);toc
Elapsed time is 0.053710 seconds.
>> tic;b=memoizedFcn(a,1,20);toc
Elapsed time is 0.001782 seconds.

Answer 1

[Fangjun Jiang]

The code can be simplified first.

I see the potential saving. f(a,1,20) depends on f(a,2,20) and f(a,1,19), which depends on f(a,3,20), f(a,2,19), f(a,2,19), f(a,1,18). Here, f(a,2,19) is probabaly calculated twice. Might be able to save it.

function result = f1(a,i,j)
    n=j-i+1;
    if (i==j)
        result = a(i);
    else 
        result = max(a(i)*(2*n-1)+f1(a,i+1,j), a(j)*(2*n-1)+f1(a,i,j-1));
    end

>> a=ones(1,20);

>> tic;b=f(a,1,20);toc

Elapsed time is 0.310610 seconds.

>> tic;b=f1(a,1,20);toc

Elapsed time is 0.045661 seconds.

[Harto Saarinen]

There are only N(N+1) needed calculations of the function f to calculate f(a,1,N) where N is the lenght of the vector a. And Matlab calculates everything again and again which makes the algorithm to run in 2^N steps!

Mathematica is more clever here for some reason.

[Harto Saarinen]

If you continue your thought there you can see that Matlab calculates many of the function calls many times! Actually so many that it makes the recursion exponential even though it should not be.

Answer 2

[Fangjun Jiang]

You might not need to use recursive function to get the solution. Instead, just calculate the result as a nxn matrix, but, it needs to be in a particular order.

For a=ones(1,20), i=1, j=20, Initially, R=nan(20)

first, calculate R(1,1),R(2,2), ..., R(20,20)

second, calculate R(1,2), R(2,3), R(3,4), ... R(19,20)

third, calculate R(1,3), R(2,4), R(3,5),...

function result = f2(a,i,j)
    n=j-i+1;
    result=nan(n);
    for k=i:j
        result(k,k)=a(k);
    end
    
    for increase=1:j-i
        for k=i:j-increase
            result(k,k+increase)=max(a(k)*(2*increase+1)+result(k+1,k+increase), a(k+increase)*(2*increase+1)+result(k,k+increase-1));
        end
    end
end

[Harto Saarinen]

Yes this is something that I might have to do. It just suprises me that Matlab does not allow to compute the recursion fast enough but Mathematica is much more clever in this. Memoize does not seem to help too much.

But thank you for writing this up and for your time on the matter. The answer can be quite easily read from the upper right corner of the result now!

This is still a bit more unituitive for me. I can't so easily figure out which argument in the maximum we picked in each step now. In the recursion you can easily mark which part of the max(A,B) we took in each round so that we know at the end which side of the vector we 'consumed' at each step to achieve the maximum.

[Fangjun Jiang]

I think at issue is the fact that your recursive function called the function TWICE in one recursion. In next recursion, the function will be called 4 times and so on. I don't have knowledge why Mathematica is smarter about this. I am glad I understood the question and figured out ways to get a solution.