双栈滑动窗口聚合

聚合是对一组数据的某种汇总，可以是求和、长度、最小值等等。我们经常需要反复计算这种汇总，比如为噪声检测器计算“过去 30 秒内的最大噪声分贝值”。这种情况下，我们说数据上有一个滑动窗口，而我们想在窗口内做聚合。

如果聚合操作是一个带逆元的二元运算符，比如整数加法，那么用双端队列就能非常简单地解决：

from collections import deque

class SlidingWindowSum:
 def __init__(self):
 self.sum = 0
 self.elems = deque()

 def push(self, x):
 self.sum += x
 self.elems.append(x)

 def pop(self):
 self.sum -= self.elems.popleft()

 def eval(self):
 return self.sum

但如果运算符没有逆元呢？ 大多数有意思的汇总都是这种情况，比如最小值、分位数、近似去重计数（例如用 HyperLogLog）。事实上，就连浮点数求和这样简单的操作也受此困扰，因为浮点加法不可逆。比如，只要输入数据里出现过一个 NaN，用上面那个朴素算法算出的和就会永远停在 NaN，哪怕那个坏值早已滑出窗口。

六年前，我为滑动窗口中最小值/最大值的维护想出了一个算法，发在了 cs.stackexchange 上。现在我觉得那个算法没什么意义了，因为存在一个简单高效的算法，能对一大类聚合解决这个问题。写这篇博客就是为了让更多人知道它——我觉得它应该更广为人知。

民间传说

我是在读一篇更高级的论文时碰到这个算法的：Tangwongsan 等人的 Low-Latency Sliding-Window Aggregation in Worst-Case Constant Time。这篇论文为什么叫低延迟？因为它做的事和我下面要讲的相同，只是每一步都是 O(1)。不过，文中还描述了一个“双栈”算法，均摊 O(1)，而且简单得多，简单得不是一点半点。

有意思的是，那篇论文把这个算法归功于 2011 年 Stack Overflow 帖子里的“adamax”。而后者又说是受 D. Sleator 的一份 2001 年讲义启发。但那份讲义讲的并不是滑动窗口聚合，而是实现 FIFO 队列的经典双栈算法及其均摊分析。说到底，这个算法如此简单又如此巧妙，就算哪天发现它早在 1970 年代某篇冷门论文里出现过，我也不会意外。

两个栈

和论文作者一样，我把双栈算法推广到了任意满足结合律的聚合函数。把聚合抽象成一组函数 empty()、unit(x)、combine(x, y) 和 finalize(x)，就能描述许多种聚合方式，比如求均值：

empty = lambda: (0, 0)
unit = lambda x: (x, 1)
combine = lambda x, y: (x[0] + y[0], x[1] + y[1])
finalize = lambda x: x[0] / x[1] if x[1] else None

这里要说明一下，这些函数的签名是这样的：

fn empty() -> Agg;
fn unit(x: Value) -> Agg;
fn combine(x: Agg, y: Agg) -> Agg;
fn finalize(x: Agg) -> Out;

我特意区分 Value、Agg 和 Out，因为虽然对于整数求和这类场景，它们表面上看起来差不多，但如果是对字符串做近似去重计数，你会得到 (Value, Agg, Out) = (String, HyperLogLogSketch, u64)——三个截然不同的类型。

废话不多说，算法如下：

class TwoStackAgg:
 def __init__(self):
 self.values = []
 self.values_agg = empty()
 self.cum_aggs = []

 def push(self, x):
 self.values.append(x)
 self.values_agg = combine(self.values_agg, unit(x))

 def pop(self):
 if not self.cum_aggs:
 cum_agg = empty()
 while self.values:
 cum_agg = combine(unit(self.values.pop()), cum_agg)
 self.cum_aggs.append(cum_agg)
 self.values_agg = empty()
 self.cum_aggs.pop()

 def eval(self):
 return finalize(
 combine(self.cum_aggs[-1], self.values_agg)
 if self.cum_aggs else self.values_agg
 )

就这些，整个算法就这么多。两个栈（values 和 cum_aggs），外加一个聚合值 values_agg。任意时刻，values_agg 保存的是 values 的聚合结果，而 cum_aggs 保存的是窗口中所有不在 values 里的值的累积聚合，顺序是反的。由此，只要把 cum_aggs 的最后一个值和 values_agg 合并，就能在常数时间内得到整个窗口的聚合结果。

巧妙之处在于（假设 w 是窗口大小），每进行 w 次操作，我们就把 values 全部倒空，一边维护一个滚动聚合值，一边把部分累积聚合压入 cum_aggs。正是这一点让它达到摊还 O(1)：每 w 次 pop 做 O(w) 次内部操作，把每个元素的总工作量限制在 O(1)，尽管单次操作未必是常数时间。

我觉得画出来最清楚。假设我们对 [1, 2, ..., 10] 求和，滑动窗口固定为四个元素，那么每次调用 eval() 后的状态是这样的（values_agg 没画出来，因为它就是这些值的聚合结果）：

 cum_aggs values out
 [] [] = 0
 [] [1] = 1
 [] [1, 2] = 1 + 2
 [] [1, 2, 3] = 1 + 2 + 3
 [] [1, 2, 3, 4] = 1 + 2 + 3 + 4
[4, 3 + 4, 2 + 3 + 4] [5] = 2 + 3 + 4 + 5
 [4, 3 + 4] [5, 6] = 3 + 4 + 5 + 6
 [4] [5, 6, 7] = 4 + 5 + 6 + 7
 [] [5, 6, 7, 8] = 5 + 6 + 7 + 8
[8, 7 + 8, 6 + 7 + 8] [9] = 6 + 7 + 8 + 9
 [8, 7 + 8] [9, 10] = 7 + 8 + 9 + 10
 [8] [9, 10] = 8 + 9 + 10
 [] [9, 10] = 9 + 10
 [10] [] = 10
 [] [] = 0

总内存占用是 O(w)，其中 w 是窗口的最大尺寸。注意，为了分析方便和举例简单，我假设窗口大小固定为 w，但双栈算法本身并不要求这一点。每次 eval() 之间，你可以随意调用 push(x) 和 pop()，按需扩大或缩小窗口。

浮点数不满足结合律

注意，前面我们要求聚合操作 combine 满足结合律，也就是说：

combine(combine(x, y), z) = combine(x, combine(y, z))

严格来说，浮点加法并不满足这一性质。尽管如此，上面的算法依然很有用，因为结果与预期相当接近——如果使用 Kahan 求和这类补偿求和算法，接近程度还会更高。

不过，上面的算法还有第二个很有用的性质：每个聚合值严格只由窗口内的元素组合而成，不涉及窗口之外的任何元素。这意味着，如果你的窗口里出现了 NaN 或无穷大（或者其他离群值），这个值只会污染包含它的那些窗口，而不会影响其余的计算。

即便没有 NaN 或无穷大，这一性质也有用，因为它不会让误差无休止地传播下去。比如，如果你的滑动窗口以 [1e20, 1] 开头，用朴素的滚动求和会得到这样的结果：

>>> 1e20 + 1 - 1e20 - 1
-1.0

补偿求和能减轻这些影响，但让结果不依赖于窗口之外的值，才能彻底消除长期误差累积。