Stably computing large quantites through recursion

Question

I have two quantities a & b that are defined by recursion and through reference to another list of values x = [ x_1, x_2, ... x_N ], which will be an input to the program. The program will iterate over all the values in x and update a & b according to:

for n in range(1,N)
    a[n] = a[n-1] * exp(+x[n]) + b[n-1] * exp(-x[n])  
    b[n] = b[n-1] * exp(+x[n]) + a[n-1] * exp(-x[n])

and starting values

a[0] = exp(+x[0])
b[0] = exp(-x[0])

The values in x are not big numbers (always <10) but N will be in the hundreds, and because of all the exponentials the final values of a & b will be very large. I'm worried that because of the form of the recursion where we are constantly multiplying exponentially large numbers with exponentially small ones and adding them this scheme will become quite numerically unstable.

Ideally I would calculate log(a) and log(b) instead to stop the values becoming too large. But because of the recursion scheme that's not possible, unless I compute something much messier like

log_a[n] = x[n] + log_a[n-1] + log( 1 + exp(-2*x[n] + log_b[n-1]-log_a[n-1] ) )

Is numerical stability something I am right to be concerned about here? And if so would something like the log based scheme help to stabilise it?

Your messier expression can be improved by using `log1p(x) == log(1+x)`, which is more numerically stable — Eric, Feb 05 '16 at 19:34
I'm assuming your starting values use `x[0]`, and not `x[n]`? — Eric, Feb 05 '16 at 19:45
Yes, that's right. It should be x[0] in the starting conditions. Thanks for the tip about log1p too — RGWinston, Feb 05 '16 at 19:46

Eric · Accepted Answer · 2016-02-05T19:44:30.300

1

We can rewrite that first as:

for n in range(1,N)
    a[n] = exp(log(a[n-1]) + x[n]) + exp(log(b[n-1]) - x[n])
    b[n] = exp(log(b[n-1]) + x[n]) + exp(log(a[n-1]) - x[n]))

Then change our iteration variables:

for n in range(1,N)
    log_a[n] = log(exp(log_a[n-1] + x[n]) + exp(log_b[n-1] - x[n]))
    log_b[n] = log(exp(log_b[n-1] + x[n]) + exp(log_a[n-1] - x[n]))

Which can be computed more stably using np.logaddexp:

for n in range(1,N)
    log_a[n] = np.logaddexp(log_a[n-1] + x[n], log_b[n-1] - x[n])
    log_b[n] = np.logaddexp(log_b[n-1] + x[n], log_a[n-1] - x[n])

The implementation of logaddexp can be seen here

edited Feb 05 '16 at 19:44

answered Feb 05 '16 at 19:38

Eric

95,302
53
242
374

This is great, thanks! Are there any disadvantages at all to using logaddexp, and do you think numerical stability would have been a problem with the original code? – RGWinston Feb 05 '16 at 19:52
@RGWinston: I'm afraid I don't know the answer to either of those questions. – Eric Feb 06 '16 at 00:29

score -1 · Answer 2 · answered Feb 05 '16 at 19:41

-1

As far as I'm aware, all(?) recursive problems can be solved through dynamic programming. For example, the Fibonacci sequence could be expressed like so:

def fibo(n):
    if n == 0:
        return 0
    elif n == 1:
        return 1
    return fibo(n-1) + fibo(n-2)

Or, iteratively:

n = 10
fibo_nums = [0, 1]
while len(fibo_nums) <= n:
    fibo_nums.append(fibo_nums[-2] + fibo_nums[-1])

Presumably if you have a recursive problem you could perform a similar unpacking.

answered Feb 05 '16 at 19:41

Wayne Werner

49,299
29
200
290

The OP has already performed such an unpacking. The formula itself is "recursive", but the code in the question is not – Eric Feb 05 '16 at 19:42
@Eric ah. I didn't see any recursion in the OP's code, and am unfamiliar with the problem domain. And since the question title talks about recursion I assumed there was something happening elsewhere in their code. – Wayne Werner Feb 05 '16 at 19:54

Stably computing large quantites through recursion

2 Answers2