Matlab recursion: inefficient code or complex recursion?

Question

I am struggling to get the solution of this recursion problem in reasonable execution times.

Here, I show the recursive function which basically computes the coefficients of a polynomial.

function [ coeff ] = get_coeff( n, k, tau, x )

if(n == 0) % 1st exit condition
    coeff = 0;
else
    if(k == 0) % 2nd exit condition
        coeff = max(0, n*tau-x)^n;
    else % Else recursion
        total = 0;
        for l = k-1:n-2
            total = total + nchoosek(l, k-1)*tau^(l-k+1)*get_coeff(n-1, l, tau, x);
        end
        coeff = (n/k) * total;            
    end
end

end

 % This symbolic summation solution gives numerical errors, probably due to rounding
 % effects.
 %           syms l;
 %           f = nchoosek(l, k-1)*tau^(l-k+1)*get_coeff(n-1, l, tau, x);
 %           coeff = (n/k) * symsum(f, l, k-1, n-2);

And this is the main script where I make use of the recursive function:

Tau = 1;
ns = [3];
%delays = 0:0.25:8;
delays = [0];
F_x = zeros(1, size(delays, 2));
rho = 0.95;
tic
for ns_index = 1: size(ns, 2)

  T = Tau*(ns(ns_index)+1)/rho;

  % Iterate delays (x)
  for delay_index = 1:size(delays, 2)
     total = 0;

     % Iterate polynomial.
     for l = 0:ns(ns_index)-1
        total = total + get_coeff(ns(ns_index), l, Tau, delays(delay_index))*(T - ns(ns_index)*Tau + delays(delay_index))^l;
     end

    F_x(1, delay_index) = T^(-ns(ns_index))*total;

  end

end
toc

I've simplified, "ns" and "delays" vectors to contain a single value so that it is easier to follow. In summary, for a fixed value of "ns", I need to compute all the coefficients of the polynomial using the recursive function and compute its final value at "delays". By increasing the number of points in "delays", I can see a curve for a fixed "ns". My question is: for any "ns" between 1 and 10, the computation is really fast, in the order of 0.069356 seconds (even for the whole "delays" vector). Conversely, for ns = [15] or [20], the computation time increases A LOT (I didn't even manage to see the result). I'm not keen on assessing computational complexity, so I don't know if there is a problem in my code (maybe nchoosek function?, or for loops?) or maybe it is the way it has to be having in mind this recursion problem.

EDIT: I see it is indeed the factorial growth of the amount of calculations, as Adriaan stated. Do you think that any kind of approximation of nchoosek could be useful to tackle this problem? Something like: en.wikipedia.org/wiki/Stirling%27s_approximation

The last formula in this paper is what I'm trying to implement (note I changed delta for tau):

Answer 1

So I've finally managed to compute the coefficients in a reasonable amount of time. Basically, I took the suggestions from Adriaan and rahnema1 and created a ns by ns matrix to store all the coefficients that I compute in a recursive way. Therefore, when a certain leaf of the recursive tree is repeated I am able to prune the tree by extracting the value from the matrix. Note that the gain is not based on precomputing the values (since I compute them on the go) but on pruning the number of recursions. Here you have some numbers:

• ns = 10 ; delay = 0: the number of calls to the old recursive function was 23713. Now, this is solved in 175 calls.
• For ns = 10; delay = [0:0.25:8]: 782529 calls with the old function and 2.74 seconds of execution time, 495 to the new one and 0.02 which is ~ 125x times faster.