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

Question

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

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

Значения в х не большие числа (всегда < 10), но N будет в сотнях, и из-за всех экспонент конечные значения & 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 можно увидеть here

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

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