一些动态几何问题的流式算法
本文为 STOC'04 Algorithms for Dynamic Geometric Problems over Data Streams 的阅读笔记。
论文作者 Piotr Indyk， 研究领域：高维几何问题, 流式算法，摘要数据结构维护, 稀疏傅立叶变换。
1 近似算法
在假设 \(\text{P}\neq\text{NP}\) 的情况下，近似算法一般针对 NP 最优化问题(NPO)，有两点要求
-
算法 \(M\) 是确定性或者随机的多项式算法。
-
算法近似界限 \(r\) 要尽可能低。
\[r(n)=\sup_{x\in X,|x|=n}\left\{\max\left\{\frac{M(x)}{M'(x)},\frac{M'(x)}{M(x)}\right\}\right\} \]其中 \(M'\) 为最优解算法，\(X\in\text{NPO}\)。
通过 \(r(n)\) 可以定义一系列近似复杂度类，例如 \(f(n)\text{-APX}\) 表示存在近似算法 \(M\) 使得 \(r(n)\in O(f(n))\) 的问题集合。
一般来说 \(O(1)\text{-APX}\)(可以简写为 \(\text{APX}\)) 算法比较常见。
- \(\texttt{MAX-3SAT}\) 问题有一个 \(r=\frac{8}{7}\) 的算法。并且在 Some Optimal Inapproximability Results(1999) 这篇论文中指出，如果存在 \(\epsilon>0\)，\(r=\frac{8}{7}-\epsilon\) 的多项式算法，则 \(\text{P}=\text{NP}\)。
还有一些问题问题，存在算法可以无限逼近最优解，但是越逼近时空开销越大，所以有 \(\text{PTAS}\subset \text{APX}\)，若对于任意 \(\epsilon>0\), \(\exists\) 时间复杂度为 \(O(f(n,1/\epsilon))\) 的算法，\(f\) 是 \(n\) 的多项式，\(r<1+\epsilon\)。若 \(f\) 为 \(n\) 和 \(1/\epsilon\) 的二元多项式，则为 \(\text{FPTAS}\)。
有时候，若得到的概率分布满足为
近似算法的算法复杂度可能涉及三个参数 \(n\), \(1/\epsilon\), \(1/\delta\)。
一般来说常见的近似算法是贪心等基本的组合算法，另一类常见的算法是规约到线性规划然后采用随机扰动(rounding)，原始对偶，线性松弛等。
另一类研究的比较少的近似算法，是针对 \(\sharp P\) 计数问题，常用的方法是马尔可夫链蒙特卡洛方法。
2 流式算法
流式算法适用的问题一般比 \(\text{NPO}\) 简单，大多数是 \(\text{PO}\) 内，但是空间不足以存下之前的输入或者要求强制在线，需要维护一个摘要或者采样的算法，所以很难得到精确解，同样需要借用用近似算法的评价标准。
早期流式算法研 究为一些简单的统计量的维护：
-
中位数. 普遍情况是求第 \(k\) 小。
流式算法最早研究 Selection and sorting with limited storage（1978）。文中介绍了一个随机数据下空间复杂度为 \(\theta(n ^{1/2})\) 错误率极低的流式随机算法。
文章还严格证明了 \(p\) 次重复读入的确定性算法空间复杂度上下界 \(\Omega(n^{1/p})\) 和 \(O(n^{1/p}\log^{2-2/p} n)\)。 -
不同元素个数 \(F_0\)。普遍地有 k 阶频率矩 \(F_k=\sum\limits_x f(x)^k\).
The space complexity of approximating the frequency moments(1996)文章证明了 \(F_0, F_1, F_2\) 存在 \(O(\log n)\) 空间的逼近算法，用通信复杂度证明了 \(F_k(k\ge 6)\) 则需要 \(n^{\Omega(1)}\) 的空间。
文章用到了均匀采样，随机变量，切比雪夫不等式，切尔诺夫界等概率论工具，以及通信复杂度，将流式随机算法的研究形式化，并且将这些工具广开来成为现在的主流研究方法，可以算这个领域的开山之作。
文章构造了 \(2\log(1/\delta)\) 个随机变量 \(Y_i\)，\(E(F_k)=E(\overline{Y_i})\)，然后 \(Y_i=\overline{X_{ij}}\) 为 \(8kn^{1-1/k}/\epsilon^2\) 个随机变量的均值，\(X\) 为互相独立的均匀采样集生成的随机变量，方法为 \(X=m(r^k-(r-1)^k)\), \(r=|\{q|q\ge p,a_p=a_p\}|\)，其中 \(p\) 为随机抽样的下标。维护这样的 \(X\) 和 \(Y\) 空间复杂度为 \(O\left(\frac{k\log(1/\delta)}{\epsilon^2}n^{1-1/k}(\log n+\log m)\right)\)。
文中用切比雪夫不等式证明了单个 \(Y_i\) 相对误差大于 \(\epsilon\) 的概率小于 \(\frac{1}{8}\)，根据切尔诺夫界，可以的到 \(2\log(1/\delta)\) 个 \(Y_i\) 的平均值的越界概率小于 \(\delta\)。文章针对 \(F_2\) 设计了专门的随机变量，得到了更优的结果。随机变量估计的一个另一个更早的例子是，Probabilistic counting(FOCS,1983) 中提出算法，\(E(\log F_0)=E(\max(\text{lowbit}(\text{hash}(x))))\)，将统计量 \(F_0\) 用 \(\max(\text{lowbit}(\text{hash}(x)))\) 进行无偏估计，然后用一些概率论的技巧计算出其错误率，然后用 Median trick 可以降低错误率。
在后续的研究中有人引入和稳态分布和随机投影。
2.1 采样
采样最早是从统计学和数据科学中发展来的方法，看起来很简单，但是合理的采样方法可以在流式算法中起到出乎意料的作用。
- 用蒙特卡罗方法计算复杂函数的积分或期望，利用均匀分布或其他简单分布生成随机变量，然后用样本均值近似真实值。
- 用接受-拒绝抽样或重要性抽样方法生成难以直接采样的分布的样本，利用一个容易采样的分布作为提议分布，然后根据一定的准则接受或拒绝样本，或者给样本赋予不同的权重。
- 用马尔可夫链蒙特卡罗方法生成难以直接采样的分布的样本，利用一个马尔可夫链逐步逼近目标分布的稳态分布，然后从稳态分布中采样。
2.2 随机投影
随机投影是高维问题中经常使用的方法。The Johnson-Lindenstrauss 引理指出高维空间中的一小部分点可以以点之间的距离几乎被保留的方式嵌入到低维空间中。
通常的方法是把一些高维的范数分解成多个一维和二维的稳态分布来近似。
3 流式问题复杂度上下界证明
3.1 通信复杂度
通信复杂度研究的是，将一个问题的输入划分到两个或多个图灵机，假设所有图灵机有无穷的计算能力且彼此互信，协议相同，求解决问题需要传输的最小比特量与输入规模之间的关系。
一个经典的协议是，两台图灵机可以在传递 \(O(\log n)\) 个比特的情况下，计算出输入的中位数。
通信复杂度在很多领域例如计算复杂性理论，量子浅层析
的采样复杂度下界，电路复杂度下界等领域有重要的作用。
在数据结构空间下界和流式算法空间下界中的应用方法一般是
一种方法是将流式问题转化为分布式问题，即将数据流分配给多个参与者，让他们各自处理自己的数据，并且在必要时进行通信，最后输出结果。这样，流式问题的内存空间限制就变成了分布式问题的通信限制。
另一种方法是使用信息论的工具，即将数据流看作是随机变量，利用熵、互信息等概念来量化流式问题中所需的信息量，从而得到通信复杂度的下界。
3.2 可压缩性分析
可压缩性分析的核心思想是将流式算法与有限状态自动机(DFA) 建立关系，利用 Myhill–Nerode 定理来寻找输入的等价类数目，确定 DFA 可压缩的下界从而刻画出流式算法空间复杂度的下界。
Stanford University CS154 的 Note 中使用这种方法重新将频率矩 \(F_k\) 下界复杂度的结论证明了一遍。
3.3 规约
规约是在计算复杂性理论和精细复杂度领域常用的非常技巧。
4 动态几何问题
4.1 离散高维度量空间
本文研究的问题是在离散 \(d\) 维度量空间 \(D=\{1\cdots\Delta\}^d\) 中动态加点或删点，维护点集的一个函数，其中中点之间的距离被定义为闵可夫斯基范数 \(d(x,y)=\left(\sum\limits_{i=1}^d|x_i-y_i|^p\right)^{1/p}\)。其中 \(\Delta\) 表示不同元素的个数。
这类问题首先在论文 The Geometry of Graphs and Some of its Algorithmic Applications(1994) 中被探讨。文章提出了一些高效算法，在低维度量空间中以很小的失真嵌入图。将一些图上难解的问题例如直径，最小割问题转化为几何问题，使用更高效的近似算法来解决。
Probabilistic approximation of metric spaces and its algorithmic applications(1996) 。这篇文章提出了任何度量空间都可以分层良分树(HST)以 polylog 失真概率近似。
在 Approximating a finite metric by a small number of tree metrics(1998)，这篇文章中使用了线性规划这种确定性算法构造了大小为 \(O(n\log n )\) 的 HST，使得每条边期望失真不超过 \(O(\log n \log \log n)\)。
本文同样使用了 HST，介绍了几个高维度量空间问题的流式算法。
4.2 最小生成树问题
Approximating the minimum spanning tree weight in sublinear time(2001) 文中给出了图上最小生成树的亚线性任意近似算法，时间复杂度为 \(O(\frac{dw}{\epsilon^2}\log\frac{d}{\epsilon})\)，接近下界 \(\Omega(\frac{dw}{\epsilon^2})\)，其中 \(d\) 为点的平均度数，\(w\) 为边权集合的大小。
The sensor spectrum: technology, trends, and requirements(SIGMOD, 2003) 中对传感器网络通信代价的研究揭示了高维几何中最小生成树问题的实用意义亚线性任意近似算法，时间复杂度为 \(O(\mathrm{polylog}(n/\epsilon^{O(1)}))\).
同年 Estimating the Weight of Metric Minimum Spanning Trees in Sublinear-Time(2003) 文中的到了高维几何(度量空间)中中最小生成树问题的
基于前面文章的工作，在 Estimating the weight of euclidean minimum spanning trees in data streams(2004) 中提出了只支持插入情况下，最小生成树问题的任意逼近流式算法。
本文对这个问题得到了在同时支持插入和删除操作的情况下以 \(r=O(d\log \Delta)\) 逼近，空间为 \(O(d\log^2\Delta)\) 的流式算法。
Chen, Jayaram, Levi, Waingarten 将这个问题改进到 \(r\le 1.10\)，空间为 \(\Omega(\sqrt{n})\)。
4.3 匹配问题
空间中匹配问题分为两种
- 无色匹配(Minimum Weight Matching)
- 双色匹配(Minimum Bi-chromatic Matching)
本文作者指出 Similarity estimation techniques from rounding(1996) 这篇文章的方法足以解决双色匹配问题。
本文给出了无色匹配的 \(r=O(d\log \Delta)\)，空间为 \(O(d\log^2\Delta\log\log\Delta)\) 的算法。
文中提到，经过仔细分析，可以做到 \(r=1+\epsilon\)，且空间复杂度只乘上 \(O(1/\epsilon^2)\) 因子的算法，但是没有在本文列出。
4.4 工厂选址问题
这类问题由在铁路上的传感器通信问题发展而来。
本文给出了一个 \(r=O(d\log^2\Delta)\)，空间为 \(O(d^2\log^2\Delta)\) 的算法。
在 Streaming Facility Location in High Dimension via New Geometric Hashing(2022) 中，提出了一种基于重要性采样的算法，改进到了双次扫描 \(r=O(1)\)，空间为 \((d\log \Delta)^{O(1)}\)，单次扫描 \(r=O(d/\log d)\)。同时这篇文章给出了一种空间为 \(O(n^{1/\epsilon})\) 的任意近似算法。
4.5 k-聚类问题
k-聚类问题是几何问题中非常经典的 NPO 问题，在机器学习等众多领域有着很重要的作用，已经有众多有效的启发式算法和近似算法。
在作者所在的时间，已经提出了 \(r=3+\epsilon\)，时间为 \(O(n^{1/2\epsilon})\) 的非流式算法。目前最优的逼近比是 \(r=2.406\)，由 Improved approximations for Euclidean k-means and k-median, via nested quasi-independent sets(2022) 文中提出。
本文使用了多种工具和算法例如互斥计数，中点成本估计，本地搜索，贪心等。分别做到了不同的复杂度。
在 The Power of Uniform Sampling for k-Median(2023) 中，作者分析了朴素的均匀采样的方法，如果要达到 \(r=O(1)\)，则采样率至少为 \(\Omega(1/\beta)\)，\(\beta\) 为均匀度。如果要达到任意近似，则至少达到 \(O(\mathrm{poly}(k/\epsilon\beta))\) 采样率。
5 科技
5.1 嵌入分层良分树(HST)
如果一个有根树 \(T\) 是 k-HST，则满足
- 每个节点到所有子节点的距离相等
- 从根到叶子节点的路径上，边权至少以 \(k>1\) 的比例持续减少。
一个 HST 上的度量空间可以定义为其上所有叶节点，而度量为叶节点之间的最短路。
本文引用了前人的结论，任意 \(\{1\cdots\Delta\}^d\) 的子集 \(P\) 上的度量可以被嵌入一系列 2-HST 度量的概率组合，且失真率小于 \(O(d\log\Delta)\)。
文中以该随机算法构建 2-HST
首先对 \(P\) 中所有点增加一个随机偏移 \(\Delta p\in [0,\Delta]^d\)
树分为 \(\log\Delta +2\) 层 \(G_0,\cdots, G_{\log\Delta+1}\)。
- \(G_0\) 为 \(P\) 中的点为叶节点。\(G_{\Delta +1}\) 表示根节点。
- \(G_i\) 表示包含了 \(P\) 中点的大小为 \(2^{i-1}\) 的网格。
- \(G_i\) 与 \(G_{i+1}\) 之间按包含关系连边，边权为 \(2^i\)。
5.2 规约到高维向量统计问题
这是本文创新提出的方法，令统计量 \(\pi_S(c)\) 表示点集 \(S\) 落在网格 \(c\) 中的数量。这种统计量前人已经研究的很多了。本文用这个统计量的组合来逼近要求的问题，将近似最优化问题转化为 \(\log\Delta\) 个近似计数问题。
- 最小双色匹配近似为 \(\sum\limits_{i=1}^{\log\Delta}2^i\left(\sum\limits_{c\in G_i}|\pi_R(c)-\pi_B(c)|\right)\)，其中 R, B 为两种颜色的点集。
- 最小生成树近似为 \(\sum\limits_{i=1}^{\log\Delta}2^i\left(\sum\limits_{c\in G_i}[\pi_P(c)>0]\right)=\sum\limits_{i=1}^{\log\Delta}2^i|G_i|\)。
- 最小匹配近似为 \(\sum\limits_{i=1}^{\log\Delta}2^i\left(\sum\limits_{c\in G_i}[\pi_P(c)\equiv 1\pmod 2]\right)\)。
- 工厂选址问题最小代价近似为 \(\sum\limits_{i=1}^{\log\Delta}\left(\sum\limits_{c\in G_i}\min\{2^i\times\pi_P(c),f\}\right)\)，\(f\) 是每个工厂自身的代价。
- 给定方案的 k-聚类问题的代价近似为 \(\sum\limits_{i=-i_0}^{\log\Delta}(r_{i}-r_{i-1})\left(\sum\limits_{c_i\cap B(Q,r_i)=\emptyset}\pi_P(c_i)\right)\)，其中 \(r_i=(1+\epsilon)^i\)，\(c_i\) 的直径为 \(\epsilon/\sqrt{2}\times(1+\epsilon)^i\)，\(i_0=O(\log(1/\epsilon)/\epsilon)\)，\(B(Q,r)\) 表示 \(Q\) 中所有点为中心，半径为 \(r\) 的高维球区域。
6 算法
6.1 MST
由 Lemma7.1 可以把问题转化为求 \(w(T)\)。通过定义和恒等变换，很容易得到 5.2 中的式子。
考虑统计每个 \(|G_i|\)，发现可以等价为统计不同元素的个数 \(F_0\)，根据前人的结论可以做到 \(O(\log(n+\Delta^d))=O(d\log\Delta)\) 空间，于是总的空间为 \(O(d\log^2\Delta)\)。
6.2 MWM
根据 Lemma7.2 可以把问题转化为 \(\log\Delta\) 个 Odd Count 问题。
维护长度为 \((\frac{\Delta}{2^i})^d\) 的数组 \(x_i\)
- \(\text{Update}(i,c)\) : \(x_i\leftarrow x_i+c,(c\in\{1,-1\})\)
- \(\text{OddCount}\) : return \(\sum [x_i\equiv 1\pmod 2]\)
文中提出了一种 \(r=O(1)\)，空间复杂度为 \(O(\log\Delta)\) 的算法，有常数概率出错。如果将算法并行运行 \(O(\log\log\Delta)\) 取均值，可以做到超出界限的概率无限减小变为 \(\frac{1}{2\log\Delta}\)。
文章首先构造了一个概率判定算法
M="On input \(OC,T\)，\(OC\) 是一个流式 Odd Count 问题，\(T\) 是参数
- 初始以 \(1/T\) 的等价概率选取 \(x\) 下标的一个子集 \(R\)。 \(s\leftarrow 0\)。
- \(\text{Update}(i, c)\)，如果 \(c\in R\), \(s\leftarrow 2 -s\)。
- \(\text{OddCount}\) : 输出 \(s\) 是否为 \(1\)。
这个算法满足（证明见 Lemma7.3）
- 如果真实 \(OC>T\) 则以不小于 \(1/3\) 的概率输出
YES。 - 如果真实 \(OC<T/10\) 则以不大于 \(1/10\) 的概率输出
YES。
通过作者本人之前的文章 pseudorandom generators, embeddings and data stream computation(2000) 中给出的方法，可以用伪随机数生成器以 \(O(d\log \Delta)\) 的空间存储这个随机集合。
本文后续没有介绍如何通过判定问题近似算法推出计数问题的算法，推测是等比构造了 \(\log \Delta\) 个 \(T\) 将近似比约束在 \(O(1)\) 的范围内，然后伪随机数可以在不同的并行运算中复用。
总空间复杂度为 \(O(d\log^2\Delta\log\log\Delta)\)，近似度为规约到 HST 的 \(O(d\log\Delta)\)。
6.3 工厂选址
根据 Lemma7.4，工厂选址问题可以以 \(O(d\log^2\Delta)\) 的近似度规约到 \(\log\Delta\) 个 Bounded Count 问题
维护 \((\frac{\Delta}{2^i})^d\) 个点集 \(S_x\)。
- \(\text{Add}(x,p)\): \(s_x\leftarrow s_x\cup\{p\}\)。
- \(\text{Delete}(x,p)\): \(s_x\leftarrow s_x/\{p\}\)。
- \(\text{BoundedCount}\): return \(\sum\limits_{x}\min\{|S_x|,T\}\)，\(T\) 是参数。
文章对 Bounded Count 给出的算法是用伪随机数生成器构造了一个值域为 \(T\) 的哈希函数 \(h\)，然后从用这个哈希函数从流中等价采样了 \(1/T\) 种不同的值，一旦任意一个值落入某个集合，则这个集合按满的答案 \(T\) 算，这样问题便转化为了求不同元素 \(F_0\)。这种转化可以保证失真在 \(O(1)\) 内，证明见 Lemma7.5。
在作者自己之前的论文中指出，生成随机哈希函数的伪随机数生成器需要消耗 \(O(\log^2 M)\)。在 \(\log\Delta\) 个并行运算中这些可以复用，总空间复杂度为 \(O(d^2\log^2\Delta)\)。
6.4 k-聚类的搜索算法
根据 Lemma7.6，存在 \(r=(1+\epsilon)\) 的近似算法，针对给定的中心点求出 k-median 的代价。
局部搜索也称为爬山法，是一种很常见的启发式算法，缺点是容易陷入局部最优解。
本文引用了 Local search heuristics for k-median and facility location problems(2002) 中的结论证明了在 k-median 问题上局部搜索可以做到 \(r=5\)。
再引入计算代价的失真，总体可以做到 \(r=5+\epsilon\)，空间复杂度 \(O(k\cdot\mathrm{polylog}(\Delta))\)。
如果换成全局搜索可以做到 \(r=1+\epsilon\)。
6.5 k-聚类的贪心算法
贪心法的思路很简单，初始让集合为空，每次选择使得代价减少最多的点。
文中为了降低复杂度，对贪心做了一些优化。
7 证明
Lemma7.1
若 MST 为 \(\mathcal{T}\)，生成的 HST 为 \(T\), \(\frac{w(T)}{w(\mathcal{T})}\le O(d\log\Delta)\)。
因为 HST 的两点距离不大于原来两点距离的 \(O(d\log\Delta)\) 倍
\[w(T)\le\sum_{u,v\in\mathcal{T}}d_T(u,v)\le O(D\log\Delta)w(\mathcal{T}) \]
Lemma7.2
HST 上最小匹配等于 \(n+\sum\limits_i^{\log\Delta}2^im_i\)，\(m_i\) 为 \(G_i\) 的 odd count.
考虑 HST 上每个节点到父节点的边经过的次数，容易发现，如果经过次数大于 \(1\) 次可以将这两组匹配交叉使得答案更优。所以最多经过一次，且仅当子树叶子节点数量为奇数时才需要。
Lemma7.3
6.2 中的判定性问题满足其中提到的概率条件。
本文引用了 Efficient search for approximate nearest neighbor in high dimensional spaces 中的 Lemma2.1 证明方法，表示两者相似。
Lemma7.4
如果 HST 上工厂选址问题最优代价是 \(C\)， \(Q=\sum\limits_{i=1}^{\log\Delta}\left(\sum\limits_{c\in G_i}\min\{2^i\times\pi_P(c),f\}\right)\) 是 \(C\) 的一个 \(O(\log\Delta)\) 近似，满足
首先，构造一个代价为 \(Q+f\) 的方案，在根节点放一个工厂，然后对所有 \(2^i\times\pi_P(c)\ge f\) 的节点放置一个工厂，这样的代价为 \(C'\)，\(C\le C'<Q+f\)。
对于等式的另一边，容易知道对于每一层 \(\sum\limits_{c\in G_i}\min\{2^i\times\pi_P(c),f\}\le C\)，所以总共 \(\log\Delta\) 层加起来小于 \(C\log\Delta\)。
Lemma7.5
6.3 中的算法是 Bounded Count 的 \(O(1)\) 近似算法。
若 \(K\) 表示有最终统计出的不同元素个数，\(p(k)=1-(1-1/T)^k\) 表示某个集合有 \(k\) 元素且没有全部被哈希函数筛出的概率。
\[E(K)=\sum_i p(|S_i|) \]\[\text{Var}(K)=\sum_i p(|S_i|)-p^2(|S_i|)\le E(K) \]根据切比雪夫不等式
\[P(|K-E(K)|\ge \epsilon E(K))\le \frac{\text{Var}(K)}{\epsilon^2E(x)^2}\le \frac{1}{\epsilon^2} \]
Lemma7.6
设给定 k-中心点 \(Q\)，k-聚类代价为 \(C(Q,P)\)。
设 \(R_i(Q)=|P-B(Q,r_i)|\)，\(R_i(Q)\) 可以用本文提出的 Exclusive Count 算法以 \(r=(1+\epsilon)\)，\(\delta=\frac{1}{2}\)，空间为 \(\Omega(k)\) 得到近似 \(\hat{R_i}(Q)\)，满足
文章使用 \(\hat{C}(P,Q)=(r_i-r_{i-1})\hat{R}_i(Q)\) 近似 \(C(P,Q)\)。
所以 \(\hat{C}(Q,P)\) 是 \(C(Q,P)\) 的一个 \((1+\epsilon)\) 近似。
