快速更新残差平方和答案

【问题标题】：Fast updating sum of squared residuals快速更新残差平方和
【发布时间】：2018-03-19 16:46:38
【问题描述】：

当我知道只有一小部分项发生变化时，我想找到一种快速更新残差平方和的方法。让我更详细地描述这个问题。

我有 N 个来自嘈杂阶跃函数数据的数据点。

N = 100000
realStepList = [200, 500, 900]

x = np.zeros(N)
for realStep in realStepList:
    x[realStep:] += 1
x+=np.random.randn(len(x))*0.1 #Add noise

我想计算此数据的残差平方和以及任意步骤位置列表。这是我的做法。

a = [0, 250, 550, N] 
def Q(x, a):
    q = np.sum([np.sum((x[ai:af] - i)**2) for i, (ai,af) in enumerate(zip(a[:-1],a[1:]))])
    return q

a 是我的潜在步骤列表。使用始终将0 作为第一个元素，将N 作为最后一个元素的列表更容易。

这是相对较慢的，因为它是N 平方的总和。但是，我意识到，如果我将a 更改相对较小的量，这些N 术语中的大部分将保持不变，这意味着我不必再次计算它们。

假设我已经像上面那样计算了Q(x,a)。我现在有另一个列表

b = [aa + dd for aa, dd in zip(a, d)]

其中d 是两个列表之间的差异。我不想像上面那样计算Q(x,b)（N 元素的另一个总和），我想找到

deltaQ(x, a, d) 这样

Q(x, b) = Q(x,a) + deltaQ(x, a, d)

我写过这样的函数，但是速度慢而且草率。实际上，它比Q 慢！

def deltaQ(x, a, d):
    z = np.zeros(len(x))
    J = np.zeros(len(x))
    s = 0
    for j, [dd, aa] in enumerate(zip(d, a[1:-1])):
            if dd >= 0:
                    z[aa:aa+dd] += 1
                    s += sum(x[aa:aa+dd])
            if dd < 0:
                    z[aa+dd:aa] += -1
                    s += -sum(x[aa+dd:aa])
            J[aa:] += 1
    dq = 2*s - sum((J**2 - (J-z)**2))
    return dq

这个想法是识别x 中将受到影响的所有点。例如，如果原始列表是a = [0, 5, 10] 和b = [0, 7, 10]，那么只有与x[5:7] 对应的项会在总和中发生变化。我使用列表z 跟踪这一点。然后我根据这个计算变化。

我不认为我是世界上第一个遇到这个问题的人。所以我的问题是：

有没有一种快速的方法来计算残差平方和的差异，因为这通常比从头重新计算新的总和要少得多的元素？

【问题讨论】：

次要注意，b = [aa + dd for zip(a, d)] 语法无效。
在进行这些小改动时，a 的形状是否保持不变？换句话说，只有值会改变还是步数也会改变？另外，你能举一个d的例子吗？
@Graipher 最终我想改变形状，但我对答案很感兴趣，如果它也没有改变。在我给出的示例中，d = [2]，因为只有一步，它从 5 -> 7
更改了合成器，@droooze。感谢您的提醒。
我读了几篇文章，才意识到您并没有要求找到真正的台阶位置。如果你是，我也可以帮助你:)

标签： python numpy sum

【解决方案1】：

首先，我能够使用原始代码运行Q，只修改N，以在相当标准的笔记本电脑上获得以下时间（没什么花哨的）：

N = 1e6: 0.00236s per loop
N = 1e7: 0.0260s per loop
N = 1e8: 0.251 per loop

进程在 N = 1e9 处进入交换状态，但假设您有足够的可用 RAM，我会发现 2.5 秒的时间对于该大小是完全可以接受的。

话虽如此，通过将调用 np.power 的结果中的内部 np.sum 更改为 np.ndarray.sum，我能够获得 10% 的加速：

def Q1(x, a):
    return sum(((x[ai:af] - i)**2).sum() for i, (ai, af) in enumerate(zip(a[:-1], a[1:])))

现在这是一个慢三倍的版本：

def offset(x, a):
    d = np.zeros(x.shape, dtype=np.int)
    d[a[1:-1]] = 1
    # Add out=d to make this run 4 times slower
    return np.cumsum(d)

def Q2(x, a):
    return np.sum((x - offset(x, a))**2)

为什么会有帮助？好吧，注意offset 做了什么：它将x 重新调整到您选择的基线。从长远来看，这有两件事。首先，您获得的解决方案比您目前提出的解决方案更加矢量化。其次，它允许您根据您选择的不同 b 数组重写您的 delta 函数，而不必计算 d，如果 len(a) != len(b) 甚至可能无法实现。

增量为(x - i)<sup>2</sup> - (x - i)<sup>2</sup>。如果你展开所有的混乱，你会得到(j - i)(j + i - 2x)。 j 和 i 是步骤的值，由 offset 返回。这不仅大大简化了计算，而且j - i 是您需要计算增量的掩码：

def deltaQ1(x, a, b):
    i = offset(x, a)
    j = offset(x, b)
    d = j - i
    mask = d.astype(np.bool)
    return (d[mask] * (j[mask] + i[mask] - 2 * x[mask])).sum()

此函数的运行速度比原始实现快 10 到 15 倍以上（但请记住，它需要 a 和 b 而不是 a 和 d 作为输入）。调用Q1(x, b) - Q1(x, a) 仍然快两倍。新函数还创建了一堆临时数组，但这些数组的数量很容易减少。

时间安排

除了上面显示的之外，这里是我计算机上的一些示例时序（使用提供的数据，以及a = [0, 250, 550, N]、b = [0, 180, 565, N] 和因此d = [0, -70, 15, 0]，如果相关：

原始残差：

Q:  147µs per loop
Q1: 135µs per loop <-- Use this one!
Q2: 453µs per loop

残差增量：

deltaQ: 8363µs per loop
deltaQ1: 656µs per loop
Q(x, b) - Q(x, a): 297µs per loop
Q1(x, b) - Q1(x, a): 275µs per loop  <-- Best solution?

最后说明：我有一个明显的印象，即您最初的 delta 函数实现是不正确的。它与Q(x, b) - Q(x, a) 的结果不一致，但deltaQ1(x, a, b) 一致。

TL;DR

请不要过早优化。如果你做得对，当然可以编写一个专门的 C 函数来为你在内存中保存 i - j 和 i + j，这会更快，但我怀疑你会从矢量化管道中获得很多里程。部分原因是您最终会花费大量时间来弄清楚一组复杂的索引是如何相互交织的，而不仅仅是将数字相加。

【讨论】：