Стабильно вычислять большие количества через рекурсию

Question

У меня есть две величины a & b, которые определены рекурсией и ссылкой на другой список значений x = [ x_1, x_2, ... x_N ], который будет входом в программу. Программа будет перебирать все значения в x и обновлять a & b в соответствии с:

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])  

и начальные значения

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

Значения в x не являются большими числами (всегда <10), но N будут в сотнях, и из-за всех экспонент конечные значения a & b будут очень большими. Я обеспокоен тем, что из-за формы рекурсии, когда мы постоянно умножаем экспоненциально большие числа на экспоненциально малые и добавляем их, эта схема станет довольно численно нестабильной.

В идеале я бы рассчитал log (a) и log (b), чтобы значения не становились слишком большими. Но из-за схемы рекурсии это невозможно, если я не вычислю что-то более запутанное

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

Является ли численная стабильность чем-то, что меня беспокоит? И если так, то что-то вроде схемы, основанной на логах, поможет стабилизировать ее?

Answer 1

Мы можем переписать это сначала как:

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]))

Затем измените наши итерационные переменные:

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]))

Который может быть вычислен более стабильно, используя 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])

Реализация logaddexp можно увидеть здесь

Answer 2

Насколько я знаю, все (?) Рекурсивные проблемы могут быть решены с помощью динамического программирования. Например, последовательность Фибоначчи может быть выражена так:

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

Или, итеративно:

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

Предположительно, если у вас есть рекурсивная проблема, вы можете выполнить аналогичную распаковку.