通过递归稳定计算大数量

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的实现可见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])

大概如果您有递归问题，您可以执行类似的解包。

通过递归稳定计算大数量

Stably computing large quantites through recursion

python

numpy

scipy

numerical-stability