压缩折扣率计算

Question

金融学和强化学习中的一个常见术语是基于原始报酬C[i]时间序列的贴现累积报酬R[i]。给定一个数组R，我们希望用C[-1] = R[-1]计算满足递归C[i] = R[i] + discount * C[i+1]的C[i] (并返回完整的数组C)。

在使用numpy数组的python中，一种数值稳定的计算方法可能是：

import numpy as np
def cumulative_discount(rewards, discount):
    future_cumulative_reward = 0
    assert np.issubdtype(rewards.dtype, np.floating), rewards.dtype
    cumulative_rewards = np.empty_like(rewards)
    for i in range(len(rewards) - 1, -1, -1):
        cumulative_rewards[i] = rewards[i] + discount * future_cumulative_reward
        future_cumulative_reward = cumulative_rewards[i]
    return cumulative_rewards

但是，这依赖于python循环。考虑到这是一种很常见的计算方法，当然有一个现有的矢量化解决方案，它依赖于其他一些标准函数，而不需要进行细胞化。

注意，任何使用np.power(discount, np.arange(len(rewards))之类的解决方案都是不稳定的。

Answer 1

您可以使用scipy.signal.lfilter来解决递归关系：

def alt(rewards, discount):
    """
    C[i] = R[i] + discount * C[i+1]
    signal.lfilter(b, a, x, axis=-1, zi=None)
    a[0]*y[n] = b[0]*x[n] + b[1]*x[n-1] + ... + b[M]*x[n-M]
                          - a[1]*y[n-1] - ... - a[N]*y[n-N]
    """
    r = rewards[::-1]
    a = [1, -discount]
    b = [1]
    y = signal.lfilter(b, a, x=r)
    return y[::-1]

此脚本测试结果是否相同：

import scipy.signal as signal
import numpy as np

def orig(rewards, discount):
    future_cumulative_reward = 0
    cumulative_rewards = np.empty_like(rewards, dtype=np.float64)
    for i in range(len(rewards) - 1, -1, -1):
        cumulative_rewards[i] = rewards[i] + discount * future_cumulative_reward
        future_cumulative_reward = cumulative_rewards[i]
    return cumulative_rewards

def alt(rewards, discount):
    """
    C[i] = R[i] + discount * C[i+1]
    signal.lfilter(b, a, x, axis=-1, zi=None)
    a[0]*y[n] = b[0]*x[n] + b[1]*x[n-1] + ... + b[M]*x[n-M]
                          - a[1]*y[n-1] - ... - a[N]*y[n-N]
    """
    r = rewards[::-1]
    a = [1, -discount]
    b = [1]
    y = signal.lfilter(b, a, x=r)
    return y[::-1]

# test that the result is the same
np.random.seed(2017)

for i in range(100):
    rewards = np.random.random(10000)
    discount = 1.01
    expected = orig(rewards, discount)
    result = alt(rewards, discount)
    if not np.allclose(expected,result):
        print('FAIL: {}({}, {})'.format('alt', rewards, discount))
        break

Answer 2

你所描述的计算被称为霍纳法则或Horner的多项式计算方法。它是用NumPy polynomial.polyval实现的。

但是，您需要整个cumulative_rewards列表，即Horner规则的所有中间步骤。NumPy方法不返回这些中间值。您的函数，用Numba的@jit装饰，可能是最佳的选择。

作为理论上的可能性，我将指出，如果给定一个系数的polyval，Hankel矩阵可以返回整个列表。这是矢量化的，但最终效率低于Python循环，因为cumulative_reward的每个值都是从头开始计算的，独立于其他值。

from numpy.polynomial.polynomial import polyval
from scipy.linalg import hankel

rewards = np.random.uniform(10, 100, size=(100,))
discount = 0.9
print(polyval(discount, hankel(rewards)))

的输出与

print(cumulative_discount(rewards, discount))

Answer 3

如果您想要一种只适用于numpy的解决方案，那么可以这样做(从unutbu的答案中借用结构)：

def alt2(rewards, discount):
    tmp = np.arange(rewards.size)
    tmp = tmp - tmp[:, np.newaxis]
    w = np.triu(discount ** tmp.clip(min=0)).T
    return (rewards.reshape(-1, 1) * w).sum(axis=0)

下面的证据。

import numpy as np

def orig(rewards, discount):
    future_cumulative_reward = 0
    cumulative_rewards = np.empty_like(rewards, dtype=np.float64)
    for i in range(len(rewards) - 1, -1, -1):
        cumulative_rewards[i] = rewards[i] + discount * future_cumulative_reward
        future_cumulative_reward = cumulative_rewards[i]
    return cumulative_rewards

def alt2(rewards, discount):
    tmp = np.arange(rewards.size)
    tmp = tmp - tmp[:, np.newaxis]
    w = np.triu(discount ** tmp.clip(min=0)).T
    return (rewards.reshape(-1, 1) * w).sum(axis=0)

# test that the result is the same
np.random.seed(2017)

for i in range(100):
    rewards = np.random.random(100)
    discount = 1.01
    expected = orig(rewards, discount)
    result = alt2(rewards, discount)
    if not np.allclose(expected,result):
        print('FAIL: {}({}, {})'.format('alt', rewards, discount))
        break
else:
    print('success')

然而，这个解决方案不能很好地扩展到大奖励数组，但是您仍然可以使用大步技巧，正如这里所指出的。