Re: efficient calculatin method of finite product series

On 16 Sep., 10:22, angelo <nicol...@xxxxxxxxxxxxxx>
wrote:
Does anyone have an efficient computational method,
preferably recursive, to calculate the following
finite product series without compromising accuracy,

\Prod_{j=1}^{k-1} (1-exp(-(t_k-t_{k-j}))

where the number of t_k's increases over time of
the simulation,
t_k's are real positive numbers, and t_k > t_{k-1}
(where k is an integer greater than 2).
Each time a new t_{k+1} is generated then the above
finite product needs
to be re-calculated. One could use a brute force
method but this is computational inefficient since
you have an ever increasing number of exponential to
evaluate.

In suggestions or solutions will be greatly
appreciated

Saving the values u_k := exp(-t_k) and using
1 - exp(-(t_k - t_{k-j})
= 1 - u_k/u_{k-j}
= (u_{k-j} - u_k)/u_{k-j}
suggests itself.

Defining
P[k] := \Prod_{j=1}^{k-1} (1-exp(-(t_k-t_{k-j}))
for your desired product and
Q[k] := \Prod_{j=1}^{k-1} u_j
we now have
P[k] = 1/Q[k] * \Prod_{j=1}^{k-1} (u_{k-j} - u_k)
and a simple recursion
Q[k+1] = u_k * Q[k]
The step from k to k+1 is of course still O(k)
(and unavoidably so), but it costs only
- one evaluation of exp (instead of k-1)
- k-1 differences
- k+1 multiplications
- one division
Admittedly, this does not save very much yet.

Since the t_k are increasing, the u_k are decreasing
towards 0.
During the calculation of
\Prod_{j=1}^{k-1} (u_{k-j} - u_k)
you observe for some j that u_k << u_{k-j}.
This allows you to replace the remaining factors
u_{k-j}-u_k by u{k-
j}, and this product has already been saved in
Q[k-j].
(IOW, you let
P[k] = Q[m]/Q[k] * \Prod_{j=1}^{k-m} (u_{k-j} - u_k)

The relative error is at most
1-(1-u_k/u_{k-j})^{k-j} or approximately
(k-j)*u_k/u_{k-j}.
This seems to suggest that you have to save all
values of Q[i], but in
fact you can throw away all older values of Q[i] (and
also of u_i)
once u_k has become small enough.

Additional precision can be obtained (back) by
keeping
track of the derivative of
f_k(x) := \Prod_{j=1}^{k-1} (u_{k-j} - x)
at x=0:
f'_{k+1}(x) = f'_k(x) (u_k - x) - f_k(x),
hence (with R[k] := -f'_k(0))
R[k+1] = R[k] * u_k + Q[k].
Now we may simplify to
P[k] = (Q[m] - u_k*R[m])/Q[k]
* \Prod_{j=1}^{k-m} (u_{k-j} - u_k)
where m should be chosen such that u_k*R[m]<<Q[m].

Note that the usefulness of these two suggestions
depends on how fast
the t_k are growing.


Dear Hagman,

Thank you for your prompt suggestion, it was quite comprehensive. I'm impressed! It seems that one is
still forced to save some subset of u_k in an array,
which may not be such a bad compromise. One thing is
for sure .. one does not want to save all u_k's into
an array where the size increases. Ideally, if there
was no need to keep some of the u_k in an array, then
the algorithm would be extremely fast.

Anyway, thanks again for your suggestions, I am sure that they will come in handy.

If anyone else would like to add their own suggestion(s) then I would be most grateful.
.

Quantcast