How to vectorize this type of function?

I have a function of the form:

f(a) = sum((-1)^n/factorial(a-n))

where the summation limits are n = 0 to a. If a is a scalar then I can just substitute 0:a for n in the expression. But if a is a vector, then the limits on the sum are different for every element of a. One way to compute it would be:

for i = 1:numel(a)
    n = 0:a(i);
    f(i) = sum((-1).^n./factorial(a(i)-n));
end

However, I would like to avoid using a for loop. I am wondering if anyone has any ideas for how to vectorize this?

Matt J provided a great solution to a simpler problem HERE. Is there any way to do something similar?

[Matt J]

For starters, how big can 'a' get? Are you sure Stirling's approximation to factorial(a) won't be required?

[Oliver]

In practice, I'm using exp(gammaln(a(i)-n+1)) to avoid overflow. Just for argument's sake, lets say 'a' is never larger than 32. This is a simplified version of the real sum that I'm interested in, which is similar in form, but has 6 arguments and factorials depending on all of them (i.e. f(a,b,c,d,e,f)).

[Oliver]

Actually, more realistic is a < 10.

Answer 1

[Matt J]

   nn=0:max(a);
    T=[1,cumprod(1:nn(end))];
     T=cumsum((-1).^nn./T);

     f=T(a+1).*(-1).^mod(a,2),

[Oliver]

Hi Matt, I tried your solution and received an error on the last line because 'T(a+1)' is not the same size as 'a' so the elementwise multiplication doesn't work:

a = randi(30,1e3,1);

nn = 0:max(a);
T = [1,cumprod(1:nn(end))];
T = cumsum((-1).^nn./T);
f2 = T(a+1).*(-1).^mod(a,2);

Error using  .* 
Matrix dimensions must agree.

[Oliver]

Nevermind, it was just confusion of row vs. column vectors. Works great now. So that I can improve my coding, how did you come up with this? It seems like with this solution and the other one that you recently posted, a lot of it is based on computing the worst case scenario and using cumsum, and then cleverly indexing into the right portions of the vector of cumulative sums. Is this in general how you try to solve these types of problems?

[Matt J]

The general idea is to create a look-up table, T, of all the things you are trying to evaluate and then index into the table with look-up values like a, nzmin, nzmax. Look-up tables are things that can often be filled sequentially/incrementally, using efficient vectorized MATLAB commands like cumsum or cumprod. Indexing into these tables is, of course, also vectorizable.

The strategy might not be so optimal if your look-up values were not reasonably dense throughout the table or if the cost of creating the table was high. For that, you might make several sub-tables or fall back to a for-loop.