UCB-CS294-数据流算法笔记-全-

UCB CS294 数据流算法笔记(全)

001:课程介绍与近似计数

概述

在本节课中,我们将学习数据流算法的基本概念,并从一个非常简单的例子——近似计数问题——开始。我们将探讨为什么需要近似算法,如何设计它们,以及如何分析其内存使用和准确性。


课程安排与物流

本课程的所有信息,包括教学大纲、评分标准和联系方式,均已在线发布。如果你在课程注册或获取资料方面遇到任何问题,请通过课程网站或电子邮件联系我。


数据流算法简介

数据流算法处理的是以连续“流”形式到达的数据。我们无法存储整个数据集,因此需要设计能够在有限内存下实时处理查询的算法。

让我们从一个具体的例子开始:计算当天的总销售额。

问题:有一个销售交易数据库,包含数百万条交易记录。我们只想知道一个查询的结果:今天的总销售额(以美元计)。

简单解法:维护一个计数器,初始为0。每发生一笔销售,就将金额加到计数器上。这是一个简单的数据流算法。

类似地,如果我们只想知道交易笔数(而非总金额),也可以维护一个计数器。如果有 n 笔交易,一个精确的计数器需要 log n 比特的内存。


近似计数问题

上一节我们介绍了精确计数器的概念。本节中,我们来看看如何通过“近似”来做得更好。

我们将问题形式化为一个数据结构问题,需要支持三种操作:

  1. 初始化:将计数器设为0。
  2. 递增:将计数器加1。
  3. 查询:返回计数器值的一个近似值

精确算法的局限性

如果我们要求算法是确定性的且输出精确值,那么可以证明,当计数器值 n0T 之间时,算法至少需要 log T 比特的内存。

证明思路:如果使用少于 log T 比特的内存,根据鸽巢原理,必然存在两个不同的递增次数 AB,使得算法最终的内存状态完全相同。因此,对于查询操作,算法会对 AB 返回相同的答案,而其中至少有一个答案是错误的。

引入近似与随机化

由于精确解需要对数级内存,我们考虑放松要求:

  1. 允许近似:我们接受答案在真实值 n(1 ± ε) 范围内(例如,1%误差)。
  2. 允许随机化:算法可以有一定的失败概率 δ,即答案不满足误差要求的概率。

我们的目标是设计一个随机化的近似算法,在保证 (1 ± ε) 近似精度的前提下,使用尽可能少的内存。

理想的内存下界

即使允许近似,我们也需要一定的内存来区分不同的可能值。

假设真实值 n0T 之间。为了达到 (1 ± ε) 的近似,我们实际上只需要知道 n 最接近 (1+ε) 的哪个幂次。例如,当 ε=0.01 时,我们需要区分的可能值是:1, 1.01, 1.01^2, 1.01^3, ..., T

这些可能值的数量大约是 log_(1+ε) T。存储其中某一个的索引所需的内存是:
log( log_(1+ε) T ) ≈ log(log T) + log(1/ε) 比特。

这给出了一个理想的内存下界:任何 (1±ε) 近似算法至少需要 Ω( log(log n) + log(1/ε) ) 比特的内存。我们后续的算法将努力接近这个下界。


Morris 计数器算法

上一节我们讨论了问题的下界。现在,我们来看一个具体的算法——Morris 计数器。

基础想法

我们不直接存储计数器 n,而是存储另一个变量 x

  • 初始化x = 0
  • 递增:以概率 1/(2^x)x 加1,否则 x 保持不变。
  • 查询:返回估计值 ñ = 2^x - 1

为什么这可能有效?

  • 内存x 的增长非常缓慢。可以证明,x 超过 O(log log n) 的概率极低,因此存储 x 所需的内存约为 log log n 比特。
  • 无偏性:可以证明,估计值 ñ 的期望恰好是 n,即 E[2^x] = n + 1

算法分析(期望)

以下是 E[2^x] = n + 1 的归纳证明:

  • 基础情况n=0 时,x=02^0 = 1 = 0+1
  • 归纳步骤:假设 E[2^(X_n)] = n+1 成立。考虑第 n+1 次递增后:
    E[2^(X_(n+1))] = Σ_j Pr[X_n = j] * E[2^(X_(n+1)) | X_n = j]
    X_n = j 时,有概率 1/(2^j)x 增至 j+1,概率 1 - 1/(2^j) 保持 x=j
    因此,E[2^(X_(n+1)) | X_n = j] = (1/(2^j)) * 2^(j+1) + (1 - 1/(2^j)) * 2^j = 2^j + 1
    代入求和式可得:E[2^(X_(n+1))] = 1 + E[2^(X_n)] = 1 + (n+1) = n+2
    由归纳法,结论成立。

因此,E[ñ] = E[2^x - 1] = (n+1) - 1 = n。我们的估计量是无偏的。

算法的不足

基础的 Morris 计数器方差较大。通过分析(此处省略详细计算,可参考课程笔记),其方差 Var[ñ] ≈ n^2 / 2。根据切比雪夫不等式,其误差超过 εn 的概率上界为 1/(2ε^2)。当 ε=0.75 时,失败概率高达约 8/9,成功率仅约 11%。这并不理想。


改进:Morris+ 与 Morris++ 算法

上一节的 Morris 计数器方差太大。本节我们介绍两种通用的改进技术,它们可以应用于许多随机化近似算法。

Morris+:取平均值降低方差

一个降低方差的标准方法是运行多个独立副本并取结果的平均值。

算法

  1. 并行运行 r 个独立的 Morris 计数器副本,得到估计值 ñ_1, ñ_2, ..., ñ_r
  2. 输出最终估计值 n̂ = (ñ_1 + ... + ñ_r) / r

分析

  • 期望E[n̂] = n(仍是无偏估计)。
  • 方差:由于副本独立,Var[n̂] = Var[ñ] / r ≈ n^2 / (2r)
  • 误差概率:根据切比雪夫不等式,Pr[|n̂ - n| > εn] ≤ Var[n̂] / (ε^2 n^2) ≈ 1/(2rε^2)
    为了使失败概率小于 η,我们只需设置 r ≥ 1/(2ε^2η)

代价:内存使用增加了 r 倍,总内存约为 O( (log log n) / (ε^2 η) )。对 εη 的依赖不够好(我们目标是 log(1/ε)log(1/η))。

Morris++:取中位数降低失败概率

Morris+ 对失败概率 η 的依赖是 1/η 倍的。我们可以通过“取中位数”的技术将其改进为对 log(1/η) 的依赖。

算法

  1. 并行运行 s 组 Morris+ 算法,每组设定其失败概率 η’ 为一个常数(例如 1/3)。每组输出一个估计值 n̂_i
  2. 输出最终估计值 ñ = median( n̂_1, n̂_2, ..., n̂_s ),即所有 n̂_i 的中位数。

分析

  • 每组 Morris+ 的失败概率为常数 η’ = 1/3
  • 最终输出 出错,仅当超过一半的 n̂_i 出错。
  • X_i 为指示变量,当 n̂_i 出错时为1。则 X_i 独立,且 E[Σ X_i] ≤ s/3
  • 我们需要 Σ X_i > s/2 才会失败,这要求其和偏离期望至少 s/6
  • 根据 切尔诺夫界:对于独立的 {0,1} 随机变量,其和偏离期望 εμ 倍的概率上界为 2 * exp(-ε^2 μ / 3),其中 μ 是期望和。
  • 在此设置中,μ ≤ s/3, εμ ≥ s/6,可得失败概率上界为 2 * exp(-Ω(s))
  • 因此,只需选择 s = O( log(1/δ) ),即可将最终失败概率降至 δ

总内存O( (log log n) * (1/ε^2) * log(1/δ) )。相比 Morris+,我们将对 δ 的依赖从 1/δ 改进为了 log(1/δ),但对 ε 的依赖仍是 1/ε^2


最优的 Morris 计数器变体

上一节介绍的“取中位数”方法是通用技术,但并非对本问题最优。Morris 在其原始论文中提出了更优的变体。

核心思想:调整递增概率

回顾基础 Morris 计数器:以概率 1/(2^x) 递增 x

  • 概率 1/(2^x) 较小,导致 x 增长慢,内存省,但方差大。
  • 如果我们以概率 1/(1.01^x) 递增呢?x 增长更快,内存消耗稍多,但方差会显著减小。
  • 极端情况:以概率 1 递增,则 x 始终等于 n,方差为0,但内存变为 log n

可见,在内存方差之间存在权衡。

最优算法描述

更一般地,我们可以设计算法,以概率 1/( (1+α)^x ) 递增 x,其中 α 是一个与目标精度 ε 和失败概率 δ 相关的参数。

经过恰当的参数设置(例如,α = Θ(ε^2 δ))和细致分析(可作为课后练习),该算法可以达到以下内存上界:
O( log log n + log(1/ε) + log(1/δ) ) 比特。

这已经比 Morris++ 的结果更优,因为对 ε 的依赖从 1/ε^2 改进为了 log(1/ε)

当前最优界

最新的研究(2020年夏季)表明,近似计数问题的最优内存界限实际上是:
min( log n, O( log log n + log(1/ε) + log log(1/δ) ) ) 比特。

解释

  1. log n 是确定性精确计数器的成本。
  2. 括号内的项是随机化近似计数器的成本。关键改进在于对失败概率 δ 的依赖从 log(1/δ) 进一步优化为了 log log(1/δ)
  3. min 操作意味着:当要求的精度 ε 极高或失败概率 δ 极低时,直接使用确定性精确计数器可能更省内存;否则,使用优化后的 Morris 类算法更优。

这个上界是紧的,存在匹配的下界证明。


总结

本节课我们一起学习了数据流算法的入门知识,并深入探讨了近似计数这个经典问题。

  1. 问题定义:我们希望在有限内存下,维护一个计数器,支持递增和近似查询操作。
  2. 下界:任何 (1±ε) 近似的算法至少需要 Ω(log log n + log(1/ε)) 比特内存。
  3. Morris 计数器:通过存储对数级的变量 x 并以递减的概率更新它,实现了 O(log log n) 的内存使用,并给出了无偏估计。
  4. 通用改进技术
    • Morris+(取平均):通过并行多个副本取平均来降低方差,代价是内存乘以 O(1/ε^2)
    • Morris++(取中位数):通过对多个“常数失败概率”的估计器取中位数,将对失败概率的依赖从 O(1/δ) 指数级改进到 O(log(1/δ))
  5. 最优变体:通过调整 Morris 计数器中递增概率的底数,可以直接得到更优的内存界限 O(log log n + log(1/ε) + log(1/δ))
  6. 最新进展:该问题的最优内存界限为 min( log n, O(log log n + log(1/ε) + log log(1/δ)) ),并且上下界已匹配。

近似计数是数据流算法中最古老、最基础的问题之一,其中蕴含的“概率更新”、“方差-内存权衡”、“中位数技巧”等思想,是设计更复杂数据流算法的基石。

002:独立元素计数与k-wise独立性

在本节课中,我们将学习如何高效地统计数据流中不同元素的数量。我们将从一种理想化的算法开始,然后探讨如何通过使用伪随机哈希函数(特别是k-wise独立哈希函数)来移除理想化假设,从而在有限内存下实现近似计数。


概述

本节课我们将重点学习“独立元素计数”问题。给定一个由整数组成的数据流,我们的目标是估算流中不同整数的数量,同时尽可能减少内存使用。我们将首先介绍一个基于理想化随机哈希函数的简单算法(FM算法),然后分析其局限性。接着,我们将探讨如何通过使用k-wise独立哈希函数来移除“理想化”的假设,从而得到一个在实际中可实现的算法(KMV算法)。核心在于理解如何用有限的随机性(少量随机比特)来模拟完全随机的行为,并保证算法的正确性。


输入、查询与朴素解法

输入流是一个整数序列:I1, I2, ..., Im,其中每个整数 Ii 的取值范围是 1n

查询的目标是:流中不同整数的数量。

例如,对于流 [1, 1, 5, 1, 3, 7, 4, 4],不同元素是 {1, 3, 4, 5, 7},因此答案是 5

存在两种朴素的解决方案:

  1. 位向量法:维护一个长度为 n 的位向量。当看到整数 i 时,将第 i 位设为 1。查询时,统计位向量中 1 的个数。这需要 n 比特的内存。
  2. 存储整个流:直接记住所有 m 个整数。每个整数需要 log₂n 比特,总共需要 m * log₂n 比特的内存。

我们的目标是设计一个算法,其内存消耗远低于 min(n, m),最好是与 log nlog m 成正比。


理论下界与FM算法(理想化版本)

理论研究表明,如果不允许近似和随机化,就不存在远低于 min(n, m) 内存的算法。因此,我们将设计一个随机化的近似算法。我们的输出 需要满足:真实值 T 与估计值 之间的相对误差超过 ε 的概率最多为 δ

我们首先介绍一个由Flajolet和Martin在20世纪80年代提出的理想化算法,称为FM算法

算法描述(理想化FM)

  1. 初始化
    • 选择一个完全随机的哈希函数 h,它将域 {1, ..., n} 均匀地映射到连续区间 [0, 1]
    • 初始化一个计数器 X = 1
  2. 更新(看到整数 i
    • X = min(X, h(i))
  3. 查询(估计不同元素数 T
    • 返回 T̃ = 1/X - 1

为什么这个算法是“理想化”的?
因为它假设我们可以“免费”存储并使用一个完全随机的哈希函数 h。实际上,描述这样一个函数需要存储 n 个随机实数,这比朴素位向量法的 n 比特还要昂贵得多。我们稍后会解决这个问题。

算法原理与期望分析

假设流中真实的独立元素集合是 {I1, I2, ..., IT}。那么,X 存储的就是这些元素哈希值的最小值:X = min{h(I1), h(I2), ..., h(IT)}

由于 h 是均匀随机的,h(I1), ..., h(IT)T 个在 [0,1] 区间上独立同分布的均匀随机变量。

关键观察:这 T 个随机变量在 [0,1] 区间上大致均匀分布。最小值 X 的期望值大约是 1/(T+1)

我们可以精确计算 X 的期望和方差:

  • 期望E[X] = 1/(T+1)
  • 方差Var[X] ≤ (E[X])^2 = 1/(T+1)^2

因此,如果我们用 1/X - 1 来估计 T,在期望上是正确的。但方差的存在意味着单次估计可能不准确。


提升精度与鲁棒性:FM+ 与 FM++

与Morris计数器类似,我们可以通过平均和取中位数的方法来提升FM算法的精度(减小 ε)和鲁棒性(减小失败概率 δ)。

FM+ 算法:通过平均减少方差

我们独立运行 r 份FM算法副本,得到估计值 X1, X2, ..., Xr
查询时,我们计算它们的平均值 X_avg = (X1 + ... + Xr) / r,然后输出 T̃ = 1/X_avg - 1

由于方差随着副本数量 r 的增加而线性减少(Var[X_avg] = Var[X] / r),根据切比雪夫不等式,我们可以通过设置 r = O(1/(ε²δ)) 来保证 (1±ε) 近似,且失败概率不超过 δ

FM++ 算法:通过中位数控制失败概率

我们独立运行 s 份FM+算法副本,每份的失败概率设为 1/3
查询时,我们收集每份副本的估计值 T̃1, T̃2, ..., T̃s,然后输出它们的中位数。

根据切尔诺夫界,只要 s = O(log(1/δ)),中位数给出错误答案的概率就能被控制在 δ 以内。

综合空间消耗:FM++ 需要存储 O(r * s) = O( (1/ε²) * log(1/δ) ) 个实数。这虽然比朴素解法好,但仍依赖于理想化的哈希函数。


移除理想化假设:KMV算法与k-wise独立性

上一节我们介绍了基于理想化哈希函数的FM算法族。本节中,我们来看看如何移除“哈希函数完全随机且免费”这个不切实际的假设。核心思想是使用k-wise独立哈希函数族

KMV 算法思想

KMV算法(K Minimum Values)由Bar-Yossef、Jayram、Kumar、Sivakumar和Trevisan提出。它做出了一个不同的权衡:不再只存储最小的哈希值,而是存储k个最小的哈希值

设我们存储的k个最小哈希值为 X = {x1, x2, ..., xk}(已排序)。
查询时:

  • 如果看到的独立元素数 T < k,我们可以精确回答(因为存储了所有元素的哈希值)。
  • 否则,我们返回 T̃ = k / xk - 1xk是第k小的哈希值)。

直观理解:如果 T 个均匀随机变量大致均匀分布在 [0,1],那么第 k 小的值期望约为 k/(T+1)。因此,k/xk 期望约为 T+1,减去1就得到了 T。存储更多的最小值(使用更高阶的统计量)比只存储最小值(FM算法)具有更小的方差,但需要更多内存。这是一种在内存和估计方差之间的权衡。

伪随机性与k-wise独立性

我们无法存储一个完全随机的哈希函数(需要 O(n) 内存)。伪随机性的目标是:用很少的随机比特(种子)来生成一个表现如同随机的哈希函数,足以“欺骗”我们的算法。

k-wise独立哈希函数族的定义:
一个哈希函数族 Hk-wise独立 的,如果对于任意 k 个不同的输入 x1, x2, ..., xk,以及任意 k 个输出 y1, y2, ..., yk,都有:
Pr[ h(x1)=y1 ∧ h(x2)=y2 ∧ ... ∧ h(xk)=yk ] = 1 / M^k
其中概率来源于从 H 中均匀随机选择哈希函数 hM 是输出范围的大小。

这意味着:任意最多 k 个输入点的哈希值,联合起来看是完全独立且均匀分布的

一个重要例子:多项式哈希
设定义域和值域大小均为一个质数 q。哈希函数族由所有次数小于 k 的多项式构成:
h(x) = a0 + a1*x + a2*x^2 + ... + a_{k-1}*x^{k-1} mod q
其中系数 a0, ..., a_{k-1} 在有限域 F_q 中随机选择。
这个族是 k-wise独立的。存储这样一个函数只需要存储 k 个系数,即 O(k log q) 比特。当 k 为常数时,这仅为 O(log n) 比特,远小于完全随机函数所需的 O(n)

将KMV算法变为现实

现在,我们使用一个 2-wise独立(即成对独立)的哈希函数族(如多项式哈希,k=2)来实现KMV算法。

  1. 选择一个 2-wise独立 的哈希函数 g: {1,...,n} -> {1,...,M},其中 M 是一个很大的数(例如 n^3)。
  2. 定义我们的“哈希值”为 h(i) = g(i) / M。这样 h(i)[0,1) 区间内分母为 M 的有理数,具有有限精度。
  3. 运行KMV算法,但使用 h(i) 代替理想化的均匀随机实数。

需要解决的两个问题:

  1. 哈希碰撞:由于 h 现在只有有限精度,两个不同的元素 ij 可能哈希到同一个值(h(i) = h(j)),这会干扰“精确计数”阶段。我们需要确保碰撞概率很低。通过设置足够大的 M(如 n^3),并利用生日悖论原理,可以证明在最多 n 个不同元素的情况下,发生任何碰撞的概率极低。
  2. 分析的正确性:原始FM算法的分析依赖于 T 个哈希值的完全独立性来计算 min 的分布。现在我们只有 2-wise独立性,无法直接沿用该分析。

新的分析策略(针对KMV)
我们不再直接分析 xk 的分布,而是分析“估计值 k/xk 偏离真实值 T 太远”的概率。这可以转化为以下两个“坏事件”:

  • 估计值过高:这意味着 xk 太小。而 xk 太小意味着至少有 k 个哈希值小于某个阈值 α
  • 估计值过低:这意味着 xk 太大。而 xk 太大意味着少于 k 个哈希值小于某个阈值 β

对于每个坏事件,我们定义指示随机变量 Y_j,表示第 j 个独立元素的哈希值是否小于阈值。那么:

  • Y = Σ Y_j 表示哈希值小于阈值的元素数量。
  • 我们可以计算 E[Y]Var[Y]
    • E[Y] 的计算只用到线性期望,不需要独立性。
    • Var[Y] = E[Y²] - (E[Y])²。计算 E[Y²] 涉及到形如 E[Y_i * Y_j] 的项,这只依赖于两个随机变量 Y_iY_j 的联合分布。由于我们的哈希函数是 2-wise独立的,h(I_i)h(I_j) 是独立的,因此 Y_iY_j 也是独立的。这意味着 2-wise独立性 足以让我们精确计算方差!
  • 一旦知道了 E[Y]Var[Y],我们就可以使用切比雪夫不等式来 bound “Y 偏离其期望值很远”的概率。通过巧妙设置阈值 αβ,我们可以证明这两个坏事件发生的概率都很小,从而保证KMV算法在使用 2-wise独立哈希函数时,仍然能以高概率给出 (1±ε) 近似。

总结

本节课我们一起学习了数据流中独立元素计数的问题。

  1. 我们首先了解了问题的定义和朴素解法。
  2. 接着,我们学习了FM算法,一个基于理想化完全随机哈希函数的简单估计器,并分析了其期望和方差。
  3. 然后,我们通过引入FM+(求平均)和FM++(取中位数)技术,提升了算法的精度和鲁棒性。
  4. 最后,我们探讨了如何移除不切实际的理想化假设。核心是使用k-wise独立哈希函数族(如多项式哈希)作为伪随机源。我们介绍了KMV算法,它通过存储k个最小哈希值来获得更好的方差特性。我们概述了如何利用 2-wise独立性 来分析KMV算法,确保其在有限内存下依然能提供可靠的近似估计。

通过本节课,我们看到了如何将概率工具(期望、方差、集中不等式)与计算理论概念(伪随机性、k-wise独立性)相结合,设计出高效实用的数据流算法。

003:几何采样

在本节课中,我们将要学习一种称为“几何采样”的技术,并将其应用于“不同元素计数”问题。这是一种在许多问题中都会用到的通用技术。我们还将讨论如何证明算法的下界,即如何证明不存在使用少于特定内存量的流算法或草图算法。

概述

几何采样的核心思想是,通过一个精心设计的哈希函数将数据流划分到多个子流中,每个子流对应一个更简单的问题实例(例如,不同元素数量较少)。然后,我们并行运行多个针对简单实例设计的数据结构,并在查询时,巧妙地组合这些子结构的结果,以获得对整个原始问题的近似解。

几何采样的核心思想

我们有一个数据流,其中的更新项来自一个大小为 N 的宇宙 {1, ..., N}

几何采样的主要思想是选取一个哈希函数 h,将宇宙映射到范围 {0, 1, ..., log N} 中,并且满足概率 P(h(i) = j) ≈ 1 / 2^j

一种实现方式是:首先选取另一个将宇宙映射到自身的哈希函数 (假设 N 是2的幂)。然后定义 h(i)h̃(i) 的二进制表示中最低有效位1的索引(从0开始计数)。例如,如果 h̃(i) 的二进制是 ...1000,则 h(i) = 3。这个事件发生的概率恰好是 1/2^4 = 1/16,大致符合 1/2^j 的规律。

这样设计的原因是,我们可能有一个数据结构 D,只能解决某个简化版本的问题(例如,当不同元素数量 F0 不超过某个小值 k 时,可以精确计算 F0)。几何采样允许我们利用这个简单的数据结构来解决更一般的问题。

具体操作是,我们根据哈希函数 h 的值,将原始流 S 划分成 log N + 1 个子流:S_0, S_1, ..., S_{log N},其中 S_j 包含所有满足 h(i) = j 的更新项 i。然后,我们为每个子流 S_j 实例化一个对应的简化版数据结构 D_j

在查询时,算法需要决定使用哪个 D_j 的查询结果,并对其进行调整(例如,乘以一个缩放因子 2^j),以得到对整个流的估计值。这个选择与调整的策略是问题特定的。

应用于不同元素计数

现在,我们具体看看如何将几何采样用于不同元素计数问题。

简化问题:有界不同元素计数

首先,我们定义一个简化问题:我们承诺流中不同元素的数量 F0 至多为 k。在这个承诺下,我们希望精确计算 F0

一个简单的解决方案是:在内存中维护所有出现过的不同元素。由于最多只有 k 个,我们最多需要 k * log N 比特的内存(每个元素名需要 log N 比特)。如果承诺被违反(即不同元素超过 k 个),我们可以检测到并报错。我们将这个简单的数据结构称为“平凡数据结构”。

完整算法框架

我们的目标是估计流中不同元素的总数 F0,并得到一个 (1 ± ε) 的近似值,且不使用太多内存。

以下是算法的高层框架:

  1. 设置参数:令 k = C / ε^2,其中 C 是一个稍后确定的常数。
  2. 初始化数据结构
    • D':一个上述的“平凡数据结构”,参数为 k
    • D_0, D_1, ..., D_{log N}:每个都是上述的“平凡数据结构”,参数也为 k
    • 选择一个哈希函数 h,满足 P(h(i) = j) ≈ 1 / 2^j,且来自一个2-wise独立族。
  3. 处理更新:对于流中每个元素 i
    • 将其插入 D'
    • 计算 j = h(i),然后将 i 插入 D_j
  4. 回答查询
    • 首先查询 D'。如果 D' 报告不同元素数 ≤ k,则直接返回 D' 的精确结果。
    • 否则(即 F0 > k),我们需要利用 D_0, ..., D_{log N}
      • 理想情况下,我们希望找到一个层级 R*,使得 F0 / 2^{R*} 大约是 Θ(1/ε^2)
      • 然后,我们查询 D_{R*},得到估计值 X
      • 最终返回的估计值为 2^{R*} * X

算法原理分析

为什么这样做是有效的?

  • 期望:在层级 R,期望看到的不同元素数量是 E = F0 / 2^R
  • 方差与误差:由于哈希函数的2-wise独立性,实际数量 Y 的方差 Var(Y) ≤ E。根据切比雪夫不等式,Y 偏离其期望的程度大约在 O(√E) 范围内。
  • 选择 R*:我们选择 R* 使得 E ≈ C/ε^2。那么,标准差 √E ≈ √C / ε。相对误差(标准差除以期望)就是 (√C / ε) / (C / ε^2) = ε / √C。通过设置常数 C 足够大,我们可以使这个相对误差小于 ε
  • 缩放:因此,YE 的一个 (1 ± O(ε)) 近似。而 E = F0 / 2^{R*},所以 2^{R*} * Y 就是 F0 的一个 (1 ± O(ε)) 近似。

关键问题:如何找到 R*

我们并不知道 F0,因此无法直接计算哪个 R 能满足 F0 / 2^R ≈ 1/ε^2

解决方案是:先获得一个 F0 的常数因子近似

我们可以运行另一个并行的、更粗糙的几何采样过程来获得这个常数近似。

  • 常数因子近似算法:这个算法结构与主算法类似,但更简单。
    • 我们维护另一组数据结构 D''_0, ..., D''_{log N},但这里的“平凡数据结构”参数 k 设置为 1(即只能记录0或1个不同元素)。
    • 处理流时,我们也用同一个哈希函数 h 将元素送入对应的 D''_j
    • 在查询时,我们找到最大的 R,使得 D''_R.query() ≠ 0(即该层级至少看到了一个元素)。
    • 返回 2^R 作为 F0 的估计值 F0_tilde

原理F0_tilde 大概率落在 [F0/16, 16F0] 范围内,即它是一个常数因子近似。证明思路是,远低于 log F0 的层级几乎不可能有元素,远高于 log F0 的层级几乎肯定有元素,因此找到的 R 会集中在 log F0 附近。

有了 F0_tilde 后,我们就可以计算一个候选的 R* = log2( ε^2 * F0_tilde )。由于 F0_tildeF0 的常数倍近似,这个 R* 也能满足 F0 / 2^{R*}Θ(1/ε^2)

最终算法流程总结

  1. 并行运行两个几何采样过程:
    • 过程A(常数近似):使用 k=1 的平凡数据结构,用于估计 F0_tilde
    • 过程B(精细近似):使用 k = C/ε^2 的平凡数据结构,并附带一个 D' 处理小 F0 的情况。
  2. 对于每个流元素 i,同时送入过程A和过程B。
  3. 查询时:
    • 检查过程B中的 D'。如果 F0 ≤ k,直接返回精确结果。
    • 否则,用过程A得到的 F0_tilde 计算 R*
    • 查询过程B中的 D_{R*},得到 X
    • 返回 2^{R*} * X 作为最终估计。

空间复杂度与优化

初始分析显示,空间复杂度约为 O( (1/ε^2) * log^2 N ) 比特。这可以优化:

  1. 优化平凡数据结构:存储 k 个元素不需要 k log N 比特,可以通过完美哈希用 O(k log k) 比特实现。
  2. 优化层级存储:不需要同时维护所有 log ND_j。可以只维护一个数据结构,当它“满”时,丢弃哈希值最低位为0的元素,剩下的视为下一层的数据,如此递推。对于常数近似过程,甚至只需要记录遇到的最大 h(i) 值。
  3. 最终空间:经过优化,空间复杂度可降至约 O( (1/ε^2) * (log(1/ε) + log log N) + log N ) 比特。这几乎是最优的。

与HyperLogLog的关联

实践中著名的 HyperLogLog 算法也基于几何采样的思想。其核心结构可以看作一个 (log N) × m 的位矩阵(m ≈ 1/ε^2):

  • 行对应几何采样的层级。
  • 列由另一个哈希函数决定。
  • 每个元素根据哈希函数确定一个行和一个列,并将该位置设为1。
  • 查询时,通过分析每列第一个被置1的行号(占 log log N 比特),并使用调和平均数等技巧进行估计,最终得到 F0
  • 其空间约为 O( (1/ε^2) * log log N ) 比特。

总结

本节课我们一起学习了几何采样这一强大技术。我们首先了解了其核心思想:利用一个概率呈几何级数衰减的哈希函数将流划分,从而创造出许多“简化版”的子问题。然后,我们深入探讨了如何将这一技术应用于不同元素计数问题,构建了一个完整的 (1±ε) 近似算法。我们分析了算法的原理,包括如何通过另一个并行的采样过程获得常数因子近似以确定合适的采样层级,并讨论了算法的空间复杂度及其优化路径。最后,我们指出了几何采样思想与实践中广泛应用的高效算法(如HyperLogLog)之间的紧密联系。

004:分位数

在本节课中,我们将要学习数据流中的分位数问题。我们将介绍问题的定义、几种确定性算法(包括Q-digest和MRL算法)以及随机化算法(如KLL草图)的核心思想。我们还将探讨相关的空间下界。内容将尽可能简单直白,以便初学者理解。

问题定义

我们有一个由可比较项组成的数据流。例如,这些项可以是来自某个范围 [1, U] 的整数,或者是字符串。流的长度为 n

我们希望能够回答以下两种查询:

  1. Rank(x):查询项 x 在所有已见流项排序后的顺序中的排名。
  2. Select(r):查询在排序后列表中排名为 r 的项。

由于内存限制,我们无法存储整个排序列表。因此,我们接受近似解:

  • Rank(x):返回的排名误差在 ± εn 以内,其中 ε 是预先给定的参数。
  • Select(r):返回的项,其真实排名在 r ± εn 以内。

确定性算法:Q-digest

上一节我们定义了分位数问题,本节中我们来看看第一个确定性算法:Q-digest。该算法假设流中的项是 1U 之间的整数。

算法原理

Q-digest 的核心思想是构建一棵覆盖整个值域 [1, U] 的完美二叉树。树的高度为 log₂ U。每个树节点 v 关联一个计数器 c(v),用于记录“归属于”该节点值域范围内的流元素数量。关键在于,每个流元素只被一个节点计数。

初始时,元素在其对应的叶子节点计数。为了节省空间(即减少非零计数器的数量),我们有一个合并规则:当某个节点 x、其父节点 p(x) 及其兄弟节点 s(x) 的计数器之和不超过阈值 εn / log₂ U 时,我们将子节点的计数合并到父节点。

以下是合并操作的伪代码描述:

if c(x) + c(p(x)) + c(s(x)) <= εn / log U:
    c(p(x)) += c(x) + c(s(x))
    c(x) = 0
    c(s(x)) = 0

空间与误差分析

我们可以证明,在任何时刻,非零计数器的节点数 K 最多为 3 log U / ε。因此,我们只需要存储这些节点的ID和计数器值。

对于排名查询 Rank(x),误差来源于从根到 x 对应叶子的路径上所有节点的计数器。因为这些节点中的元素可能来自 x 的左侧或右侧,我们无法区分。路径上的节点数最多为 log U,每个节点的计数值最多为 εn / log U,因此总误差最多为 εn

确定性算法:MRL(比较排序法)

上一节我们介绍了基于值域树的Q-digest算法,本节中我们来看看一种基于比较和压缩的算法:MRL算法。该算法是比较排序的,适用于任何可比较项。

算法原理

算法维护 L = log₂(n/k) 个“压缩器”,每个压缩器可以存储最多 k 个项。我们将它们编号为 0L-1

流程如下:

  1. 新流元素总是插入到第 0 级压缩器。
  2. 当第 j 级压缩器存满 k 个项时,我们将其中的项排序,然后提升所有奇数索引项到第 j+1 级压缩器,接着清空第 j 级压缩器。

查询与误差

要回答 Rank(x) 查询,我们计算所有压缩器中小于等于 x 的项,但每项根据其所在压缩器级别 j 具有权重 2^j(因为它代表了 2^j 个原始项)。

误差产生于压缩过程。每次在第 j 级发生的压缩,最多可能引入 2^(j-1) 的排名误差。通过分析,所有压缩引入的总误差上限为 (n/k) * log(n/k)。为了使这个误差小于 εn,我们需要设置 k = O( (1/ε) * log(εn) )

因此,总空间复杂度为 k * L = O( (1/ε) * log²(εn) ) 个单词。

随机化算法:引入随机性

MRL算法是确定性的。我们可以通过引入随机性来改进空间复杂度。核心改进在于压缩步骤:当压缩器满时,排序后我们随机地选择提升所有奇数索引项所有偶数索引项。

误差分析

此时,每次压缩引入的误差不再是固定的 2^(j-1),而是随机的 ±2^(j-1)(有时误差为0)。总误差 E 是一个随机变量,形式为许多带随机符号 σ 的项之和:
E = Σ_j Σ_r η_{j,r} * σ_{j,r} * 2^(j-1)

我们可以利用 Khinchin不等式 来约束这个随机和。该不等式指出,对于独立的随机 ±1 变量 σ_i 和实数 x_i,有:
Pr[ |Σ σ_i x_i| > λ ] ≤ 2 * exp( -λ² / (2 * Σ x_i²) )

在我们的设定中,Σ x_i² 可以证明至多为 O(n² / k²)。为了使误差超过 εn 的概率小于 δ,我们需要设置 k = O( (1/ε) * √log(1/δ) )

结合级别数 L = log(n/k),总空间为 O( (1/ε) * √log(1/δ) * log(εn) )

进阶优化:KLL草图

为了达到更优的空间复杂度 O( (1/ε) * log log(1/δ) ),KLL草图采用了多重优化:

  1. 几何增长的压缩器容量:不再让所有压缩器容量相同,而是从底层开始容量按几何级数增长(例如,比例为 2/3)。这改变了空间构成,从乘以 log n 因子变为加上 log n 因子。
  2. 采样替代小型压缩器:最底层的多个小型压缩器(容量为2)的功能等价于一个采样器,可以从连续块中均匀采样。我们可以用常数空间的采样器来实现这一部分。
  3. 混合确定性算法:将最高层的几个压缩器替换为一个最优的确定性分位数算法(如Greenwald-Khanna算法)。
  4. 误差和分离分析:将总误差随机和分为两部分:对低层级部分使用Khinchin不等式进行紧致界分析;对高层级部分使用简单的确定性界(即忽略随机抵消)。通过精细设置参数,最终得到 O(1/ε * log log(1/δ)) 的空间上界。

下界与总结

对于确定性的、基于比较的算法,存在一个空间下界为 Ω( (1/ε) * log(εn) ) 个单词。Greenwald-Khanna算法匹配了这个下界,但其分析和实现较为复杂。

本节课中我们一起学习了数据流分位数问题的近似解法。我们从简单的确定性算法Q-digest和MRL开始,理解了压缩和合并的基本思想。接着,我们看到了如何通过引入随机性(随机提升策略)并利用概率工具(Khinchin不等式)来分析误差,从而改进算法。最后,我们概述了目前最先进的KLL草图算法所采用的一系列优化技术,使其能达到近乎最优的空间复杂度。这些算法在性能监控等实际场景中有着广泛的应用。

005:基于压缩论证的下界证明

在本节课中,我们将学习一种证明数据流算法空间下界的核心技巧——基于压缩的论证。我们将通过两个经典问题——元素唯一性问题和分位数问题——来具体展示这种技巧的应用。课程内容将循序渐进,从简单的确定性精确算法下界开始,逐步深入到更复杂的确定性和随机性近似算法的下界证明。

概述:压缩论证的核心思想

基于压缩的论证是一种证明下界的通用技术。其核心思路如下:

假设存在一个针对某个问题(如元素唯一性)的素描算法 A,它使用 S 比特的内存。那么,我们可以构造一个单射 F,将某个大集合 Q 映射到一个有界集合(例如所有长度为 S 的比特串集合 {0,1}^S)。

由于单射意味着集合 Q 的大小不能超过目标集合的大小,即 |Q| ≤ 2^S,因此我们可以推导出 S ≥ log₂(|Q|)

证明的艺术在于如何巧妙地构造这个单射,以及如何定义集合 Q。我们将看到,对于不同的问题和算法类型,构造方法会有所不同。

热身:确定性精确算法的下界

上一节我们介绍了压缩论证的基本框架。本节中,我们首先来看一个简单的热身例子:针对元素唯一性问题的确定性精确算法的下界。

定理:假设存在一个确定性且精确的元素唯一性算法 A,在长度为 Ω(n) 的流上使用 S 比特内存,那么 S ≥ n

证明:我们将证明,如果存在这样的算法 A,那么我们可以构造一个从大小为 2^n 的集合到大小为 2^S 的集合的单射。由于 2^n ≤ 2^S 必须成立,因此 S ≥ n

以下是构造单射(编码函数 ENC)和解码函数 DEC 的方法:

编码过程 (ENC)
给定一个比特向量 x ∈ {0,1}^n,我们将其解释为一个集合 X = {i | x_i = 1}

  1. 创建一个包含 X 中所有元素的流。
  2. 在流上运行算法 A。
  3. 编码输出就是算法 A 完成后的内存内容(一个 S 比特的字符串)。

解码过程 (DEC)
给定一个 S 比特的内存状态 M,我们需要恢复原始的比特向量 x

  1. 将算法 A 的内存初始化为 M
  2. 查询算法 A,得到当前流中不同元素的数量 s
  3. 初始化返回值 x 为全零向量。
  4. 对于 i 从 1 到 n
    • 向算法 A 更新元素 i
    • 再次查询算法 A,得到新的结果 r
    • 如果 r == s,则意味着 i 在原始流中已经出现过(否则不同元素数应增加),因此设置 x_i = 1
    • 否则,设置 x_i = 0,并更新 s = r
  5. 返回 x

由于算法 A 是确定且精确的,上述解码过程能完美地恢复出原始的 x。这意味着编码函数 ENC 是一个单射。因为 ENC2^n 个可能的 x 映射到 2^S 个可能的内存状态,所以必须有 2^n ≤ 2^S,即 S ≥ n

这个简单的论证表明,对于元素唯一性问题,任何确定性精确算法都无法实现亚线性空间。

关键工具:纠错码简介

在证明更复杂的下界(如确定性近似算法的下界)之前,我们需要引入一个关键工具:纠错码。

定义:一个参数为 (Q, L, D) 的码 C 是集合 {1, ..., Q}^L 的一个子集。其中 Q 是字母表大小,L 是块长度(码字长度),D 是码的最小汉明距离(任意两个不同码字之间不同位置的最小数量)。

性质:如果一个码的最小距离为 D,那么任意一个码字在经历少于 D/2 个位置的错误(任意修改)后,得到的字符串仍然可以唯一地解码回原始的码字。这是因为以每个码字为中心、半径为 D/2 的汉明球是互不相交的。

存在性:通过概率方法可以证明,对于给定的 Qn,存在一个码 C,其块长度 L = O(Q log n),包含 n 个码字,且具有相对距离至少为 1 - O(1/Q)

基于这样的码,我们可以推导出一个关于集合的推论,它将用于后续的下界证明:

推论:给定 QL,存在一个集合族 B_{Q,L},其中每个集合是 {1, ..., n} 的子集(n = Q*L),满足以下条件:

  1. 每个集合的大小恰好为 L
  2. 任意两个不同集合的交集大小至多为 6L/Q
  3. 集合族的大小至少为 exp(Ω(L/Q))

这个推论可以通过将码 C 中的每个符号(1 到 Q)扩展为一个长度为 Q 的独热编码(one-hot encoding)比特向量来构造。

确定性近似算法的下界

上一节我们引入了纠错码这个强大工具。现在,我们利用它来证明针对元素唯一性问题的确定性近似算法的下界。

定理:假设算法 A 是一个确定性算法,对于元素唯一性问题,它总是输出估计值 ,满足 T ≤ T̃ ≤ 1.1 * T,并使用 S 比特内存,那么 S = Ω(n)

证明思路:我们将证明,这样的算法 A 蕴含了一个从推论中构造的大集合族 B_{Q,L}(取 Q=100, L=n/100)到小集合 {0,1}^S 的单射。由于 |B_{Q,L}| = exp(Ω(n)),这将迫使 S = Ω(n)

编码 (ENC):与热身例子类似。对于一个集合 F ∈ B_{Q,L},创建包含 F 所有元素的流,运行算法 A,取其最终内存状态作为编码。

解码 (DEC):解码器接收内存状态 M,目标是恢复 F。它采用如下“猜测-验证”策略:

  1. 对于 B_{Q,L} 中的每一个候选集合 Q
    • 将算法 A 的内存重置为 M(此时 A“认为”它刚处理完流 F)。
    • 将候选集合 Q 中的所有元素依次输入算法 A 进行更新。
    • 查询算法 A,得到新的估计值
    • 判断:如果 T̃ ≤ 1.1 * L,则输出 Q 作为解码结果。

解码正确性分析

  • 如果 Q == F,那么算法 A 第二次看到的流仍然是 F(因为 F ∪ Q = F)。不同元素数仍为 L,因此其输出 满足 L ≤ T̃ ≤ 1.1L。解码器会接受这个 Q
  • 如果 Q ≠ F,那么算法 A 第二次看到的流是 F ∪ Q。根据推论,|F ∩ Q| ≤ 6L/Q = 0.06L。因此,|F ∪ Q| = |F| + |Q| - |F ∩ Q| ≥ L + L - 0.06L = 1.94L。算法 A 的输出 至少为 1.94L,这大于 1.1L。解码器会拒绝这个 Q

由于 B_{Q,L} 中任意两个集合的交集都很小,保证了错误候选集 Q 会导致显著更大的并集。因此,解码器能够唯一且正确地识别出原始的集合 F。这证明了 ENC 是一个单射,从而完成了下界证明。

随机性精确算法的下界

最后,我们探讨随机性算法是否能够避免线性的空间下界。答案是否定的,即使允许随机性和小概率错误,精确计算元素唯一性仍然需要 Ω(n) 空间。

定理:假设算法 A 是一个随机化算法,以至少 2/3 的概率精确输出不同元素的数量,并使用 S 比特内存,那么 S = Ω(n)

证明思路:核心思想是通过多次独立运行并取中位数,先将算法 A 强化为一个错误率极低(如 10^{-6})的算法 A‘。然后,我们为算法 A’ 的每个固定随机种子 R 定义其对应的确定性编码/解码函数。虽然对于单个集合,解码可能失败,但我们可以证明存在一个“好”的随机种子 R*,使得基于 R* 的确定性编解码函数能够正确处理一个相当大的集合族(大小为 |B_{Q,L}|/2)。这就构造出了从一个大小为 exp(Ω(n)) 的集合到 {0,1}^S 的单射,从而证得下界。

关键步骤

  1. 错误率放大:通过并行运行 O(1) 份 A 的副本并输出中位数,得到新算法 A‘,其失败概率 ≤ 10^{-6},空间使用为 O(S)。
  2. 定义随机编解码:对于算法 A‘,定义依赖于随机串 R 的编码函数 ENC_R 和解码函数 DEC_R(解码过程类似确定性下界中的“逐元素探测”)。
  3. 分析期望性能:对于固定的集合 F ∈ B_{Q,L},由于 A‘ 的错误率极低,DEC_R(ENC_R(F)) 输出错误结果的概率 < 1/2(通过马尔可夫不等式,并利用纠错码属性处理少量错误)。
  4. 寻找确定性子函数:定义指示变量 Z_F(R),当 DEC_R(ENC_R(F)) 正确时为 1。我们有 E_R[Σ_F Z_F(R)] > |B_{Q,L}|/2。因此,存在一个固定的“好”随机种子 R*,使得 Σ_F Z_F(R*) > |B_{Q,L}|/2。令 B 为被 R* 正确处理的那些 F 的集合,则 |B| ≥ |B_{Q,L}|/2
  5. 完成单射构造:固定 R*,函数 ENC_{R*} 就是从大集合 B{0,1}^{O(S)} 的一个单射。因此,O(S) ≥ log|B| = Ω(n),即 S = Ω(n)

这个论证巧妙地通过概率方法,从一个随机化算法中“提取”出了一个对大量输入都有效的确定性单射,从而将随机化算法归约到确定性情况来证明下界。

总结与展望

本节课中我们一起学习了基于压缩论证来证明数据流算法空间下界的强大技巧。

  • 我们首先看到了一个简单的例子,证明了元素唯一性问题的任何确定性精确算法都需要 Ω(n) 空间。
  • 接着,我们引入了纠错码作为构造具有“大距离”的集合族的工具。
  • 利用纠错码,我们证明了即使对于常数因子近似(如 1.1 倍)的确定性算法,求解元素唯一性问题仍然需要 Ω(n) 空间。
  • 最后,我们通过更复杂的概率方法论证,证明了即使是允许随机性和小概率错误的精确算法,也逃不过 Ω(n) 的空间下界。

这些结果深刻地说明了,为了在亚线性空间内解决元素唯一性问题,随机化和近似两者都是必不可少的。这正好对应了我们在课程中学到的 FM 算法、HyperLogLog 等算法的特性。

基于压缩的论证技巧并不仅限于元素唯一性问题。在课程笔记中,它还被应用于证明分位数问题等其他问题的下界。这种通过算法构造编码,进而利用信息论界限来证明空间下界的思路,是数据流算法理论中的一个核心工具。

006:通过通信复杂度证明下界

在本节课中,我们将结束关于下界的讨论,并开始介绍流式处理和素描算法中的“翻转式”模型。我们将重点讨论该模型下的素描技术。

概述

上一节我们介绍了通过压缩论证来证明下界的方法。本节中,我们将探讨另一种方法:通信复杂度。我们将看到如何利用通信复杂度中的困难问题,来证明流式算法所需内存的下界。

通信复杂度简介

通信复杂度是一个计算模型。在这个模型中,通常有两个参与方,称为Alice和Bob。Alice有一个输入X,Bob有一个输入Y。他们的目标是计算某个关于XY的函数F(X, Y)

模型的核心是:Alice和Bob可以通过互相发送消息来进行交流。最后收到消息的一方输出答案。目标是在保证正确计算函数F的前提下,最小化通信过程中传输的总比特数

通信复杂度有几个重要的变体模型:

  • 确定性通信复杂度:协议是确定性的。
  • 随机性通信复杂度
    • 公开随机性:Alice和Bob共享一个公开的随机字符串。
    • 私有随机性:Alice和Bob各自私下抛硬币生成随机性。
  • 分布通信复杂度:输入XY是从某个分布中随机抽取的,协议可以是确定性的,目标是高概率成功。

通信复杂度与流式下界的联系

通信复杂度与流式算法下界的基本联系思想是构造一个归约

核心思想:如果计算某个函数F在通信模型中很困难(需要大量通信),并且我们可以利用一个解决流式问题P的算法A,以黑盒方式构造一个解决F的通信协议,那么P也必定是困难的(需要大量内存)。

具体来说,归约通常这样操作:

  1. Alice将她的输入视为流的一部分,在本地运行流式算法A
  2. Alice将算法A内存状态作为消息发送给Bob。
  3. Bob接收到内存状态后,将其作为初始状态,继续运行算法A处理他自己的输入部分。
  4. Bob最终查询算法A,得到答案,从而推断出F(X, Y)的结果。

在这个协议中,传输的比特数就等于算法A的内存使用量。因此,如果F需要Ω(S)的通信量,那么算法A就至少需要Ω(S)的内存。

示例:利用相等函数证明精确判重下界

让我们看一个具体例子:利用相等函数的下界来证明精确计算不同元素数量需要大量内存。

相等函数 EQ(X, Y):当X = Y时输出1,否则输出0。已知结论:确定性通信复杂度 CC(EQ) ≥ n,其中n是输入字符串的长度。

定理:任何确定性的、能精确计算不同元素数量的流式算法A,至少需要使用n比特内存(假设全集大小为n)。

证明(通过归约):

  1. Alice有一个集合X(表示为n位指示向量),Bob有一个集合Y
  2. Alice在本地运行算法A,处理X中的所有元素。
  3. Alice将算法A内存状态发送给Bob。这消耗了Space(A)比特。
  4. Bob首先检查|X|是否等于|Y|(Alice可以随内存一起发送|X|,这只需log n比特,不影响主导项)。如果不相等,则X ≠ Y,输出False
  5. 如果大小相等,Bob将接收到的内存状态作为初始状态,继续运行算法A,处理Y中的所有元素。
  6. Bob查询算法A,得到当前不同元素数量T
    • 如果X = Y,则流中元素全部来自XT = |X| = |Y|
    • 如果X ≠ Y|X| = |Y|,则Y中必然存在至少一个不在X中的元素,因此T > |Y|
  7. 通过比较T|Y|,Bob可以判断XY是否相等。

这样,我们就用算法A构造了一个解决EQ问题的通信协议,其通信成本为Space(A)。由于CC(EQ) ≥ n,因此Space(A) ≥ n

随机算法的下界与私有随机性

上面的例子针对确定性算法。我们实际使用的流式算法(如不同元素估计)通常是随机化的。能否为随机近似算法证明下界呢?答案是肯定的,这需要用到私有随机性模型下的通信复杂度。

对于相等函数EQ

  • 公开随机性、单向通信模型中,只需要O(1)比特即可高概率判断相等性(例如,通过随机向量点积)。
  • 但在私有随机性、单向通信模型中,其通信复杂度是Ω(log n)。这是一个指数级的差距。

关键点:在流式算法中,算法使用的随机哈希函数等,其描述是需要存储在内存中的,这类似于私有随机性模型——随机比特不是“免费公开”的,需要成本。因此,为私有随机性模型证明的通信下界,可以对应到随机流式算法的内存下界。

一个关键难题:集合不相交性

一个在通信复杂度中众所周知的难题是集合不相交性问题。

  • Alice和Bob分别有集合ST(均为{1,...,n}的子集)。
  • 函数DISJ(S, T) = 1 当且仅当 S ∩ T = ∅
  • 定理:即使允许使用公开随机性和多轮交互,DISJ的随机通信复杂度也是Ω(n)

我们可以利用这个结论证明另一个流式问题的下界。

定理:在数据流中,以至少2/3的概率,输出L∞范数(即最大频率)的1.1倍近似值,需要Ω(n)比特内存。

证明(归约):

  1. Alice有集合S,Bob有集合T
  2. Alice运行L∞近似算法A,处理S中元素,然后将内存状态发送给Bob。
  3. Bob继续运行算法A,处理T中元素,然后查询得到估计值Z
  4. 情况分析:
    • 如果ST不相交,则流中任何元素最多出现1次,L∞范数 ≤ 1。
    • 如果ST相交,设a ∈ S ∩ T,则a出现了2次,L∞范数 = 2。
  5. 由于算法A提供1.1倍近似,它可以可靠地区分≤ 1.1≥ 2/1.1 ≈ 1.818,从而判断ST是否相交。

因此,算法A的内存大小必须至少是DISJ问题的通信下界,即Ω(n)

扩展到多玩家与Lp范数下界

通过考虑t个玩家的不相交性问题变体(t-玩家唯一相交性问题),我们可以得到依赖于近似比α的下界。

结论:在数据流中,以常数概率计算L∞范数的α倍近似,需要Ω(n / α^2)比特内存。这意味着对于常数近似比(如α=1.1),内存需求是线性的Ω(n),没有非平凡的解。

进一步,我们可以将这个结论推广到Lp范数(p > 2)的估计上。

定理:对于p > 2,在数据流中以常数概率计算Lp范数的(1+ε)倍近似,需要Ω(n^{1 - 2/p})比特内存。

证明思路:通过将t-玩家不相交性问题归约到Lp范数估计。选择合适的玩家数t ≈ n^{1/p},使得在不相交情况下Lp^p ≤ n,而在唯一相交情况下Lp^p ≥ t^p ≈ n,从而产生常数因子差距。再利用t-玩家问题的通信下界Ω(n/t^2),即可得到内存下界Ω(n^{1 - 2/p})。这个下界是紧的,存在匹配的算法。

证明一个通信下界:索引问题与汉明距离间隙问题

最后,我们简要介绍如何直接证明一个通信下界,并将其链接到不同元素估计问题。我们遵循以下归约链:
索引问题间隙汉明距离问题不同元素(1+ε)近似问题

索引问题

  • Alice有一个n比特字符串x
  • Bob有一个索引i ∈ [n]
  • 目标:Bob想知道x_i的值。只允许Alice向Bob发送一次消息。
  • 定理:在公开随机性、单向通信模型中,以失败概率δ < 1/2求解索引问题,需要至少(1 - H_2(δ)) * n比特通信,即Ω(n)

间隙汉明距离问题

  • Alice和Bob各有n比特字符串ab
  • 承诺两种情况之一:
    1. ab的汉明距离 ≥ n/2 + √n
    2. ab的汉明距离 ≤ n/2 - √n
  • 目标:区分属于哪种情况。
  • 定理:该问题的随机通信复杂度也是Ω(n)

归约:索引 → 间隙汉明距离

  1. 设要解决索引问题实例:Alice有x∈{0,1}^n,Bob有索引i
  2. 利用公开随机性,生成N个随机±1向量R1, ..., RNN足够大)。
  3. Alice将x转换为x'0→+1, 1→-1),并计算字符串A,其中A_j = sign(〈x', Rj〉)
  4. Bob计算字符串B,其中B_j = (Rj)_i(即第i个分量)。
  5. 可以分析,AB的汉明距离的期望值取决于x_i是0还是1:
    • x_i = 0,期望汉明距离 ≈ N/2 + c√N
    • x_i = 1,期望汉明距离 ≈ N/2 - c√N
      (其中c为某个常数)
  6. 通过设置N = Θ(n),这个期望差距约为√N = Θ(√n)。利用集中不等式,可以保证以高概率,汉明距离落在n/2 ± Θ(√n)的不同区间。
  7. 因此,任何能解决间隙汉明距离问题(间隙为Θ(√n))的协议,都可以用来解决索引问题。由于索引问题需要Ω(n)通信,间隙汉明距离问题也需要Ω(n)通信。

归约:间隙汉明距离 → 不同元素(1+ε)近似

  1. 设要解决间隙汉明距离实例,字符串长度为m = 1/ε^2
  2. Alice和Bob将各自的m比特字符串ab视为集合的指示向量。
  3. Alice运行不同元素(1+ε)近似算法A,处理她集合中的元素,并将内存状态以及她集合的大小发送给Bob。
  4. Bob继续运行算法A,处理他集合中的元素,然后查询得到不同元素数量的估计值。
  5. 可以证明,(1+ε)近似不同元素数量,允许我们以±εm的误差估计ab对应集合的对称差大小,进而推断其汉明距离(因为Ham(a,b) = |sym_diff(S_a, S_b)|)。
  6. 由于间隙汉明距离问题的间隙是√m = 1/ε,而εm = 1/ε,因此该估计足以区分两种承诺情况。
  7. 通信成本主要是算法A的内存M,加上O(log(1/ε))比特。由于间隙汉明距离问题需要Ω(m) = Ω(1/ε^2)通信,因此M必须为Ω(1/ε^2)

这就得到了一个重要的下界:任何在数据流中提供(1+ε)近似的不同元素估计算法,至少需要Ω(1/ε^2)比特内存。这个下界是紧的,存在算法(如HyperLogLog的变种)使用O(1/ε^2 + log n)比特内存。

总结

本节课我们一起学习了如何利用通信复杂度来证明数据流算法的下界。

  1. 我们介绍了通信复杂度的基本模型及其变体(确定性、随机性公开/私有硬币)。
  2. 我们理解了通信复杂度与流式下界之间的核心归约思想:将流式算法用作通信协议的黑盒子。
  3. 我们通过相等函数集合不相交性的例子,展示了如何为精确算法和近似算法证明内存下界。
  4. 我们探讨了将多玩家不相交性问题下界应用于L∞Lp范数估计,得到了依赖于近似比α和范数指数p的下界。
  5. 最后,我们概述了通过索引问题间隙汉明距离问题不同元素估计问题的归约链,证明了不同元素(1+ε)近似需要Ω(1/ε^2)比特内存的关键下界。

这些下界结果告诉我们,对于许多重要的流式统计问题,我们所看到的高效素描算法在内存使用上本质上是最优的。

007:线性素描、旋转门流模型与频繁项

在本节课中,我们将学习一种设计流算法的通用技术——线性素描。我们将探讨其在旋转门流模型中的应用,并重点解决频繁项(或称“大热门”)问题。课程内容将涵盖线性素描的基本思想、旋转门流模型的几种变体,以及如何使用线性素描技术(特别是Count-Min和Count Sketch算法)来高效解决点查询和频繁项发现问题。

线性素描与旋转门流模型概述

线性素描是一种设计流算法和素描算法的通用技术,适用于数据流处理以及分布式应用等场景。其核心思想是将数据表示为一个高维向量 x,然后选取一个行数远小于列数的矩阵 Π,在内存中存储 y = Πx。由于 y 的维度远小于 x,因此可以节省大量内存。关键在于,矩阵 Π 虽然庞大,但通常可以简洁地表示(例如,通过确定性算法或使用少量种子的随机哈希函数生成),而无需显式存储整个矩阵。

在流处理中,这与旋转门流模型密切相关。在该模型中,我们需要维护一个初始为零的向量 x,并处理形如 update(i, Δ) 的更新操作,该操作将 xᵢ 增加 Δ。向量 x 可以看作是一个直方图,记录了流中各个元素的频率或计数。

旋转门流模型有几个重要的子模型:

  1. 仅插入模型Δ 始终为 +1。这是我们目前主要讨论的模型。
  2. 严格旋转门模型Δ 可为任意实数,但保证 x 的每个分量始终非负。这模拟了允许插入和删除,但不会删除不存在的项的场景。
  3. 通用旋转门模型Δ 可为任意实数,x 的分量可能为负。这适用于需要计算差异的场景,例如检测两天数据流之间的变化。

在线性素描框架下,处理更新 update(i, Δ) 非常高效:它等价于将素描 y 更新为 y + Δ * Πᵢ,其中 Πᵢ 是矩阵 Π 的第 i 列。研究表明,在旋转门流模型中,线性素描在某种意义上具有“普适性”,即任何高效的旋转门流算法本质上都可以表示为某种线性素描。

频繁项(大热门)问题

频繁项问题,即近似地找出流中出现频率最高的前 k 个元素。然而,精确求解前 k 个频繁项需要 Ω(n) 的内存,这可以通过从多集合相交问题的通信复杂度下界推导得出。因此,我们转向一个可近似求解的问题:L₁ 频繁项

定义:给定参数 k,如果元素 i 的频率 xᵢ 满足 xᵢ > (1/k) * ||x||₁,则称其为 k-频繁(或称 k-大热门)。这里 ||x||₁ 是向量 x 的 L₁ 范数(即所有频率之和)。满足此条件的元素最多有 k-1 个。

目标:设计一个算法,当查询时,能输出一个大小至多为 O(k) 的列表 L,并保证所有 k-频繁 项都包含在 L 中。

一个相关的问题是 L₁ 点查询:当查询索引 i 时,返回一个值 x̃ᵢ,满足 |x̃ᵢ - xᵢ| ≤ (1/k) * ||x||₁

声明:如果能以空间 S 和每个查询失败概率至多 δ/n 解决 3k-L₁ 点查询问题,那么就可以构造一个算法,以空间 S + O(k) 和总失败概率至多 δ 来解决 k-L₁ 频繁项问题。

证明思路:使用点查询数据结构作为黑盒。查询时,对全集中的每个元素 i 都进行点查询,得到估计值 x̃ᵢ,然后输出估计值最大的 3k 个元素作为列表 L。可以证明,在所有的点查询都成功的条件下,所有真正的 k-频繁 项都会包含在这个列表 L 中。

Count-Min Sketch 算法

接下来,我们介绍第一个能在旋转门流模型中解决上述问题的线性素描算法——Count-Min Sketch,它适用于严格旋转门模型(x 的分量非负)。

数据结构

  • 维护一个 LB 列的计数器数组。
  • 每一行 r 关联一个独立的二元哈希函数 h_r,将元素映射到 {1, 2, ..., B}
  • 设置 B = 2kL = ⌈log₂(1/δ)⌉,其中 δ 是期望的失败概率。

更新操作 update(i, Δ)
对于每一行 r,计算 j = h_r(i),然后将第 r 行第 j 列的计数器增加 Δ

查询操作 query(i)
对于元素 i,查看它在每一行对应的计数器:c_r = C[r, h_r(i)]。返回所有 c_r 中的最小值作为 xᵢ 的估计值 x̃ᵢ

原理分析
对于某一行 r,计数器 c_r 的值是 xᵢ 加上所有其他哈希到同一位置的元素 jxⱼ 之和(即噪声)。由于在严格旋转门模型中所有 xⱼ ≥ 0,因此 c_r ≥ xᵢ。噪声的期望是 ∑_{j≠i} xⱼ * Pr[h_r(i)=h_r(j)] = (1/B) * ∑_{j≠i} xⱼ ≤ ||x||₁ / B。根据马尔可夫不等式,噪声超过 2||x||₁ / B 的概率小于 1/2。由于 B=2k,这意味着噪声超过 (1/k)*||x||₁ 的概率小于 1/2。

因为所有 L 行的哈希函数是独立的,所以所有行的噪声都超过阈值的概率小于 (1/2)^L ≤ δ。因此,以至少 1-δ 的概率,最小值估计器返回的 x̃ᵢ 满足 xᵢ ≤ x̃ᵢ ≤ xᵢ + (1/k)*||x||₁

空间复杂度O(k log(1/δ)) 个字。

扩展到通用旋转门模型

在通用旋转门模型中,x 的分量可能为负,此时最小值估计器可能严重低估。我们需要修改算法:

  1. 设置 B = 3k
  2. 查询时,不再取最小值,而是取所有 L 个计数器值的中位数作为估计值 x̃ᵢ

分析:对于每一行,噪声的绝对值期望仍为 ≤ ||x||₁ / B。根据马尔可夫不等式,噪声绝对值超过 3||x||₁ / B = (1/k)*||x||₁ 的概率小于 1/3。因此,每行给出“好”估计(误差在阈值内)的概率至少为 2/3。通过切尔诺夫界,可以证明当 L = O(log(1/δ)) 时,超过半数的行给出好估计的概率很高,从而中位数估计器的误差能以高概率控制在要求范围内。

这个改进后的算法通常被称为 Count Sketch

提升频繁项查询效率:Dyadic Trick

前面提到的通过点查询解决频繁项的方法,查询时需要扫描整个全集,耗时 O(n),即使输出列表大小仅为 O(k)。这在实践中不可接受。Dyadic Trick 可以显著提升查询效率。

基本思想

  1. 在概念上为全集构建一棵完全二叉树,树高为 log n
  2. 为树的每一层维护一个独立的 Count-Min Sketch(或 Count Sketch)数据结构,用于对该层聚合的向量(每个内部节点是其子树叶子节点之和)进行点查询。
  3. 更新一个叶子节点 i 时,沿着从叶子到根的路径,更新路径上每个节点在对应层 sketch 中的值。
  4. 查询频繁项时,从根节点开始进行广度优先搜索。在每一层,只对那些点查询估计值超过阈值的节点,继续查询其子节点。由于在严格旋转门模型中,一个叶子节点是频繁项意味着其所有祖先都是频繁项,因此搜索空间被有效限制。

效果:查询时间从 O(n) 降为 O(k log n),而空间和更新时间仅增加一个 log n 因子。

尾部保证与 Lp 频繁项

尾部保证:我们可以定义更强的点查询和频繁项问题,其中误差界或频繁性判断不是基于总的 L₁ 范数 ||x||₁,而是基于“尾部”范数 ||x_tail||₁,即去掉频率最高的 k 个分量后的 L₁ 范数。Count Sketch 算法经过稍加修改的分析,同样能满足这种尾部保证,从而在数据倾斜严重时提供更精确的结果。

Lp 频繁项:我们可以将频繁性的定义从 L₁ 范数推广到 Lp 范数。研究结果表明:

  • 对于 p > 2,求解 Lp 频繁项需要多项式空间,难度很大。
  • 对于 p 越大,所能提供的误差界限越好(更紧)。因此,在内存允许的情况下,我们希望解决 p=2 的问题。
  • 存在 Count Sketch 算法(由 Charikar、Chen、Farach-Colton 提出)能以 O(k log n) 空间解决 L₂ 点查询和频繁项问题,其误差界与 L₂ 范数相关,通常比 L₁ 版本更优。

总结

本节课我们一起学习了线性素描这一设计流算法的强大框架。我们首先了解了线性素描如何与旋转门流模型结合,并通过维护一个低维素描向量来高效处理更新。接着,我们聚焦于频繁项发现问题,认识到精确求解的困难,从而转向可近似的 L₁ 频繁项和点查询问题。我们详细分析了 Count-Min Sketch 算法如何在严格旋转门模型中解决点查询,并通过取中位数将其推广到通用的旋转门模型(即 Count Sketch)。为了改进频繁项查询的效率,我们介绍了 Dyadic Trick 方法。最后,我们简要探讨了更强大的尾部保证以及将问题推广到 Lp 范数下的情况。这些线性素描技术是许多现代流处理系统的核心组件。

008:图数据流算法

在本节课中,我们将继续学习线性数据流草图技术。具体来说,我们将探讨一个图数据流算法的例子。这是一种由An、Guh和McGregor在2012年引入的技术。今天,我们将看到他们开发的一种算法,该算法使用线性数据流草图,用于在动态更新的数据流中寻找图的生成森林。自那时起,图数据流草图已被用于解决众多图问题,如匹配、生成树等。这是一个非常巧妙且出人意料的思路。

动态生成森林问题

我们考虑一个动态更新的无向图。图中的边可以被插入删除。我们的目标是,当收到查询时,能够输出当前图的一个生成森林。当然,这个生成森林可以用来回答连通性问题,例如判断两个顶点S和T之间是否存在路径,或者在无向图中找到从S到T的路径。

这个问题可以建模为严格数据流模型中的问题。我们可以定义一个向量 x,其维度为所有可能的边(即 n choose 2)。其中,x[uv] 表示边 (u, v) 在图中出现的次数(重数)。在简单图中,这个值只能是0或1。插入边 (u, v) 对应于更新操作 x[uv] += 1,删除边则对应于 x[uv] -= 1。通过本节课,你将看到如何在数据流中,使用低内存解决这个问题。

如果图中只有插入操作,这个问题是简单的。我们可以使用并查集算法。更简单的方法是:每当插入一条新边时,如果它连接了两个不同连通分量中的顶点,我们就存储这条边;否则,如果它连接了同一连通分量中的顶点,我们就可以安全地忽略它,因为在没有删除的情况下,这两个顶点将始终连通。一个生成森林最多包含 n-1 条边,因此我们最多只需要存储 O(n) 条边,每条边需要 O(log n) 位来存储端点,总内存为 O(n log n) 位。这可以被证明是插入流模型下的最优解。

然而,在严格数据流模型(允许删除)中,情况就复杂了。我们不能简单地只记住生成森林中的边。考虑一个简单的例子:先插入边1,再插入边2,然后插入边3。由于边3连接的两个顶点已经连通,我们可能决定不存储它。接着,如果我们删除了边1和边2,图就只剩下边3,但我们因为没有存储它,会误以为图是空的。因此,当删除操作发生时,我们需要能够找到“替代边”来维持连通性。

令人惊讶的是,存在一种随机化算法,即AGM草图,它仅使用 O(n log^3 n) 位内存,就能以至少 1 - 1/poly(n) 的概率成功解决动态生成森林问题。这个 n log^3 n 的内存界限被证明是最优的。有趣的是,即使只是回答图是否连通这一个比特的问题,也需要至少 Ω(n log^3 n) 位内存。

AGM草图不仅适用于流模型,也适用于分布式计算模型。在该模型中,每个顶点只知道自己的邻居,一个中央服务器需要基于各顶点发送的短消息来重建生成森林。AGM草图表明,每个顶点只需发送长度为 O(log^3 n) 位的消息即可,这远优于发送整个邻接表(需要 Ω(n) 位)的朴素方法。

所需的基本子程序

在深入AGM草图的工作原理之前,我们需要介绍两个基本的子程序。这两个子程序本身也是基于线性数据流草图构建的,并且在其他领域(如压缩感知)也有应用。

1. K-稀疏恢复

问题定义:在数据流模型中维护一个向量 x。当收到查询时,如果 x 的非零项数量(即支撑集大小)最多为 K,则必须精确地返回 x;否则,可以输出任意值。

解决方案思路:我们需要设计一个线性映射 Π(一个矩阵),使得对于所有支撑集大小不超过 2K 的向量 z,都有 Πz ≠ 0。这意味着,从 K-稀疏向量到其草图 Πx 的映射是单射。

一种构造方法是使用转置的范德蒙德矩阵。设我们有一个素数 p > 2n,并在有限域 F_p 上运算。定义矩阵 Π2K 行,第 i 列为 [1, a_i, a_i^2, ..., a_i^{2K-1}]^T,其中所有的 a_i 互不相同。对于任意 2K 列组成的子矩阵,其行列式是一个非零的乘积(因为 a_i 互异),因此该子矩阵是满秩的,没有非零向量在其零空间中。这就保证了所需的单射性质。

查询算法:给定草图 y = Πx,我们遍历所有大小为 2K 的列索引子集 S。对于每个 S,我们求解线性系统 Π_S * z = y,其中 Π_S 是只保留 S 中列的矩阵。如果解出的 z 的支撑集大小不超过 K,我们就将其补零后返回作为 x。由于映射是单射,只有一个 S 会给出正确的 K-稀疏解。

内存:存储 Πx 需要 O(K log p) = O(K log n) 位。
缺点:查询时间是指数级的 O(n^{2K})。虽然存在多项式时间的解码算法(如与里德-所罗门码解码相关的算法),但本节课我们只关注内存消耗。

2. 支撑集查找

问题定义:在数据流模型中维护向量 x。当收到查询时,返回任意一个满足 x[i] ≠ 0 的索引 i。算法应以至少 1 - δ 的概率成功。

解决方案思路:我们使用几何分桶采样的策略,这与我们在学期初估计不同元素数量时用到的技巧类似。

我们准备 log nK-稀疏恢复数据结构,记为 KSR_1, KSR_2, ..., KSR_{log n}。此外,我们使用一个哈希函数 h,它将每个可能的索引 i 映射到 {1, 2, ..., log n},并且满足 Pr[h(i) = j] ≈ 1/2^j(最后一级的概率会调整以使总和为1)。哈希函数只需要 log n 阶独立即可。

更新:对于每个更新 (i, Δ),我们将其送入所有 KSR_j 中满足 h(i) = j 的那个数据结构。
查询:我们从最底层(j = log n)开始向上检查每个 KSR_j。当我们第一次遇到一个返回非零向量的 KSR_j 时,我们就从该向量中任意选取一个非零项对应的索引 i 返回。如果所有 KSR_j 都返回零向量,则我们判断 x 本身可能就是零向量。

参数设置与正确性:我们将每个 KSR_jK 设置为 O(log(1/δ))。对于底层(j 较大)的桶,期望落入的元素数量指数级减少,切尔诺夫界保证了它们以极高概率是空的或元素极少(不超过 K)。对于中间某个“红色”层级,我们期望有大约 Θ(log(1/δ)) 个元素落入,切尔诺夫界同样保证它既非空,元素数量也不超过 K 的概率很高。通过设置合适的常数并联合所有层级的失败概率,我们可以使总失败概率低于 δ

内存:每个 KSR_j 使用 O(K log n) = O(log(1/δ) log n) 位。共有 log n 个,总内存为 O(log(1/δ) log^2 n) 位。对于支撑集查找问题,这个空间复杂度在大多数情况下是最优的。

AGM草图:动态生成森林算法

现在,我们利用上述子程序来构建解决动态生成森林的AGM草图算法。核心思想是模拟一个多轮次的离线算法,每轮中每个(超级)顶点尝试选择一条连向其他组件的边进行合并。

离线算法框架

  1. 初始时,每个顶点自成一个组件。
  2. 对于每一轮:
    • 每个组件(超级顶点)尝试选择一条离开本组件的边(如果存在)。
    • 将所有被选中的边加入生成森林。
    • 沿着这些边,将相连的组件收缩成新的超级顶点。
  3. 重复此过程,直到没有新的边可以被选中(即所有顶点都已连通或成为孤立点)。

可以证明,由于每轮至少能使非最大组件的数量减半,所以总轮数 R = O(log n)

在线流算法模拟

我们需要在数据流中模拟这个过程,且只存储线性草图。

顶点向量表示:对于每个顶点 u,我们定义一个长度为 n choose 2有符号邻接向量 x_u。对于每条可能的边 (u, v)

  • 如果 (u, v) 不是图中的边,则 x_u[(u,v)] = 0
  • 如果 (u, v) 是图中的边:
    • u < v,则 x_u[(u,v)] = -1
    • u > v,则 x_u[(u,v)] = +1
      (符号约定只要一致即可)。这个设计的关键在于,当两个顶点 uv 合并时,x_u + x_v 中连接 uv 的边会因符号相反而抵消,结果向量的支撑集恰好是离开这个新组件的所有边。

草图存储:我们并不存储完整的 x_u,而是存储 R 个独立的线性草图。对于 r = 1 到 R,我们为每个顶点 u 存储 Π_r * x_u,其中 Π_r 是一个支撑集查找数据结构的草图矩阵。因此,总存储空间为:n(顶点数) * R(轮数) * O(log(1/δ) log^2 n)(每个草图大小)。通过精细设置 δ = 1/poly(n) 并利用“已知失败/静默失败”两种模式的支撑集查找变体,可以将总空间优化至 O(n log^3 n) 位。

算法执行

  1. 第1轮:对于每个顶点 u,查询 Π_1 * x_u 对应的支撑集查找数据结构,得到一条与 u 关联的边(如果存在)。这模拟了每个顶点选择一条边。
  2. 合并与更新:对于本轮中被选中的边 (u, v),我们将顶点 uv 合并。在草图层面,这意味着对于后续所有轮次 r >= 2,我们将 Π_r * x_uΠ_r * x_v 相加,作为代表这个新组件的超级顶点的草图。实际上,我们维护每个当前组件的草图,它是组内所有顶点草图的
  3. 第 r 轮:对于当前每个组件(超级顶点)A,我们查询其对应的草图 Π_r * (sum_{u in A} x_u)。根据之前的性质,这个草图的支撑集查找结果,就是一条离开组件 A 的边(如果存在)。我们选择这些边进行新一轮的合并。
  4. 重复步骤2和3,共进行 R = O(log n) 轮。

为什么需要多个独立草图? 一个关键点是,每一轮选择的合并边依赖于该轮所用草图 Π_r 的随机性。如果我们只在所有轮次中重复使用同一个 Π,那么后续轮的查询输入(即上一轮合并后的组件)将与 Π 的随机性相关。这违反了蒙特卡洛随机算法要求“输入独立于随机硬币”的前提,导致正确性保证失效。使用独立的 Π_1, Π_2, ..., Π_R 可以确保每一轮的查询与当前轮的随机性独立,从而满足算法分析的条件。从下界角度看,这也确实是必需的,因为如果单个草图就够用,将得到 O(n log^2 n) 的算法,而这与已知的 Ω(n log^3 n) 下界矛盾。

总结

本节课我们一起学习了图数据流算法中的核心方法——AGM草图。我们首先明确了动态生成森林问题的定义及其在严格数据流模型下的挑战。接着,我们引入了两个基础构建块:K-稀疏恢复支撑集查找,并了解了如何用线性草图实现它们。最后,我们详细阐述了如何利用多个独立的支撑集查找草图,以多轮合并的方式模拟离线算法,从而在 O(n log^3 n) 位内存内解决动态生成森林问题。这个算法优美地结合了图论、概率分析和线性草图技术,是数据流算法中的一个经典成果。

009:解耦、Hanson-Wright证明、ℓp范数估计、Nisan的PRG

在本节课中,我们将要学习几个核心内容。首先,我们将证明Hanson-Wright不等式,这是我们周一用来证明分布式Johnson-Lindenstrauss引理的关键工具。完成之后,我们将探讨其他范数(ℓp范数,其中p≠2)的估计问题,特别是0 < p < 2的情况。最后,我们将讨论用于空间有界计算的伪随机生成器,这是一种比k-wise独立性更通用的去随机化工具。

Hanson-Wright不等式的证明

上一节我们概述了本节课的目标。本节中,我们来看看如何证明Hanson-Wright不等式。该不等式给出了当随机变量为±1时,二次型偏离其期望值的尾概率界。

Hanson-Wright不等式陈述如下:设σ₁, ..., σₙ为独立均匀的±1随机变量,A是一个n×n矩阵。那么对于所有λ > 0,有:

Pr[ |σᵀAσ - E[σᵀAσ] | > λ ] ≤ C * exp( -c * min( λ²/||A||_F², λ/||A|| ) )

其中,||A||_F是A的Frobenius范数,||A||是A的算子范数。这个界是高斯尾(e{-λ²})和指数尾(e)的混合。

尾界和矩界是等价的。证明一个形如上式的尾界,等价于证明该随机变量的所有p阶矩都有界。我们将通过界定所有p范数来证明Hanson-Wright不等式。

矩与尾的等价性

首先,我们简要说明尾界和矩界的等价关系。对于一个实值随机变量Z,其p范数定义为:

||Z||_p = ( E[ |Z|^p ] )^{1/p}

如果已知尾概率界 Pr[|Z| > λ] ≤ exp(-cλ²/α²),那么可以通过积分推导出矩界 ||Z||_p ≤ O(α√p)。反之,如果已知矩界 ||Z||_p ≤ O(α√p),那么通过选择适当的p并应用马尔可夫不等式,可以得到尾概率界 Pr[|Z| > λ] ≤ exp(-c‘λ²/α²)。因此,证明Hanson-Wright不等式等价于证明二次型 σᵀAσ - E[σᵀAσ] 的p范数满足 O(√p ||A||_F + p ||A||) 的界。

Khintchine不等式

在证明二次型之前,我们先回顾一个关于线性型的结果——Khintchine不等式。设σ₁, ..., σₙ为独立均匀的±1随机变量,x是一个固定向量。那么对于p ≥ 1,有:

|| ∑ σ_i x_i ||_p ≤ O(√p) * ||x||_2

这个不等式表明,Rademacher线性型的矩增长不会快于高斯线性型。其证明思路是将Rademacher随机变量的矩与高斯随机变量的矩进行比较,利用高斯矩的已知结果。

解耦技术

现在,我们转向二次型。证明的关键技术是“解耦”。解耦定理指出,对于非对角矩阵A(即对角线元素为0),有:

|| ∑_{i≠j} A_{ij} σ_i σ_j ||_p ≤ 4 * || ∑_{i≠j} A_{ij} σ_i σ_j‘ ||_p

其中,σ和σ‘是两组独立的±1随机变量。这个定理允许我们将一个二次型中的一个σ因子替换为一个独立的副本,代价只是一个常数因子。

以下是证明的概要:

  1. 引入独立的伯努利随机变量η_i ~ Bernoulli(1/2)。
  2. 将原表达式重写为 4 * E_η[ || ∑_{i≠j} A_{ij} η_i (1-η_j) σ_i σ_j ||_p ]
  3. 利用詹森不等式将期望移到范数外。
  4. 由于期望是平均值,存在一个固定的η‘使得该表达式值至少达到平均值。
  5. 这个固定的η‘定义了一个集合S = {i: η‘_i = 1},从而将和式限制在i∈S, j∉S上。
  6. 将j∉S的σ_j替换为独立的σ‘_j,分布不变。
  7. 再次应用詹森不等式,将期望(现在是对σ‘_S的期望)移到范数外,最终得到形如 || ∑_{i,j} A_{ij} σ_i σ‘_j ||_p 的表达式,即完成了证明。

Hanson-Wright的证明

现在,我们利用解耦和Khintchine不等式来证明Hanson-Wright的矩界。令 Z = σᵀAσ - E[σᵀAσ] = ∑_{i≠j} A_{ij} σ_i σ_j

  1. 应用解耦||Z||_p ≤ 4 * ||σᵀAσ‘||_p,其中σ‘独立于σ。
  2. 条件于σ,应用Khintchine:将σ视为固定,对σ‘应用Khintchine不等式:
    ||σᵀAσ‘||_p ≤ √p * || Aσ‘ ||_2 的p范数。
    更准确地说,是 (E_σ‘[ |σᵀAσ‘|^p ])^{1/p} ≤ √p * (E_σ‘[ ||Aσ‘||_2^p ])^{1/p}
  3. 处理 ||Aσ‘||_2:定义 E = (E_σ‘[ ||Aσ‘||_2^p ])^{1/p}。我们的目标是界定E。
    注意到 ||Aσ‘||_2^2 = σ‘ᵀ (AᵀA) σ‘
  4. 中心化并再次解耦:将 ||Aσ‘||_2^2 写成其期望(即 ||A||_F^2)加上一个中心化的二次型。对这个中心化的二次型再次应用解耦技术。
  5. 再次应用Khintchine:对解耦后的表达式再次应用Khintchine不等式,最终得到一个关于E的不等式:
    √p * E² ≤ O( √p * ||A||_F + p^{3/4} * ||A||^{1/2} * E )
  6. 求解二次不等式:上述不等式可以视为关于E的二次不等式 E² - αE - β ≤ 0。这意味着E必须小于该二次方程的正根,从而得到:
    E ≤ O( ||A||_F + √p * ||A|| )
  7. 代回完成证明:将E的界代回最初的表达式,最终得到:
    ||Z||_p ≤ O( √p * ||A||_F + p * ||A|| )
    根据矩与尾的等价性,这就证明了Hanson-Wright不等式。

这个证明展示了如何通过解耦技术将二次型问题约化为线性型问题,并利用已知的矩不等式进行推导。

ℓp范数估计(p ≠ 2)

上一节我们证明了Hanson-Wright不等式。本节中,我们来看看如何估计数据流中向量的ℓp范数,其中p不等于2。

p稳定分布

核心思想是使用“p稳定分布”。一个分布D被称为p稳定的(0 < p ≤ 2),如果满足以下性质:设Z₁, ..., Zₙ i.i.d. ~ D,对于任意固定向量x ∈ ℝⁿ,有:

∑_{i=1}^n x_i Z_i ~ ||x||_p * Z

其中,Z ~ D,||x||_p = (∑ |x_i|^p)^{1/p}。这意味着,用p稳定随机变量对向量x进行加权和,结果的分布等同于x的ℓp范数乘以一个p稳定随机变量。

  • 当p=2时,高斯分布是2稳定的。
  • 当p=1时,柯西分布是1稳定的。
  • 对于其他0<p<2,也存在p稳定分布,但其概率密度函数可能没有闭式解,其特征函数(傅里叶变换)具有形式 exp(-|t|^p)

Indyk的算法

Indyk在2000年提出了以下算法来估计ℓp范数:

  1. 设Π是一个m×n的矩阵,其每个元素独立采样自一个p稳定分布D_p,并且我们对该分布进行缩放,使得 Pr[|Z| ≤ 1] = 1/2(即中位数为1)。
  2. 在数据流中,维护草图 y = Π x。这是一个m维向量。
  3. 查询时,输出 median(|y₁|, |y₂|, ..., |y_m|) 作为 ||x||_p 的估计值。

算法分析

为什么这个算法有效?考虑y的第r个分量 y_r = ∑_{i=1}^n Π_{ri} x_i。根据p稳定性,y_r 的分布与 ||x||_p * Z 相同,其中Z ~ D_p。

由于我们缩放D_p使得 Pr[|Z| ≤ 1] = 1/2,因此对于每个r,有 Pr[|y_r| ≤ ||x||_p] = 1/2

现在,考虑稍微放宽的区间:

  • Pr[|y_r| ≤ (1+ε)||x||_p] = 1/2 + Θ(ε) (概率质量增加了约ε)
  • Pr[|y_r| ≤ (1-ε)||x||_p] = 1/2 - Θ(ε) (概率质量减少了约ε)

S_+ = (1/m) * #{ r : |y_r| ≤ (1+ε)||x||_p }S_- = (1/m) * #{ r : |y_r| ≤ (1-ε)||x||_p }

  • E[S_+] = 1/2 + Θ(ε)
  • E[S_-] = 1/2 - Θ(ε)

由于y的各分量独立,S_+S_-的方差均为 O(1/m)。根据切比雪夫不等式,只要 m = Ω(1/ε²)S_+S_-就会以高概率分别大于1/2和小于1/2。这意味着超过一半的 |y_r| 小于 (1+ε)||x||_p,同时少于一半的 |y_r| 小于 (1-ε)||x||_p,因此中位数必然落在 [(1-ε)||x||_p, (1+ε)||x||_p] 区间内。

去随机化

Indyk的原始算法需要完全独立的p稳定随机变量,这需要大量随机位。去随机化的思路是使用伪随机生成器。

  1. 两两独立性:算法分析中,仅在对S_+S_-进行方差分析时使用了独立性。方差可加性只需要两两独立性。因此,我们可以让矩阵Π的不同行由两两独立的种子生成。
  2. k-wise独立性:要保证p稳定性在k-wise独立下近似成立,是一个更深入的结果。Kane、Nelson和Woodruff在2010年证明:如果使用 k = Ω(1/ε^p)-wise独立的p稳定随机变量,那么其累积分布函数与完全独立的情形仅相差ε。
  3. Nisan的PRG:在Indyk的原始论文(2000年)中,他使用了Nisan的伪随机生成器,这是一种适用于空间有界计算的通用去随机化工具。

Nisan的伪随机生成器

上一节我们讨论了ℓp估计及其去随机化需求。本节中,我们来看看Nisan的伪随机生成器,它是一种用于空间有界计算的通用去随机化工具。

分支程序模型

我们使用只读一次分支程序来建模空间有界的流算法。

  • 程序有L层,对应处理L个输入块(每个块包含S比特,代表内存状态的所有可能)。
  • 每层有 2^S 个节点,代表所有可能的内存配置(因为内存大小为S比特)。
  • 从一层到下一层,根据当前读入的S比特输入,决定转移到下一个内存状态。
  • 程序从初始状态开始,经过L步后到达最终状态。

Nisan定理

Nisan在1992年证明了以下定理:
对于一个宽度为 2^S、层数为L的只读一次分支程序P,存在一个伪随机生成器G。G使用 Q = O(S * log L) 个真正的随机比特,并输出 L * S 个“看起来随机”的比特。对于该分支程序P,有:

|| Distribution(P(Uniform(L*S bits))) - Distribution(P(G(Uniform(Q bits)))) ||₁ ≤ ε

其中ε可以很小(例如 2^{-S})。也就是说,用G生成的伪随机串喂给程序P,其最终状态分布与使用完全随机串时几乎不可区分。

应用意义

这意味着,如果一个流算法使用S比特内存,处理长度为L的输入,并且其正确性分析依赖于输入的随机性,那么我们可以用Nisan的PRG来生成这些随机比特。我们只需要存储 O(S log L) 个真正的随机比特(以及PRG的小状态),而不是 L * S 个比特。当S很小而L很大时,这带来了巨大的空间节省。

在Indyk的ℓp估计算法中,可以将矩阵Π的生成过程视为一个分支程序。使用Nisan的PRG可以显著减少生成整个随机矩阵所需存储的随机种子大小,从而实现算法的去随机化。

总结

本节课中我们一起学习了以下几个核心内容:

  1. Hanson-Wright不等式的证明:我们通过解耦技术将二次型的矩界问题转化为线性型问题,并利用Khintchine不等式完成了证明。该不等式是分析Rademacher二次型尾分布的重要工具。
  2. ℓp范数估计(p≠2):我们了解了利用p稳定分布的特性来估计数据流向量的ℓp范数。Indyk的算法通过维护一个由p稳定随机变量构成的线性草图,并输出其中位数作为范数估计。
  3. 去随机化技术:我们讨论了如何减少上述算法对完全独立随机性的依赖,包括使用两两独立性、k-wise独立性,以及通用的Nisan伪随机生成器。Nisan的PRG为空间有界的计算提供了一种强大的去随机化方法,只需消耗少量真实随机性即可模拟大量随机比特的效果。

这些工具和技术在数据流算法设计和分析中非常重要,它们帮助我们在资源受限的环境下进行高效的随机化计算和去随机化。

010:嘉宾讲座 - Ilya Razenshteyn

概述

在本节课中,我们将学习Ilya Razenshteyn关于范数空间嵌入与草图之间关系的研究。核心定理表明,如果一个有限维范数空间存在一个高效的通信协议(草图)用于距离估计,那么该空间可以线性嵌入到一个更简单的空间(如L1-ε空间)中,且失真度可控。这为证明通信复杂度下界提供了一种新方法。


通信模型与草图简介

上一节我们概述了本节课的核心内容。本节中,我们来看看本讲座所基于的通信模型。

我们有一个度量空间 (X, dist)。通信问题涉及Alice和Bob,他们各自拥有该空间中的一个点。他们希望通过交换尽可能少的比特信息,来判断两点间的距离是“近”(≤ 1)还是“远”(≥ D)。

该模型有几个关键点:

  • 协议是随机的,成功概率至少为2/3。
  • 允许使用公共随机硬币,且数量是可数无穷的。这使得协议能够生成如高斯分布这样的连续随机变量样本。

我们将协议交换的比特数称为草图大小通信复杂度 S,将区分阈值 D 称为近似比


已知的草图构造

在深入核心定理前,我们先回顾一些已知的高效草图构造。以下是两个经典例子:

欧几里得空间 (L2) 草图

对于d维欧几里得空间,存在一个高效的草图协议。

  1. 生成随机向量:Alice和Bob利用公共随机性,独立采样 g 个标准高斯随机向量 r_i ∈ R^d
  2. 计算投影:Alice计算 x̃_i = <x, r_i>,Bob计算 ỹ_i = <y, r_i>
  3. 离散化与通信:双方将实数轴随机划分为大小为 w 的区间,并为每个区间随机分配 +1-1 标签。Alice发送 x̃_i 所在区间的标签,Bob发送 ỹ_i 所在区间的标签。
  4. 判断:如果所有 g 次比较中,双方标签均相同,则判定为“近”;否则判定为“远”。

通过设置 g = O(1/ε²) 和合适的 w,该协议能以高概率实现 (1+ε) 的近似比,草图大小为 O(1/ε²)

Lp 空间草图 (p ≤ 2)

对于 p ≤ 2 的Lp范数,构造是类似的。只需将高斯随机向量替换为 p-稳定分布 的随机向量,其余步骤和分析与L2情况完全相同。最终也能获得 (1+ε) 近似比和 O(1/ε²) 的草图大小。


核心定理陈述

现在,我们来看Ilya Razenshteyn等人证明的核心定理。

定理:设 (X, ||·||) 是一个有限维范数空间。假设存在一个公共随机硬币通信协议,能够以草图大小 S 和近似比 D 解决 X 上的距离估计问题。那么,对于任意 ε > 0,存在一个从 XL_{1-ε} 空间的线性嵌入 F,其失真度至多为 O(S * D / ε)

核心概念解释

  • 线性嵌入:映射 F: X -> L_{1-ε} 是一个线性映射,即满足 F(x+y) = F(x) + F(y)
  • 失真度:存在常数 c,使得对于所有 x, y ∈ X,有 (1/c) * ||x-y|| ≤ ||F(x)-F(y)||_{L_{1-ε}} ≤ c * ||x-y||。定理中的 O(S*D/ε) 就是这个 c 的上界。
  • L_{1-ε} 空间:这不是一个范数空间(不满足三角不等式),但对于它,我们仍然有高效的草图(使用p-稳定分布,p=1-ε)。

定理的意义:该定理建立了草图存在性可嵌入到简单空间之间的等价关系。它表明,对于范数空间,获得高效距离估计草图的唯一通用方法,本质上就是先将其嵌入到 L1L_{1-ε},然后使用基于稳定分布的草图。这为证明通信复杂度下界提供了强大工具:要证明某个范数空间没有小尺寸草图,只需证明它不能低失真地嵌入到 L_{1-ε} 即可。


定理证明概览

证明分为两大步骤,我们重点概述第二步的几何部分。

第一步:从协议到“弱嵌入”

首先,假设存在大小为 S、近似比为 D 的协议。目标是构造一个到希尔伯特空间(近似为 L2)的映射 f,它满足:

  • ||x-y|| ≤ 1,则 ||f(x)-f(y)|| ≤ 1
  • ||x-y|| ≥ C*S*D(C为大常数),则 ||f(x)-f(y)|| ≥ 10

这个映射 f 只在一个距离尺度上保存信息,并且不一定是线性的。这一步的证明使用了信息论和通信复杂度的工具,特别是通过构造“困难分布对”并利用反证法来完成。由于涉及较多专业背景,此处不展开。

第二步:从“弱嵌入”到“强线性嵌入”

第二步是证明的几何核心。它接收第一步产生的“弱嵌入” f,并将其转化为最终所需的、适用于所有距离尺度的线性嵌入 F: X -> L_{1-ε}

这一步又分为两个子步骤:

子步骤2.1:正则化映射
首先,将 f 转化为一个性质更好的映射 f̃: X -> HH 为某希尔伯特空间)。 满足:

  1. 上半赫尔德条件||f̃(x)-f̃(y)|| ≤ O( sqrt(||x-y||) )
  2. 下半界:当 ||x-y|| 很大时,||f̃(x)-f̃(y)|| 也有一个常数下界。

这个构造使用了两个几何工具:

  • 事实一:欧几里得度量的平方根仍然是欧几里得度量(可等距嵌入到希尔伯特空间)。这是Schoenberg定理的结果。
  • 事实二:定义在度量空间子集上的 1/2-赫尔德映射,可以延拓到整个空间上,且保持赫尔德常数。这是Kirszbraun定理在赫尔德映射上的推广。

子步骤2.2:利用调和分析完成线性嵌入
现在,我们有一个满足上述条件的 。目标是构造线性嵌入 F

  1. 转化为核函数:利用另一个Schoenberg定理,可以构造一个从希尔伯特空间到其单位球面的映射 g,使得 <g(u), g(v)> = exp(-||u-v||²)。将 g 复合,得到映射 h(x),其点积性质为:<h(x), h(y)> ≈ exp(-||x-y||)
  2. 平移不变性与Bochner定理:通过对称化,可以将上述点积核函数转化为一个平移不变的正定函数 φ(x-y) = <h(x), h(y)>。根据Bochner定理,这样的函数必然是某个概率分布 μ 的傅里叶变换的特征函数。即存在分布 μ,使得 φ(z) = E_{v~μ}[exp(i <v, z>)]
  3. 分析分布 μ 的性质:可以证明,从这个分布 μ 中采样出的随机向量 v,与固定向量 x 的点积 <v, x> 具有良好的性质:
    • 以常数概率,|<v, x>| 不小于 Ω(||x|| / D)
    • 具有柯西型的尾部:P(|<v, x>| > t * ||x||) ≤ O(1/t)
  4. 构造最终嵌入:定义线性映射 F: X -> L_{1-ε}(μ)[F(x)](v) = <v, x>。这里,L_{1-ε}(μ) 范数定义为 ||F(x)||_{1-ε} = (∫ |<v, x>|^{1-ε} dμ(v))^{1/(1-ε)}。利用上一步得到的点积性质,可以验证 F 确实是一个失真度为 O(S*D/ε) 的嵌入。其中,1-ε 次幂的引入正是为了处理柯西尾部,确保积分收敛且范数受控。

总结

本节课我们一起学习了Ilya Razenshteyn关于范数空间草图与嵌入的深刻定理。

  • 我们首先回顾了基于随机投影和稳定分布的距离估计草图。
  • 然后,我们学习了核心定理:有限维范数空间存在高效草图,当且仅当它可以低失真地线性嵌入到 L_{1-ε} 空间。失真度上界为 O(草图大小 * 近似比 / ε)
  • 最后,我们概览了定理的证明思路,其核心是将通信协议的信息论下界,转化为几何上的嵌入障碍,并利用Schoenberg定理、延拓定理和Bochner定理等经典分析工具,构造出所需的线性嵌入。

这一定理不仅在哲学上揭示了距离估计草图的内在结构,也为证明通信复杂度下界提供了强大而优雅的新工具。

011:Johnson-Lindenstrauss下界 🎯

在本节课中,我们将学习Johnson-Lindenstrauss引理的下界。我们将回顾其历史发展,理解证明下界的基本思想,并最终探讨如何证明一个与JL上界匹配的最优下界。

概述

Johnson-Lindenstrauss引理表明,任何n个点的集合都可以嵌入到O(ε⁻² log n)维的欧几里得空间中,同时保持任意两点间的距离在(1±ε)因子内。一个自然的问题是:这个上界是否是最优的?是否存在一个点集,使得任何低失真嵌入都必须使用至少Ω(ε⁻² log n)维的目标空间?本节课我们将证明答案是肯定的,JL引理确实是最优的。

历史回顾

首先,我们回顾一下关于JL下界的研究历史:

  • 1984年 (JL原始论文):证明了当ε小于某个常数(例如1/3)时,存在点集要求目标维度至少为Ω(log n)。这个论证对ε不敏感。
  • 2003年 (Noga Alon):证明了存在点集要求目标维度至少为Ω(ε⁻² log n / log(1/ε))。这几乎是最优的,但与JL上界相比,分母上多了一个log(1/ε)因子。
  • 2016年 (Kasper Green Larsen, Jelani Nelson等人):证明了对于线性映射,存在点集要求目标维度至少为Ω(ε⁻² log n)。这个下界仅针对线性嵌入。
  • 2017年 (Kasper Green Larsen, Jelani Nelson等人):证明了对于任意映射(包括非线性映射),存在点集要求目标维度至少为Ω(ε⁻² log n)。这是最终的最优下界。

我们将要学习的证明基于压缩论证,其核心思想是:假设一个大的点集族中的每一个点集都存在一个“优于JL”的低维嵌入,那么我们可以利用这些嵌入构造一个从大集合到小集合的单射,从而产生矛盾。

预备知识

在进入核心证明之前,我们需要了解一些凸几何和度量空间的基本概念。

凸体与范数

一个凸体是ℝᵈ中的一个紧致、凸且具有非空内部的子集。紧致在ℝᵈ中等价于闭且有界。

一个凸体K是对称的,如果x ∈ K意味着 -x ∈ K。对称凸体必然包含原点。

关键事实:对称凸体与范数一一对应。

  • 给定一个范数||·||,其单位球{ x : ||x|| ≤ 1 }是一个对称凸体。
  • 给定一个对称凸体K,可以定义一个范数:对于x ∈ ℝᵈ,||x||_K = inf { t > 0 : x/t ∈ K }。

覆盖数与熵数

设(X, d)是一个度量空间,T ⊆ X。T的一个ε-网是一个子集T‘ ⊆ T,使得对于T中的每一个点x,在T’中都存在一个点x‘满足d(x, x’) ≤ ε。

T的熵数 N(T, d, ε) 是T的最小ε-网的大小。

填充数与体积论证

T的填充数 P(T, d, ε) 是可以在T中放置的、半径为ε、两两不相交的球的最大数量(球心必须在T内)。

关键关系:对于任何度量空间,有 N(T, d, ε) ≤ P(T, d, ε/2)。证明思路是:取一个最大的半径为ε/2的填充,然后将这些球的中心作为ε-网的候选点。如果有点未被覆盖,则可以添加一个新的半径为ε/2的球,与最大性矛盾。

体积论证:如果我们在某个范数(例如ℓ₂范数)的单位球中进行填充,那么所有半径为ε/2的填充球都包含在一个半径为(1 + ε/2)的大球内。通过比较体积,我们可以得到填充数(从而熵数)的上界:
N(单位球, ||·||, ε) ≤ P(单位球, ||·||, ε/2) ≤ ( (1 + ε/2) / (ε/2) )^d = (1 + 2/ε)^d

JL原始下界证明(1984)

原始JL论文中的下界证明是一个经典的体积论证。他们考虑的点集X由零向量和标准基向量组成:{0, e₁, e₂, ..., e_n}。假设存在一个嵌入f: X → ℝᵐ,失真为(1±ε),且ε < 1/3。通过平移,可以假设f(0) = 0。

令ẽ_i = f(e_i)。根据失真条件,我们有:

  • ||ẽ_i|| ≈ 1 ± ε (因为d(e_i, 0)=1)
  • ||ẽ_i - ẽ_j|| ≥ (1 - ε)√2 (因为d(e_i, e_j)=√2)

现在考虑围绕每个ẽ_i、半径为 r = (1 - ε)√2 / 2 的球。由于任意两个ẽ_i之间的距离至少为2r,这些球是两两不相交的。同时,所有这些球都包含在一个以原点为中心、半径为 R‘ = 1 + ε + r 的大球内。

设V(r)和V(R‘)分别表示半径为r和R‘的m维球的体积。由于n个小球体积之和不能超过大球的体积,我们有:
n * V(r) ≤ V(R‘)
=> n ≤ V(R‘) / V(r) = (R‘ / r)^m
=> m ≥ log n / log(R‘/r)

由于当ε为常数时,R‘/r也是一个常数,因此我们得到 m = Ω(log n)。这个论证的局限性在于它对ε不敏感,无法得到ε⁻²的依赖关系。

最优下界证明思路(2017)

现在,我们转向证明最优下界Ω(ε⁻² log n)的核心思路。证明采用压缩论证法。

构造困难点集族

首先,我们构造一个庞大的点集族A。令d为原始维度,k = Θ(1/ε²)。对于任意一个大小为k的子集S ⊆ [d],定义向量:
y_S = (1/√k) * 1_S (即,在S中下标处为1/√k,其余为0)

我们的点序列(允许重复)形式如下:
X = (0, e₁, e₂, ..., e_d, y_{S₁}, y_{S₂}, ..., y_{S_t})
其中 t = n - d - 1,n是总点数。

点集族A由所有这样的序列X构成,其中S₁, ..., S_t是[d]中所有可能的大小为k的子集。因此,族A的大小为:
|A| = (C(d, k))^t ≈ exp( t * k * log(d/k) )

我们假设d ≈ n / log(1/ε),并且ε足够大(大于1/√n和1/√d),这使得 log(d/k) = Θ(log n)。因此,log|A| = Θ( t * k * log n ) = Θ( n/ε² * log n )。

关键观察:低失真嵌入近似保持内积

假设点集X包含0,且其他点具有单位ℓ₂范数。如果嵌入f: X → ℝᵐ具有(1±ε)失真,并且我们平移使得f(0)=0,那么f能近似保持任意两点间的内积。

证明:对于单位向量x, y ∈ X,有:
||x - y||² = ||x||² + ||y||² - 2⟨x, y⟩ = 2 - 2⟨x, y⟩
同时,由于失真,
||f(x) - f(y)||² = (1 ± ε) ||x - y||² = (1 ± ε)(2 - 2⟨x, y⟩)
另外,直接计算:
||f(x) - f(y)||² = ||f(x)||² + ||f(y)||² - 2⟨f(x), f(y)⟩ ≈ 2 - 2⟨f(x), f(y)⟩ (因为||f(x)|| ≈ 1 ± ε)
比较两式,可得 ⟨f(x), f(y)⟩ ≈ ⟨x, y⟩ ± O(ε)。

具体到我们的构造,对于e_i和y_S:

  • 如果 i ∈ S,则 ⟨e_i, y_S⟩ = 1/√k = Θ(ε)
  • 如果 i ∉ S,则 ⟨e_i, y_S⟩ = 0
    因此,在嵌入后,⟨f(e_i), f(y_S)⟩ 在两种情况下会有明显的差距(例如≥ 8ε 与 ≤ 2ε),这使我们能够从嵌入向量中解码出集合S的信息。

压缩论证

矛盾假设:假设对于族A中的每一个点集X,都存在一个嵌入 f_X: X → ℝᵐ,失真为(1±ε),且目标维度m远小于c ε⁻² log n(c为某个小常数)。

目标:利用这些嵌入{f_X},构造一个从大族A到一个小集合的单射,从而产生矛盾。

初步尝试(有问题):对于每个X ∈ A,其编码就是所有嵌入向量的列表:Enc(X) = (f_X(0), f_X(e₁), ..., f_X(e_d), f_X(y_{S₁}), ..., f_X(y_{S_t}))。这是一个长度为n的序列,每个元素是ℝᵐ中的一个向量。根据关键观察,我们可以从这些向量的内积中恢复出所有的S_j,从而唯一确定X。然而,这个编码的问题是:每个向量f_X(·)有m个实数坐标,需要无限精度存储,编码长度并非有限比特。

第一次改进(量化):将每个嵌入向量的每个坐标四舍五入到最接近的γ的整数倍,其中γ是一个非常小的数(例如 γ = ε/√n)。可以证明,只要γ足够小,量化后的向量˜f_X(·)之间的内积仍然能保持足够差距,以区分 i ∈ S 和 i ∉ S。现在,每个坐标只有O(1/γ)种可能取值,因此整个编码可以用大约 n * m * log(1/γ) 比特来表示。这给出了一个有限比特的编码。

通过信息论论证,编码长度必须至少为log|A| = Θ( n/ε² * log n )。因此:
n * m * log(1/γ) ≥ Θ( n/ε² * log n )
=> m ≥ Θ( log n / (ε² log(1/γ) ) )
代入 γ = ε/√n,得到 m ≥ Θ( ε⁻² log n / (log(1/ε) + log log n) )。这接近但略弱于Noga Alon的下界。

第二次改进(使用ℓ₂网):第一次改进本质上是使用了ℓ∞范数下的γ-网。更好的方法是使用ℓ₂范数下的ε-网。我们知道所有嵌入向量f_X(·)的范数大约为1,因此它们都位于某个ℓ₂球内。我们取该球的一个最小ε-网。对于每个嵌入向量,我们存储它在网中最近点的索引。

存储一个索引需要的比特数是网的大小的对数,即熵数的对数。根据之前的体积论证,对于m维ℓ₂球,熵数 N ≤ (1 + 2/ε)^m。因此,存储所有n个向量的网索引所需的总比特数约为 n * log N ≤ n * m * log(1/ε)

再次应用信息论下界:
n * m * log(1/ε) ≥ Θ( n/ε² * log n )
=> m ≥ Θ( ε⁻² * (log n / log(1/ε)) )
这正是Noga Alon在2003年证明的下界。

最终改进(达到最优):为了去除分母上的log(1/ε)因子,达到Ω(ε⁻² log n)的最优下界,需要更精巧的构造。核心思想是使用一个更高效的“编码簿”,它不是简单地存储整个ε-网,而是利用嵌入向量之间的相关性进行联合编码。具体细节涉及更复杂的概率方法,但其基本精神仍然源于压缩论证。

总结

本节课我们一起学习了Johnson-Lindenstrauss引理的下界。我们从历史发展入手,理解了证明下界的核心工具——体积论证和压缩论证。我们详细剖析了JL原始下界证明,并学习了如何通过构造一个庞大的困难点集族,并假设它们都存在“好”的嵌入,进而利用压缩论证推导出维度下界。我们看到了如何从简单的量化想法出发,逐步优化到使用覆盖网,最终触及最优下界的证明思路。这个证明过程展示了理论计算机科学中一种强大而优美的证明范式。

012:Johnson-Lindenstrauss下界总结与稀疏Johnson-Lindenstrauss变换 🎯

在本节课中,我们将完成Johnson-Lindenstrauss(JL)最优下界的证明,并探讨如何加速实际的JL嵌入过程,特别是介绍稀疏Johnson-Lindenstrauss变换(SJLT)。

Johnson-Lindenstrauss下界证明总结 📊

上一节我们介绍了JL下界证明的核心思路。本节中,我们来详细完成这个证明。

证明的组织方式不是针对单个点集,而是针对一个点集集合。我们将证明,如果该集合中的每一个点集都能以低失真嵌入到低维空间,那么就会存在一个从该点集集合到一个小集合的单射,从而产生矛盾。

点集集合的构造

我们构造的集合包含有序的点序列。每个序列由以下向量组成:

  • 零向量 0
  • d 个标准基向量 e_i(第 i 个位置为1,其余为0)。
  • n - d - 1 个向量 y_S1, y_S2, ..., y_S(n-d-1)

其中,y_S 是集合 S 的指示向量,经过缩放使其具有单位L2范数。集合 S 的大小 k = 1/(100ε²)。指示向量 y_S 在位置 i(若 i ∈ S)的值为 1/√k(即 10ε),否则为0。

可能的点序列总数是 (C(d, k))^(n-d-1),其对数约为 n log n(假设 d 约为 n / log(1/ε),稍后我们将说明如何移除这个假设)。

覆盖数与熵数

我们回顾了覆盖数(或称熵数)的概念。对于任何范数的单位球,在其自身范数下构造一个 ε-网,所需网点的数量最多约为 (1/ε)^d,其中 d 是维度。

证明的核心:构造单射

我们假设集合 X 中的每个点集 X 都存在一个低失真嵌入 f_X(映射到 m 维),使得应用 f 后,X 内所有点对的内积在加性误差 ε 内被保持。

如果 f 映射到维度 m,我们的目标是证明:这种 f 对集合中每个点集都存在的事实,意味着存在一个从该点集集合到位串(长度为 O(nm))的单射。这蕴含着 O(nm) = Ω(nk log n),从而推出 m = Ω(k log n)。代入 k ≈ 1/ε²,即得到下界 m = Ω((log n)/ε²)

现在的问题是:这个单射是什么?

方法一(离散化坐标)
嵌入后的向量 f(0), f(e1), ..., f(ed), f(y_S1), ...m 维向量,其各分量值在 [-1, 1] 区间内。我们将区间 [-1, 1] 离散化为步长为 γ(例如0.01)的网格点,并将每个向量分量四舍五入到最近的网格点。

关键在于,即使经过四舍五入,我们仍然能从这些近似向量中恢复出原始点集的信息。因为原始内积 e_i · y_Sj 要么约为 10ε(若 i ∈ S_j),要么为0。经过嵌入和四舍五入后,近似内积要么至少为 ,要么至多为 ,仍然存在可区分的间隙。通过选择 γ ≈ ε/√m,可以控制舍入误差。

然而,这种方法最终会推导出 m = Ω((log n)/(ε² log(1/ε))),分母中多了一个 log(1/ε) 因子,并非最优。

方法二(使用L2网)
我们不进行逐坐标的舍入,而是将每个嵌入后的向量 f(x) 舍入到其在一个L2范数 ε-网中的最近点 f̃(x)。可以证明,即使使用舍入后的向量计算内积,其误差仍然在 O(ε) 内,因此仍然可以区分原始的内积模式。

一个L2球的 ε-网大小约为 (1/ε)^m。因此,指定一个网点需要 O(m log(1/ε)) 比特。对于 n 个向量,总比特长度为 O(nm log(1/ε))。这必须至少为 Ω(nk log n),从而推出 m = Ω((log n)/(ε² log(1/ε)))。这仍然是次优的。

方法三(利用凸几何达到最优下界)
为了消除分母上的 log(1/ε) 因子,我们需要更精细的论证。关键思想是利用对称凸体。

  1. 定义子空间和凸体
    令矩阵 A 的行由 f̃(e1)^T, ..., f̃(ed)^T 构成(d 行,m 列)。定义子空间 EA 的列空间(E ⊆ R^d,维数至多为 m)。然后定义对称凸体 K = E ∩ B_∞^d(12ε),即子空间 E 与边长为 24εL∞ 立方体的交集。K 是一个维数至多为 m 的对称凸体。

  2. 关键观察
    对于每个 j,向量 v_j = A * f̃(y_Sj) 位于 K 中。我们不需要精确知道 v_j,只需要知道一个近似值 ṽ_j,使得 ||v_j - ṽ_j||_∞ ≤ ε,就能判断 i 是否属于 S_j

  3. 在凸体范数下取网
    我们不在整个 R^dL∞ 范数下取网,而是在凸体 K 自身的范数下取一个 (1/12)-网 K'。根据凸几何知识,覆盖一个凸体在其自身范数下所需的网点数约为 (1/ε)^(dim(K)) ≈ (12)^m

  4. 比特复杂度
    因此,指定 K' 中的一个网点只需要 O(m) 比特(因为 log(12^m) = O(m))。对于 nv_j,总比特复杂度为 O(nm)。为了解码,编码器还需要额外写出矩阵 A(用于定义子空间 E 和凸体 K)。Ad 行,每行是 m 维L2网中的一个点,需要 O(m log(1/ε)) 比特。通过设定 d = n / log(1/ε),这部分的总比特数也是 O(nm)

  5. 得到最终下界
    因此,整个编码的总比特长度为 O(nm)。这必须至少为 Ω(nk log n)。代入 k ≈ 1/ε² 并消去 n,我们得到最终的最优下界:
    m = Ω((log n)/ε²)

关于维度 d 的假设
我们之前假设了 d ≈ n / log(1/ε)。如果 d 更大,我们可以对困难的点集补零,不影响其难度。如果 d 更小(小于 Ω(1/ε²)),那么恒等映射本身就是无失真的JL嵌入,下界自然不成立。因此,我们的下界在参数范围内是最优的。

加速Johnson-Lindenstrauss嵌入 ⚡

原始的JL变换(1984年)使用随机投影,即投影到一个随机选取的 m 维子空间。实现方式通常是生成一个稠密的 m × d 随机矩阵(例如通过Gram-Schmidt正交化一组随机高斯向量)。矩阵向量乘法的时间复杂度为 O(md)。如果 d 很大,这个预处理步骤本身可能就很耗时。

因此,一个自然的问题是:能否在保持JL保证的前提下,使用更快的嵌入方法?

发展历程

  1. Achlioptas (2001): 提出了第一个稀疏化方法。他构造的矩阵中,每个条目以 1/3 的概率为 ±1,以 2/3 的概率为0。这样,矩阵在期望下有 2/3 的零元,同时保持了与稠密随机投影相同的目标维度 m(常数因子相同),将计算时间减少了约 1/3

  2. Ailon & Chazelle (2006) - 快速JL变换 (FJLT): 首次给出了渐近意义上的改进。FJLT将映射 z → Πz 的时间复杂度降至 O(d log d + m³)。其思想是,对于一般向量 z,先用一个随机矩阵(如傅里叶变换矩阵)将其“搅散”,使其质量均匀分布,然后再应用一个稀疏的采样矩阵。然而,如果原始向量 z 本身是稀疏的,FJLT的预处理步骤反而可能增加计算量。

  3. 稀疏Johnson-Lindenstrauss变换 (SJLT): 为了利用输入向量的稀疏性,Dasgupta, Kumar & Charikar (2010) 以及后续Kane & Nelson的工作提出了SJLT。他们证明,可以使用CountSketch矩阵作为JL变换矩阵。

稀疏JL变换 (SJLT) 与分析

SJLT构造
变换矩阵 Π 是一个 m × d 的矩阵,其每一列有恰好 s 个非零元(s 是一个参数)。具体来说,可以将行分成 s 个块,在每个块中,每一列有且仅有一个随机位置的 ±1 值,其余为0。这正是CountSketch矩阵的结构。

性能
可以证明,当 s = O(1/ε) 且目标维度 m = O((log(1/δ))/ε²) 时,SJLT能以至少 1-δ 的概率满足JL引理要求。计算 Πz 的时间复杂度为 O(nnz(z) * s) = O(nnz(z)/ε),其中 nnz(z)z 的非零元个数。这比稠密矩阵的 O(md) 快得多,尤其适用于稀疏向量。

下界
存在点集使得任何能成功嵌入它们的稀疏矩阵,其每列的非零元个数必须至少为 Ω((log n)/(ε log(1/ε)))。这表明SJLT在稀疏性上几乎是最优的(仅差一个 log(1/ε) 因子)。

SJLT的分析概要(基于矩方法)

我们想证明,对于单位向量 z||z||_2 = 1),随机变量 E = ||Πz||_2² - 1 以高概率很小。

  1. 误差表达式
    将矩阵 Π 的条目写为 Π_{r,i} = (1/√s) * η_{r,i} * σ_{r,i},其中 σ_{r,i} 是随机 ±1 符号,η_{r,i} 是指示该位置是否为非零的 0/1 随机变量(满足每列恰有 s 个1)。通过展开 ||Πz||_2² 并减去1,可以发现误差 E 可以写为二次型:
    E = (1/s) * Σ_{r} Σ_{i≠j} η_{r,i} η_{r,j} σ_{r,i} σ_{r,j} z_i z_j
    这可以进一步写成 E = σ^T A_η σ,其中 σ 是所有 σ_{r,i} 拼接成的向量,A_η 是一个依赖于 ηz 的块对角矩阵。

  2. 应用Hanson-Wright不等式
    我们对随机变量 σ 应用Hanson-Wright不等式。该不等式指出,对于 σ 这样的Rademacher向量和矩阵 A,其二次型的 L_p 范数有上界:
    ||σ^T A σ||_p ≤ C ( √p * ||A||_F + p * ||A||_{op} )
    其中 ||A||_F 是Frobenius范数,||A||_{op} 是算子范数。

  3. 控制范数

    • 算子范数: 由于 A_η 是块对角矩阵,其算子范数是各块算子范数的最大值。可以证明每个块的算子范数 ≤ 2/s,因此 ||A_η||_{op} ≤ 2/s
    • Frobenius范数||A_η||_F² = (1/s²) Σ_{i≠j} z_i² z_j² Q_{ij},其中 Q_{ij} = Σ_r η_{r,i} η_{r,j} 是第 i 列和第 j 列的非零元共同出现的行数。Q_{ij} 近似服从参数为 (M, s²/M²) 的二项分布(虽然不是完全独立)。通过计算其矩,可以 bound 住 ||A_η||_F 的期望。
  4. 得到尾界
    将算子范数和Frobenius范数的界代入Hanson-Wright不等式,得到 EL_p 范数上界。最后,通过马尔可夫不等式(Pr[|E| > ε] ≤ (||E||_p / ε)^p),并选择合适的 p,即可证明 Pr[|E| > ε] ≤ δ,从而完成SJLT的分析。

总结 🎓

本节课中,我们一起学习了:

  1. JL下界的最优证明:通过构造一个点集家族,并利用覆盖数/熵数以及凸几何中对称凸体的性质,我们证明了任何JL嵌入的目标维度必须至少为 Ω((log n)/ε²),这个下界是紧的。
  2. 稀疏Johnson-Lindenstrauss变换 (SJLT):为了加速嵌入计算,特别是对于稀疏数据,我们可以使用像CountSketch这样的稀疏随机矩阵。SJLT在几乎保持最优目标维度的同时,将嵌入时间复杂度从 O(md) 降到了 O(nnz(z)/ε),并且被证明在稀疏性上几乎是最优的。其分析核心是使用矩方法和Hanson-Wright不等式来控制误差。

013:快速Johnson-Lindenstrauss变换与Krahmer-Ward定理 🚀

在本节课中,我们将要学习如何通过一种称为“快速Johnson-Lindenstrauss变换”的方法来加速高维数据的降维过程。我们将探讨其核心思想、原始分析,以及后续由Krahmer和Ward提出的一个更深入的理解框架。

上一节我们介绍了通过使嵌入矩阵稀疏化来加速JL变换的方法。本节中,我们来看看另一种更早被提出的加速方法——快速Johnson-Lindenstrauss变换。

FJLT的核心思想

FJLT的基本思想是利用采样来加速。具体来说,我们考虑一个采样矩阵 S,它有 D 列(原始向量维度)和 M 行(目标维度)。S 的每一行在随机位置有一个1,其余为0。那么,S 乘以向量 z 的结果,就是从 z 中独立均匀采样 M 个分量组成的新向量。

然而,直接采样存在一个问题:如果原始向量 z 的“质量”分布非常不均匀(例如,绝大部分质量集中在一个坐标上),那么我们需要采样非常多次(M 需要接近 D)才能准确估计其范数,方差会很大。

Ailon和Chazelle的想法是,在采样之前,先对向量进行一个“预处理”。他们引入一个随机化的预条件矩阵 A,使得 A 几乎是一个正交矩阵,并且能以高概率将任何向量 z 的“质量”均匀地散布到所有坐标上。这样,后续的采样步骤就能用较小的 M 获得良好的效果。

这个预条件矩阵的构造受到了量子力学中“不确定性原理”的启发:一个向量和它的傅里叶变换不能同时在少数坐标上高度集中。因此,他们选择使用离散傅里叶变换或哈达玛变换这样的有界正交系统作为 A 的基础,并辅以随机符号翻转。

FJLT的构造

完整的FJLT嵌入矩阵 Π 构造如下:

Π = (1/√M) * S * H * D

其中:

  • S 是之前描述的采样矩阵。
  • H 是一个非归一化的有界正交系统矩阵,例如非归一化的离散傅里叶变换矩阵或哈达玛变换矩阵。关键点在于,乘以 H 可以通过快速傅里叶变换或快速哈达玛变换在 O(D log D) 时间内完成。
  • D 是一个随机对角矩阵,对角线元素是独立的 Rademacher 随机变量(以1/2概率取+1或-1)。这相当于对向量的每个分量进行随机符号翻转。

Π 作用于向量 z 的过程是:先随机翻转符号 (D z),然后通过快速变换 (H) 将能量扩散,最后进行采样 (S)。

Ailon-Chazelle 分析概要

他们的分析分为两步:

  1. 证明预处理后的向量是“良扩散”的:他们证明,以高概率(超过 1 - δ/2),经过 H D z 变换后的向量 y,其最大分量(无穷范数)被很好地控制:
    ||y||_∞ ≤ √(2 log(4D/δ))
    这意味着向量的能量没有过度集中在某个坐标上。

  2. 证明在“良扩散”条件下,采样是有效的:在第一步成立的事件下,他们使用伯恩斯坦不等式证明,采样估计的范数平方 ||Π z||² 会以高概率(超过 1 - δ/2)集中在真实范数平方 ||z||² 附近。
    综合两步,整体失败概率不超过 δ。

通过此分析,他们得到的目标维度 M 需满足:
M = O( (1/ε²) * log(1/δ) * log(D/δ) )
这与最优的JL维度 O((1/ε²) * log(1/δ)) 相比,多了一个 log(D/δ) 因子。

为了消除这个多余因子,可以采用一个技巧:先使用上述FJLT将维度从 D 降至 M' = O((1/ε²) * log(1/δ) * log(D/δ)),然后再用一个传统的、较慢的稠密随机投影(如高斯矩阵)从 M' 降至最终的 M = O((1/ε²) * log(1/δ))。由于第二阶段输入的维度 M' 已经很小,所以总的计算成本仍然主要由第一阶段的 O(D log D) 主导。

Krahmer-Ward 定理:一个不同的视角 🔍

多年后,Krahmer和Ward提出了一种完全不同的思路来分析FJLT这类构造为何有效。他们的思想来源于压缩感知领域。

他们重新组织了FJLT的括号:将 Π 视为 (1/√M) S H 这个整体矩阵先乘以随机符号对角矩阵 D。即,他们关注的是 (1/√M) S H 这个矩阵的性质。

这个矩阵的关键性质是满足限制等距性

定义(限制等距性,RIP)
一个矩阵 Φ 满足 (ε, K)-RIP,如果对于所有 K-稀疏 的向量 z(即非零元素不超过K个),都有:
(1 - ε) ||z||² ≤ ||Φ z||² ≤ (1 + ε) ||z||²

Krahmer和Ward证明了,像 (1/√M) S H 这样的矩阵以高概率满足 (ε, K)-RIP,其中 M 大约为 O(K * log D * polylog(K))

他们的核心定理如下:

Krahmer-Ward 定理
假设矩阵 Φ 满足 (ε, 2K)-RIP,其中 K = C * log(1/δ)C 为某个常数。
那么,对于任意单位范数向量 z,有:
P_σ( | ||Φ D_σ z||² - 1 | > ε ) ≤ δ
其中 D_σ 是以独立Rademacher变量为对角元的对角矩阵,概率取自随机符号向量 σ

这个定理的惊人之处在于,它将构造JL映射的问题归结为构造一个满足RIP的矩阵的问题。而我们知道,许多来自压缩感知的快速构造(如部分哈达玛矩阵)都满足RIP。

定理证明思路

证明的核心是将范数平方 ||Φ D_σ z||² 写成一个二次型:
||Φ D_σ z||² = σᵀ (D_z Φᵀ Φ D_z) σ = σᵀ X σ
其中 D_z 是以 z 的分量为对角元的对角矩阵,X = D_z Φᵀ Φ D_z

接下来,将向量 z 的分量按绝对值从大到小排列,并分成大小为 K 的块。相应地,将矩阵 XK×K 的块进行分块。然后将 X 分解为四个部分之和:
X = A + B + Bᵀ + C
其中 A 包含所有对角块,B 包含紧邻对角线上方的那些块,C 包含剩余的所有远离对角线的块。

于是,二次型分解为:
σᵀ X σ = σᵀ A σ + σᵀ B σ + σᵀ Bᵀ σ + σᵀ C σ

证明的目标是说明:

  1. σᵀ A σ 由于 Φ 的RIP性质,以概率1等于 1 ± O(ε)
  2. σᵀ B σσᵀ C σ 可以利用 σ 的随机性,通过Hanson-Wright不等式等工具,证明它们以高概率是 O(ε)

将这些部分组合起来,就得到了定理的结论。

在本次课程中,我们完成了对 σᵀ A σ 项的分析。它直接利用了RIP性质:因为 A 对应的是原向量 z 分块后的各部分,而 Φ 作用在每个(经过随机符号翻转的)块上时,输入的向量是 K-稀疏 的,所以RIP保证了其范数几乎不变。将这些块的范数平方相加,就得到了整体范数平方的近似。

本节课中我们一起学习了快速Johnson-Lindenstrauss变换的Ailon-Chazelle构造及其分析,并介绍了Krahmer-Ward定理这一从压缩感知RIP性质出发理解FJLT的深刻框架。下一讲我们将完成Krahmer-Ward定理的剩余证明部分。

014:Krahmer-Ward 证明收尾与近似矩阵乘法

在本节课中,我们将完成 Krahmer-Ward 定理的证明,并开始介绍如何使用素描技术加速线性代数计算,特别是近似矩阵乘法。

概述

本节课首先将完成上周开始的 Krahmer-Ward 定理证明。该证明包含一些有趣的思想,我们将在后续讨论压缩感知时用到。之后,我们将开始探讨随机化线性代数,即使用素描技术来加速诸如近似矩阵乘法、最小二乘回归、低秩近似和 K 均值聚类等计算。

Krahmer-Ward 证明收尾

上一节我们介绍了 Krahmer-Ward 分析的核心思想,即通过条件概率以不同顺序分析快速 Johnson-Lindenstrauss 变换矩阵 SHD。本节我们将完成证明中剩余部分的推导。

回顾与设定

回忆一下,我们试图证明对于任意单位范数向量 z,有:
||ΠDΣz||² ≈ 1 ± ε
其中 Π = SHD 是对角符号矩阵,Σ 是随机符号向量。我们假设 Π 满足 (ε, 2k)-RIP 性质,这意味着 Π 能同时保留下所有 2k-稀疏向量的范数。

我们将目标表达式 ||ΠDΣz||² 重写为二次型:
||ΠDΣz||² = Σᵀ X Σ
其中矩阵 X 被根据向量 z 的分块结构进行了划分。我们定义了四个子矩阵 A, B, Bᵀ, C,并希望证明:

  1. Σᵀ A Σ = 1 ± ε (已证明)
  2. Σᵀ B Σ = ± O(ε)
  3. Σᵀ Bᵀ Σ = ± O(ε)
  4. Σᵀ C Σ = ± O(ε)

通过联合界,这些项同时成立的概率很高,从而完成证明。

分析 Σᵀ B Σ

我们首先分析 Σᵀ B Σ。通过代数操作,可以将其表示为:
Σᵀ B Σ = vᵀ Σ₋₁
其中 v 是一个依赖于 Σ₁Σ 的前 k 个分量)的固定向量,而 Σ₋₁Σ 中剩余的分量,与 v 独立。

这是一个固定向量与随机符号向量的点积。根据 Khintchine 不等式,其尾部概率有界:
P(|vᵀ Σ₋₁| > ε) ≤ 2 exp(-ε² / (2||v||²))
因此,关键在于界定向量 v 的范数 ||v||

通过一系列推导,并利用 Π 的 RIP 性质以及称为 壳层法 (shelling) 的技术(在压缩感知中常用),我们可以证明:
||v|| ≤ O(ε / √k)

将其代入 Khintchine 不等式,并选择 k = Θ(log(1/δ)),可得失败概率最多为 δ/3

分析 Σᵀ Bᵀ Σ

由于 Σᵀ Bᵀ Σ 是一个实数,它等于其转置 (Σᵀ B Σ)ᵀ,而后者就是 Σᵀ B Σ。因此,对 Σᵀ B Σ 的界自动适用于此项。

分析 Σᵀ C Σ

最后分析 Σᵀ C Σ。这是一个更复杂的二次型。我们使用 Hanson-Wright 不等式 来界定其偏离概率。

我们需要界定矩阵 CFrobenius 范数算子范数。通过利用 Π 的 RIP 性质、算子范数与 Frobenius 范数之间的关系,以及壳层法,可以证明:

  • ||C||_F² ≤ O(ε² / k)
  • ||C||_{op} ≤ O(ε / √k)

将这两个界代入 Hanson-Wright 不等式,并再次利用 k = Θ(log(1/δ)),可得失败概率最多为 δ/3

证明完成

通过联合界,所有四项 (A, B, Bᵀ, C) 同时满足其误差界的概率至少为 1 - δ。这就完成了 Krahmer-Ward 定理的证明,表明 ΠDΣ 能以高概率实现分布式的 Johnson-Lindenstrauss 嵌入。

关于 RIP 矩阵的说明:证明中要求 Π 满足 RIP。虽然存在确定性的 RIP 矩阵构造,但已知的最佳结果需要大约 行。在实践中,通常使用随机矩阵(如高斯矩阵或快速 JL 矩阵),它们能以高概率满足 RIP,且行数约为 O(k log d)

转向线性代数的素描算法

上一节我们完成了对 Krahmer-Ward 定理的理论分析。现在,我们开始一个新的主题:使用素描技术加速线性代数计算。

近似矩阵乘法

第一个问题是 近似矩阵乘法。假设我们有矩阵 A ∈ ℝ^(d×n)B ∈ ℝ^(n×p),想要计算乘积 AᵀB。朴素算法需要 O(dnp) 时间。

是否存在更快的近似算法?已知的快速矩阵乘法算法(如 Strassen 及其后续改进)虽然能降低理论指数,但在实践中往往不实用,尤其是对于非方阵情况。

DKM 采样方法 (2006)

Drineas, Kannan, Mahoney 提出了一种基于采样的近似方法。核心思想是利用矩阵乘法的外积表示:
AᵀB = Σ_{i=1}^{n} a_i b_iᵀ
其中 a_iA 的第 i 列,b_iᵀB 的第 i 行。

以下是 DKM 方法的基本步骤:

  1. 非均匀采样:以概率 p_i 采样索引 ip_i 正比于对应行/列范数的乘积:p_i ∝ ||a_i|| * ||b_i||
  2. 构造估计量:进行 m 次独立采样(有放回),得到索引 i_1, ..., i_m。估计矩阵 C 定义为:
    C = (1/m) * Σ_{r=1}^{m} (1/p_{i_r}) * a_{i_r} b_{i_r}ᵀ
  3. 性质CAᵀB 的无偏估计。通过精心选择 p_im,可以证明以下误差界以高概率成立:
    ||AᵀB - C||_F ≤ ε * ||A||_F * ||B||_F
    其中所需样本数 m = O(1/(ε²δ)) 可确保失败概率为 δ

提升成功概率

上述方法中,失败概率 δ 直接出现在样本复杂度 m 中。若希望得到对数依赖(log(1/δ)),可以使用 中位数提升 (median boosting) 技巧。

然而,由于输出是矩阵而非标量,需要定义“矩阵中位数”。思路如下:

  1. 独立运行上述算法 T = O(log(1/δ)) 次,得到 T 个输出矩阵 C_1, ..., C_T
  2. 根据切尔诺夫界,其中至少 2T/3 个矩阵是“好”的(满足误差界)。
  3. 这些“好”的矩阵彼此接近(根据三角不等式)。
  4. 算法:遍历所有 T 个矩阵,找到第一个与大多数(> T/2 个)其他矩阵都接近的矩阵,将其作为输出。这可以在 O(T²) 次矩阵距离比较内完成。
  5. 或者,随机选取一个矩阵,检查它是否与大多数矩阵接近;若不是,则重复。期望只需常数次尝试,期望比较次数为 O(T)

基于 JL 的素描方法

与 DKM 同年,Charikar 提出了基于 Johnson-Lindenstrauss 变换的另一种方法。其核心是选择一个 JL 素描矩阵 Π ∈ ℝ^(m×n)(例如随机符号矩阵或快速 JL 矩阵),然后计算:
C = (ΠA)ᵀ (ΠB)
可以证明,以高概率有:
||AᵀB - (ΠA)ᵀ(ΠB)||_F ≤ ε * ||A||_F * ||B||_F
其原理在于,JL 性质不仅保范数,也保内积。误差矩阵的每个元素是 (Πa_i)ᵀ(Πb_j) - a_iᵀb_j,由于 JL 保证内积近似,这些误差的平方和(即 Frobenius 范数平方)也就很小。

我们将在下节课详细展开这种方法的分析。

总结

本节课中,我们一起完成了 Krahmer-Ward 定理的证明,该定理为快速 Johnson-Lindenstrauss 变换提供了严格的分析框架。随后,我们引入了使用素描技术加速线性代数计算的主题,并重点介绍了近似矩阵乘法的两种早期方法:DKM 的基于非均匀采样的方法和 Charikar 的基于 JL 变换的方法。这些方法为后续讨论更复杂的线性代数问题(如回归、低秩近似)奠定了基础。

015:JL矩性质、子空间嵌入与草图求解

在本节课中,我们将继续讨论随机化线性代数中的草图算法。我们将学习JL矩性质的定义,并了解如何利用它来分析近似矩阵乘法。接着,我们将引入一个更强大的概念——子空间嵌入,并探讨其在“草图求解”范式中的应用,特别是用于加速最小二乘回归问题。

JL矩性质与近似矩阵乘法

上一节我们介绍了基于采样的近似矩阵乘法方法,其误差以Frobenius范数衡量。该方法通过采样行并适当缩放来构造无偏估计量。然而,该方法依赖于矩阵A和B的行范数,在数据流等动态更新场景中可能不适用。

本节中,我们将探讨另一种基于Johnson-Lindenstrauss引理的方法。该方法通过一个随机投影矩阵Π来近似A和B,然后计算(ΠA)T(ΠB)来近似ATB。只要Π的分布良好且行数足够,就能以高概率保证Frobenius范数误差很小。

为了分析这种方法,我们首先引入一个核心概念。

定义:JL矩性质
一个分布D满足(ε, δ, p) JL矩性质,如果对于所有单位范数向量z,都有:
E_Π~D [ | ||Πz||_2^2 - 1 |^p ] ≤ ε^p * δ
其中p ≥ 1。

这个性质意味着,通过马尔可夫不等式,我们可以从矩的界推导出分布D能提供分布式的JL保证。我们之前讨论过的多种草图矩阵,如Count Sketch、稀疏JL、快速JL等,都满足特定参数下的JL矩性质。

接下来,我们看看这个性质如何帮助我们分析点积的近似。

引理:点积的矩界
假设分布D满足(ε, δ, p) JL矩性质。那么,对于所有单位范数向量x和y,有:
|| (Πx)^T(Πy) - x^Ty ||_p ≤ ε * δ^(1/p)

证明思路:
利用恒等式 x^Ty = 1/4 (||x+y||^2 - ||x-y||^2),将点积差表示为范数平方差的组合,然后应用三角不等式和JL矩性质即可得证。

这个引理表明,仅从保持向量范数的矩界,就能自动得到保持点积的相同质量的矩界,无需付出额外代价。

现在,我们可以回到近似矩阵乘法的分析。

定理:基于JL的近似矩阵乘法
假设分布D满足(ε, δ, p) JL矩性质,且p ≥ 2。那么,以至少1-δ的概率,有:
|| (ΠA)^T(ΠB) - A^TB ||_F ≤ ε * ||A||_F * ||B||_F

证明概要:

  1. 定义误差矩阵 M = (ΠA)^T(ΠB) - A^TB。其(i, j)元素为 (Πa_i)^T(Πb_j) - a_i^T b_j,其中a_i, b_j是A和B的列。
  2. 计算M的Frobenius范数的p次矩。通过归一化列向量并应用之前的点积引理,可以将问题转化为对一系列独立随机变量矩的求和。
  3. 利用三角不等式、柯西-施瓦茨不等式,最终将p次矩的界控制在 (ε * ||A||_F * ||B||_F * δ^(1/p))^p 以内。
  4. 应用马尔可夫不等式,即可得到关于Frobenius范数的概率界。

这个定理表明,任何满足JL矩性质的分布(对应我们已知的各种JL构造),都能自动提供近似矩阵乘法保证,且无需像朴素联合界论证那样为处理大量点积而将δ设置得非常小。

从算子范数到子空间嵌入

上一节我们讨论了Frobenius范数下的近似矩阵乘法。现在,我们考虑一个更强的误差度量——算子范数(谱范数)。特别地,我们关注B = A的情况,即近似计算A^TA。

我们希望:|| (ΠA)^T(ΠA) - A^TA ||_op ≤ ε * ||A||_op^2
根据定义,这等价于要求对于所有单位范数向量z,有:
| ||ΠAz||_2^2 - ||Az||_2^2 | ≤ ε * ||Az||_2^2

这引出了一个更强大且广泛应用的概念:子空间嵌入。

定义:子空间嵌入
对于一个线性子空间E ⊆ R^n,矩阵Π被称为E的ε-子空间嵌入,如果对于所有x ∈ E,都有:
(1 - ε) ||x||_2^2 ≤ ||Πx||_2^2 ≤ (1 + ε) ||x||_2^2

设U是一个n×d矩阵,其列构成E的一组标准正交基(即U^T U = I)。那么,Π是E的ε-子空间嵌入,当且仅当:
|| (ΠU)^T(ΠU) - I ||_op ≤ ε

这可以看作是分布JL引理在矩阵上的自然推广:从保持单个向量的范数,推广到保持一个子空间内所有向量的范数。

定义: oblivious子空间嵌入
一个分布D被称为(d, ε, δ) oblivious子空间嵌入,如果对于任意n×d矩阵U(满足U^T U = I),以至少1-δ的概率从D中抽取的矩阵Π满足:
|| (ΠU)^T(ΠU) - I ||_op ≤ ε

如何构造子空间嵌入

以下是几种构造子空间嵌入的方法:

1. 已知正交基
如果已知子空间E的标准正交基矩阵U,那么最简单的嵌入就是取Π = U^T。这将E中的向量x = Uz映射为z,完美保持范数(ε=0),并将维度降至d(这是可能的最小维度)。缺点是计算U(例如通过Gram-Schmidt或SVD)可能代价高昂。

2. 行采样
对于以矩阵A列空间形式给出的子空间,可以考虑对A的行进行采样。即构造采样矩阵Π,使得(ΠA)T(ΠA)是ATA的无偏或近似估计。关键问题在于如何设置采样概率。与Frobenius范数近似乘法不同,此时的最佳采样概率与A的“杠杆得分”有关,我们将在后续课程中详细讨论。

3. Oblivious子空间嵌入
这是我们本节重点介绍的方法。使用满足JL矩性质的随机矩阵Π(如Count Sketch、稀疏JL等)。通过更复杂的分析(通常涉及矩阵的矩不等式,如矩阵Bernstein不等式),可以证明当投影维度m约等于d/ε²时(可能带有对数因子),Π能以高概率成为任意d维子空间的嵌入。这种方法 oblivious,因为Π的选择不依赖于具体的子空间。

应用:草图求解范式

在介绍了子空间嵌入的概念后,我们来看一个重要的应用:“草图求解”范式。该范式由Sarlos在2006年提出,用于加速最小二乘回归等计算问题。

考虑最小二乘回归问题:β_ls = argmin_β ||Xβ - y||_2,其中X是n×d矩阵(n >> d)。其解析解为 β_ls = (X^T X)^† X^T y,其中†表示伪逆。直接计算需要O(nd²)时间。

草图求解法的思路是,通过一个子空间嵌入矩阵Π将问题“草图化”到低维空间,然后在草图空间中求解:

  1. 草图化:计算 X̃ = ΠXỹ = Πy,其中Π是一个m×n的随机矩阵,m远小于n。
  2. 求解:在草图空间中求解回归问题:β̃ = argmin_β ||X̃β - ỹ||_2。其解为 β̃ = (X̃^T X̃)^† X̃^T ỹ
  3. 时间复杂度:求解步骤耗时O(md²)。如果Π是稀疏的(如Count Sketch)或具有快速变换结构(如基于FFT的快速JL),则草图化步骤ΠX也可以快速完成,总时间远小于O(nd²)。

关键问题:为什么β̃是一个好的解?

定理:草图求解的近似保证
令E为X的列空间与向量y张成的子空间,其维度至多为d+1。如果Π是子空间E的一个ε-子空间嵌入,那么用β̃得到的回归残差满足:
||Xβ̃ - y||_2^2 ≤ (1+ε)/(1-ε) * ||Xβ_ls - y||_2^2

证明概要:

  1. 由于β̃是草图空间中最小化||ΠXβ - Πy||_2的解,因此有 ||ΠXβ̃ - Πy||_2^2 ≤ ||ΠXβ_ls - Πy||_2^2
  2. 因为对于任意β,向量Xβ - y都在子空间E中,而Π是E的嵌入,所以有:
    (1-ε) ||Xβ - y||_2^2 ≤ ||Π(Xβ - y)||_2^2 ≤ (1+ε) ||Xβ - y||_2^2
  3. 将上述嵌入性质应用于不等式两边,经过简单代数推导即可得到定理中的界。

这个定理表明,尽管我们求解的是一个压缩后的问题,但其解在原问题上的表现几乎与最优解一样好。通过结合子空间嵌入和近似矩阵乘法等技术,可以进一步优化所需的草图维度m,甚至可以将其降至O(d)(与ε无关),或用于为迭代法(如梯度下降)构造高效的预处理器。我们将在后续课程中深入探讨这些高级主题。

总结

本节课我们一起学习了:

  1. JL矩性质:一种通过矩来刻画分布JL保证的方式,它允许我们推导出点积和矩阵乘法的近似界,而无需使用浪费的联合界。
  2. 子空间嵌入:一个比JL更强的概念,要求随机投影能同时保持一个子空间内所有向量的范数。它是分析更复杂草图算法的核心工具。
  3. 草图求解范式:利用子空间嵌入将大规模优化问题(如最小二乘回归)压缩到一个小规模问题上求解,从而显著降低计算成本,同时提供理论上的近似保证。

下节课我们将继续探讨如何通过杠杆得分采样来构造子空间嵌入,并研究子空间嵌入在低秩近似、聚类等其他问题中的应用。

016:杠杆得分采样与无感知子空间嵌入构建

在本节课中,我们将继续探讨用于线性代数问题(如回归和低秩近似)的素描算法。我们将重点学习两种获取子空间嵌入的方法:杠杆得分采样和无感知子空间嵌入构建,并深入分析它们在回归问题中的应用。


子空间嵌入回顾

上一讲我们讨论了如何获取子空间嵌入,并提到了三种方法。其中一种理想情况是,如果我们已经拥有子空间的一个标准正交基矩阵 U,那么 U^T 本身就是最佳的子空间嵌入。然而,计算 U 的计算成本很高。因此,我们倾向于寻找计算成本更低的方法。

今天,我们将详细探讨两种这样的方法:

  1. 杠杆得分采样:一种基于行重要性采样的方法。
  2. 无感知子空间嵌入:使用与输入矩阵无关的随机矩阵分布。

杠杆得分采样

采样方法的核心思想是,通过按特定概率随机选取矩阵 A 的行,构造一个较小的矩阵 ΠA,使得 ΠA^T ΠA 能够很好地近似 A^T A

我们定义一个采样矩阵 Π。一种简单的构造是:独立地采样 M 行,每行在随机位置有一个非零元素 1/√p_i。更便于分析的一种等价形式是定义一个对角矩阵 Π,其第 i 个对角元素为 η_i / √p_i,其中 η_i 是一个伯努利随机变量,以概率 p_i 取值为1。这样,ΠA^T ΠA 的期望就是 A^T A

杠杆得分与敏感度

为了确定采样概率 p_i,我们引入敏感度的概念。第 ia_i 的敏感度 r_i 定义为:
r_i = sup_x ( (a_i^T x)^2 / ||A x||_2^2 )
这个比值衡量了第 i 行在矩阵 A 的二次型中的重要程度。

一个关键的发现是,敏感度 r_i 等于该行的杠杆得分 l_i。杠杆得分 l_i 定义为:
l_i = a_i^T (A^T A)^† a_i
其中 表示伪逆。从奇异值分解 A = U Σ V^T 的角度看,杠杆得分 l_i 就是矩阵 Ui 行范数的平方。

神奇之处在于,所有行的杠杆得分之和等于子空间的维度 d
∑_{i=1}^n l_i = ||U||_F^2 = d

采样定理

根据Spielman和Srivastava的定理,如果我们设置采样概率 p_i 满足:
p_i ≥ min(1, C * l_i * (log d / δ) / ε^2)
其中 C 是一个常数,那么以至少 1-δ 的概率,采样得到的 Π 是矩阵 A 列空间的一个 ε-子空间嵌入。

此时,期望的采样行数(即草图大小)为:
E[#rows] = ∑ p_i = O( (d log d / δ) / ε^2 )
这达到了接近最优的行数 O(d),仅多出一个对数因子。

实际应用中的关键点

你可能会问:计算杠杆得分 l_i 需要知道 U,而计算 U 本身就很昂贵,这岂不是矛盾?

关键在于,我们不需要精确的杠杆得分。只要我们的采样概率 p_i 是真实杠杆得分的一个常数倍近似,我们就可以稍微过度采样,并仍然能保证定理成立。事实上,存在快速算法(例如通过Johnson-Lindenstrauss变换)来近似计算所有杠杆得分,这将在课程作业中涉及。


无感知子空间嵌入

无感知子空间嵌入是指一个分布 Dm×n 矩阵上,使得对于任何 n×d 的标准正交基矩阵 U(满足 U^T U = I),从 D 中随机抽取的矩阵 Π 能以高概率满足:
|| (ΠU)^T (ΠU) - I ||_op ≤ ε
这意味着 Π 能几乎等距地嵌入 U 的列张成的子空间。

构造方法一:网论证

子空间单位球是无限的,无法直接对每个向量应用Johnson-Lindenstrauss引理。网论证的思路是:

  1. 在子空间的单位球上构造一个有限的 ε-网 E'。可以证明,网的大小最多为 (O(1/ε))^d
  2. 应用分布式的JL引理,确保以高概率保留网 E' 中所有向量的范数。
  3. 通过线性代数论证,可以证明这足以保证子空间中所有向量的范数都被保留。

由此所需的行数 m 为:
m = O( (d + log(1/δ)) / ε^2 )
这已被证明是最优的。标准的随机±1矩阵、稀疏JL矩阵或快速JL变换矩阵都可以作为满足条件的分布 D

构造方法二:矩方法

我们直接分析失败概率:
Pr[ || (ΠU)^T (ΠU) - I ||_op > ε ]
通过马尔可夫不等式,对于偶数 p,上式不超过:
E[ || (ΠU)^T (ΠU) - I ||_S_p^p ] / ε^p
其中 ||·||_S_p 是Schatten p-范数(即奇异值向量的 L_p 范数)。

关键观察:

  • p ≥ log d 时,Schatten p-范数与算子范数(即Schatten ∞-范数)仅相差一个常数因子。因此,用Schatten范数上界算子范数不会造成大的损失。
  • 对于实对称矩阵 M,当 p 为偶数时,||M||_S_p^p = Trace(M^p)
  • (ΠU)^T (ΠU) - I 展开为随机矩阵的和,然后可以利用非交换的Khintchine不等式来 bound Trace(M^p) 的期望。

这种方法可以用于分析稀疏JL矩阵(如Count Sketch)作为OSE的性能,有时能得到比网论证更紧的界。

构造方法三:近似矩阵乘法

注意到,若要求:
|| (ΠU)^T (ΠU) - I ||_F ≤ ε
由于Frobenius范数大于算子范数,这足以保证算子范数条件。而左边正是 || U^T U - (ΠU)^T (ΠU) ||_F,这正是一个近似矩阵乘法问题。

对于Count Sketch,已知其满足AMM性质:|| A^T B - (ΠA)^T (ΠB) ||_F ≤ ε ||A||_F ||B||_F
这里 A = B = U,且 ||U||_F^2 = d。因此,我们需要设置AMM的误差参数为 ε' = ε / d
Count Sketch 达到此误差所需行数为:m = O( 1 / (ε')^2 ) = O( d^2 / ε^2 )
这虽然行数较多,但 ΠA 的乘法计算极快,时间复杂度与 A 的非零元个数成比例。


回归问题的进一步分析

上一讲我们分析了“素描-求解”方法:若 Π 是增广矩阵 [X, y] 列空间的 ε-子空间嵌入,则通过求解草图系统 min_β ||ΠX β - Πy|| 得到的解 β̃ 满足:
||X β̃ - y|| ≤ (1+O(ε)) * OPT

更精细的分析

本节介绍一种更精细的分析方法,它只要求 ΠX 列空间的一个常数精度(如 1/4)子空间嵌入,并结合近似矩阵乘法性质。

设定

  • β* 为原问题最优解,β̃ 为草图问题最优解。
  • w = y - X β*,则 OPT = ||w||^2
  • X 进行SVD:X = U Σ V^T
  • 定义 α 使得 X β* = U α
  • 定义 γ 使得 X β̃ - X β* = U γ

目标:证明 ||X β̃ - y||^2 = OPT + ||γ||^2,并 bound ||γ||^2

推导

  1. 利用 β̃ 在草图系统中的最优性,可以推导出关键等式:
    (ΠU)^T (ΠU) γ = (ΠU)^T (Π w)
  2. 由于 ΠX 列空间(即 U 的列空间)的常数精度子空间嵌入,矩阵 (ΠU)^T (ΠU) 的最小奇异值离1不远。因此:
    || (ΠU)^T (ΠU) γ || ≥ Ω(1) * ||γ||
  3. 结合步骤1的等式,有:
    Ω(1) * ||γ|| ≤ || (ΠU)^T (Π w) || = || (ΠU)^T (Π w) - U^T w ||
    (因为 w 正交于 U 的列空间,所以 U^T w = 0
  4. 右边正是近似矩阵乘法的形式。利用AMM性质,并设置合适的误差参数,可以证明:
    || (ΠU)^T (Π w) - U^T w || ≤ O(√ε) * ||U||_F * ||w|| = O(√(ε * d)) * ||w||
  5. 综合起来,得到 ||γ||^2 ≤ O(ε) * ||w||^2 = O(ε) * OPT

因此,最终误差为 OPT + O(ε) * OPT = (1 + O(ε)) * OPT。这种分析有时允许使用更小的草图维度。

预处理迭代方法

最后简要提及另一种思路:使用子空间嵌入构造预条件子,然后结合迭代法(如梯度下降)来求解回归问题。这种方法通常不直接求解草图系统,而是利用草图来加速原问题的迭代求解过程。我们将在下一讲中继续讨论。


总结

本节课我们深入学习了两种构建子空间嵌入的实用技术:

  1. 杠杆得分采样:通过按与行重要性(杠杆得分)成比例的概率采样行,可以高效地构建子空间嵌入,所需行数接近最优的 O(d)
  2. 无感知子空间嵌入:通过网论证、矩方法或近似矩阵乘法性质,可以设计和分析多种随机矩阵(如稠密高斯矩阵、稀疏JL矩阵、Count Sketch)作为适用于任意子空间的嵌入。

我们还看到了如何将这些工具应用于回归问题的“素描-求解”框架,并学习了两种不同的分析技巧:一种基于严格的子空间嵌入性质,另一种结合常数精度嵌入和近似矩阵乘法,可能获得更好的参数。

在接下来的课程中,我们将继续探索这些线性代数素描算法在低秩近似等更多问题中的应用。

017:更快的迭代回归与低秩近似

在本节课中,我们将学习两种利用草图技术加速线性代数计算的方法:基于预处理的迭代回归求解低秩矩阵近似。我们将看到,这两种方法的核心分析技巧与我们之前讨论的回归问题分析紧密相关。

概述:两种加速方法

上一节我们介绍了分析草图最小二乘的两种方法。本节中,我们将探讨第三种方法:使用草图技术寻找一个好的预处理器,然后将其用于梯度下降等迭代求解器。接着,我们将讨论如何将草图技术应用于低秩近似问题。我们会发现,其分析过程几乎完全复用了回归问题中已讨论过的技巧。

第一部分:基于预处理的迭代回归

梯度下降与最小二乘

首先,回顾梯度下降法。对于最小化函数 F(x) 的问题,我们从初始猜测 x0 开始,迭代更新:
x_{t+1} = x_t - γ * ∇F(x_t)
其中 γ 是一个小的步长参数。

对于最小二乘问题,我们的目标是:
min_β ||Xβ - y||_2^2
其梯度为 ∇F(β) = X^T(Xβ - y)。若设置步长 γ=1,则更新规则为:
β_{t+1} = β_t - X^T(Xβ_t - y)

收敛性分析的关键

对上述迭代过程进行分析后,我们发现收敛速度与矩阵 X条件数密切相关。如果 X 的所有奇异值都接近1(即条件数很小),则算法会快速收敛。

然而,一般的输入矩阵 X 可能不具备这个性质。我们的解决方案是:使用草图技术对 X 进行预处理,使其条件数变好

使用草图作为预处理器

以下是具体步骤:

  1. 计算草图:选择一个常数质量的子空间嵌入矩阵 Π(例如,1/3-子空间嵌入)。计算 ΠX
  2. 计算预处理矩阵:对 ΠX 进行SVD分解:ΠX = U' Σ' V'^T。定义预处理矩阵 R = V' (Σ')^{-1}
  3. 转换问题:考虑等价问题 min_β ||(XR)β - y||_2^2。关键性质是,矩阵 XR 的所有奇异值都接近1。
  4. 迭代求解:在新的矩阵 XR 上运行梯度下降。由于 XR 条件良好,算法会快速收敛。
  5. 初始化:为了获得初始解 β0,我们可以对原问题 min_β ||Xβ - y||_2^2 运行一次草图求解,只需常数精度(如 ε=0.1)即可。这提供了满足 ||Xβ0 - y||^2 ≤ 1.1 * OPT 的起点。

算法优势

  • 运行时间:主要开销是每次迭代中与稀疏矩阵 X 的乘法(O(nnz(X))),加上应用小矩阵 R 的代价(O(d^2))。总运行时间为 O((nnz(X) + d^2) * log(1/ε)),实现了对高精度 ε 的高效求解。
  • 与草图求解对比:传统的草图求解方法运行时间依赖于 poly(1/ε),而本方法仅依赖于 log(1/ε),在需要高精度解时优势显著。

上一部分我们利用草图改善了矩阵的条件数,从而加速了迭代求解。接下来,我们将目光转向另一个重要问题——低秩近似。

第二部分:基于草图的低秩近似

问题定义与经典解

给定矩阵 A ∈ R^{n×d} 和参数 k,我们希望找到秩不超过 k 的矩阵 B,以最小化 AB 之间的差异(通常使用Frobenius范数或算子范数)。

根据Eckart–Young定理,最优解 A_k 可通过截断 A 的SVD分解得到:A_k = U_k Σ_k V_k^T,其中只保留前 k 个最大的奇异值及对应的奇异向量。精确计算SVD的代价约为 O(n * d^2)

草图加速方案

我们不直接计算全SVD然后截断,而是采用以下近似方案:

  1. 计算草图矩阵 ΠA^T,其中 Π ∈ R^{m×n}m 远小于 n
  2. P 为到 ΠA^T 列空间的正交投影矩阵。
  3. 计算 PA 的最佳秩 k 近似,作为我们的近似解 Ã_k

即:
Ã_k = argmin_{rank(B) ≤ k} ||PB - A||_F^2
这等价于在 ΠA^T 的列空间内寻找 A 的最佳秩 k 近似。

算法效率

直接按照定义计算 Ã_k 的步骤如下:

  1. 计算 ΠA^T 的SVD:ΠA^T = U' Σ' V'^TU' 的列构成了 ΠA^T 列空间的一组标准正交基。
  2. 计算矩阵 C = U'^T A。这是一个 m × d 的矩阵。
  3. 计算 C 的最佳秩 k 近似 C_k(例如,通过计算 C 的SVD并截断)。
  4. 输出 Ã_k = U' C_k

运行时间主要消耗在计算 U'^T AO(m * n * d))和计算 C 的SVD(O(d * m^2))。通过精心选择草图,我们可以使 m 约为 O(k/ε)O(k^2),从而相比全SVD的 O(n * d^2) 有显著加速。

更进一步的优化是,将步骤2和3中的精确最小化,替换为使用草图求解技术来近似求解一个相关的回归问题,从而将运行时间中的 d 也替换为 k

近似质量保证

Ã_k 的质量取决于草图矩阵 Π 需满足的两个性质:

  1. 子空间嵌入Π 需是某个特定 k 维子空间的常数质量(如 1/2)子空间嵌入。
  2. 近似矩阵乘法:对于Frobenius范数,Π 需提供误差约为 O(√(ε/k)) 的近似矩阵乘法保证。

满足这些条件的草图(如Count Sketch、高斯矩阵等)可以确保:
||A - Ã_k||_F^2 ≤ (1 + ε) * ||A - A_k||_F^2

证明思路概要

证明的关键是将误差 ||A - Ã_k||_F^2 分解为两部分:
||A - Ã_k||_F^2 = ||(I - P)A_k||_F^2 + ||(I - P)(A - A_k)||_F^2
其中,P 是到 ΠA^T 列空间的投影。

  • 第二部分 ||(I - P)(A - A_k)||_F^2 很容易被 ||A - A_k||_F^2(即最优误差OPT)所界定。
  • 核心在于控制第一部分 ||(I - P)A_k||_F^2。通过展开和变形,可以发现这一项可以写成一系列向量差异的平方和。每个向量差异项,都恰好对应一个以 A_k^T 为设计矩阵、以 A^T 的某一列为响应向量的最小二乘回归问题的草图求解误差

因此,我们可以直接套用之前(Charikar风格)分析草图求解回归误差的技术:依赖于某个特定子空间(此处是 V_k 的列空间)的子空间嵌入性质,以及近似矩阵乘法性质。最终可以证明,在 Π 满足前述两个条件时,第一部分误差可以被控制在 ε * OPT 以内。

总结

本节课中我们一起学习了:

  1. 迭代回归的加速:通过使用草图技术构建预处理器 R,将原矩阵 X 转换为条件数良好的 XR,从而使得梯度下降法等迭代算法能够快速收敛到高精度解。该方法运行时间仅对数依赖于目标精度 ε
  2. 低秩近似的加速:通过将矩阵 A 投影到草图矩阵 ΠA^T 的列空间,并在该子空间内寻找最佳秩 k 近似,我们可以高效地得到一个质量有保证的近似解 Ã_k。其分析巧妙地归结到了我们已经熟悉的回归问题草图求解分析框架中。

这两种方法展示了草图技术在优化和矩阵计算中的强大与灵活,其核心思想都是通过随机投影降低问题维度或改善问题性质,从而在保证近似精度的前提下大幅提升计算效率。

018:投影成本保持草图 (PCP Sketches) 📊

在本节课中,我们将学习一种称为投影成本保持草图的技术。这是一种用于加速低秩近似的草图方法,并且其应用范围更广,例如可以应用于K均值聚类。我们也会在今天讨论这一点。

概述

投影成本保持草图,有时简称为PCP草图,是Cohen、Eldar、Musco和Musco在2015年引入的概念。它允许我们通过“草图求解”范式来解决低秩近似问题:首先将高维问题投影到低维空间,在低维空间中求解,然后将找到的解提升回原始空间,并保证该解对于原始问题是近似最优的。

核心定义

A 是一个给定的矩阵。我们说矩阵 Π 是矩阵 A 的一个 ε-K 投影成本保持草图,如果对于所有秩至多为 K 的正交投影矩阵 P,都有以下不等式成立:

(1 - ε) * ||(I - P)A||_F^2 ≤ ||(I - P)AΠ^T||_F^2 ≤ (1 + ε) * ||(I - P)A||_F^2

这里,||·||_F 表示弗罗贝尼乌斯范数。直观地说,这意味着用 ΠA 的行进行草图化(即降维)后,任何秩-K投影所造成的“误差成本”都被保留在 (1 ± ε) 的因子内。

草图求解的工作原理

上一节我们介绍了PCP草图的核心定义,本节中我们来看看它如何实现“草图求解”。

假设 ΠA 的一个 ε-K PCP草图。我们考虑草图化后的问题:找到最优的秩-K投影 ,以最小化 ||(I - P̃)AΠ^T||_F^2。设 P 是原始问题 ||(I - P)A||_F^2 的最优解。

利用PCP草图的性质,我们可以推导出:

||(I - P̃*)A||_F^2 ≤ (1+ε)/(1-ε) * ||(I - P*)A||_F^2

这意味着,在草图空间中找到的最优投影 *,对于原始问题也是一个接近最优的解(误差在常数因子内)。这个论证过程与我们在回归问题中看到的“草图求解”分析类似。

扩展到约束低秩近似与K均值聚类

PCP草图不仅适用于无约束的低秩近似,还能处理约束低秩近似问题。设 C_K 是所有秩-K正交投影矩阵集合 P_K 的一个子集。我们希望找到:

P* = argmin_{P ∈ C_K} ||(I - P)A||_F^2

如果 Π 是一个PCP草图,那么它对于 C_K 中的所有投影矩阵都保持成本。因此,在草图空间中求解约束问题得到的解 *,对于原始约束问题也是近似最优的。

以下是K均值聚类如何被表述为约束低秩近似问题:

  1. 问题重述:给定n个d维数据点,K均值聚类的目标是最小化每个点到其最近聚类中心的平方距离之和。
  2. 与投影的联系:可以证明,任何K分区(聚类)的成本,都等于 ||A - X_P X_P^T A||_F^2,其中 A 的行是数据点,X_P 是一个根据分区P构造的标准化指示矩阵(其列是标准正交的)。而 X_P X_P^T 正是一个秩-K正交投影矩阵。
  3. 结论:因此,K均值聚类等价于在特定结构的秩-K投影矩阵集合(对应所有可能的分区)上,最小化 ||(I - P)A||_F^2。这是一个约束低秩近似问题。拥有PCP草图 Π 后,我们就可以在草图化的数据上运行K均值算法,找到的聚类方案对原始数据也是近似最优的。

关于近似算法的重要说明:即使原始问题(如K均值)是NP难的,我们无法精确求解草图化后的问题,但只要我们在草图空间中使用一个α-近似算法找到解 ,那么该解对于原始问题的近似比最多是 α * (1+ε)。这意味着近似最优性可以从草图问题“提升”到原始问题。

PCP草图的存在性与构造定理

那么,如何构造一个PCP草图呢?Cohen等人给出了一个定理,列出了确保 Π 成为PCP草图的几个充分条件。

定理:设 A = U Σ V^T 是矩阵 A 的奇异值分解。令 A_K = U_K Σ_K V_K^T 为其最优秩-K近似,A_{-K} = A - A_K。如果草图矩阵 Π 满足以下四个性质,那么它是 A 的一个 (O(ε)-K) PCP草图:

  1. 子空间嵌入ΠA_K 的行空间的 ε-子空间嵌入。这需要草图维度约为 O(K/ε^2)(对于高斯随机矩阵)。
  2. 范数保持| ||A_{-K} Π^T||_F^2 - ||A_{-K}||_F^2 | ≤ ε * ||A_{-K}||_F^2
  3. 近似矩阵乘法 (AMM) 1||A_{-K} A_{-K}^T - (A_{-K}Π^T)(Π A_{-K}^T)||_F ≤ (ε/√K) * ||A_{-K}||_F^2
  4. 近似矩阵乘法 (AMM) 2||V_K^T A_{-K}^T - (Π V_K)^T (Π A_{-K}^T)||_F ≤ (ε/√K) * ||A_{-K}||_F

关键点:性质1要求 Π 是一个** oblivious 子空间嵌入**。性质2、3、4可以通过要求 Π 满足 Johnson-Lindenstrauss (JL) 矩性质 来保证。具体来说:

  • 性质2是JL性质对矩阵范数的直接推广(见下方引理)。
  • 性质3和4是近似矩阵乘法的要求,而JL矩性质也能推出AMM。

因此,如果一个分布生成的随机矩阵 Π 同时是 oblivious 子空间嵌入 并 满足 JL 矩性质,那么以高概率,从该分布中抽取的 Π 将是任意固定矩阵 A 的PCP草图。

一个有用的引理:如果 Π 来自满足 (ε, δ)-p JL矩性质的分布,那么对于任意矩阵 M,有:

Pr[ | ||Π M||_F^2 - ||M||_F^2 | > ε * ||M||_F^2 ] ≤ δ

这个引理说明,JL性质不仅保持向量范数,也能自动保持矩阵的弗罗贝尼乌斯范数。

实例

  • 高斯随机矩阵:当草图行数 m = O(K/ε^2) 时,可满足条件。
  • Count Sketch 矩阵:当草图行数 m = O(K^2/ε^2) 时,可满足条件。

定理证明思路

上一节我们陈述了定理,本节中我们简要勾勒其证明思路,以理解四个条件如何共同保证PCP性质。

证明目标是:对所有秩-K投影 P,验证PCP定义中的不等式。令 Y = I - P

  1. 分解:将 A 分解为 A_KA_{-K}。目标表达式 ||Y A Π^T||_F^2||Y A||_F^2 可以展开为四项之和(分别涉及 A_KA_{-K} 的纯项和交叉项)。
  2. 误差项:定义四个误差矩阵 E1, E2, E3, E4,它们代表了草图化版本与原始版本在对应项上的差异。我们需要证明 |Trace(Y E_i Y)| 对于每个 i 都是 O(ε) * Trace(Y C Y),其中 C = A A^T
  3. 逐项控制
    • E1 (涉及 A_K):利用条件1(子空间嵌入),因为 A_K^T Y 的列在 A_K 的行空间中。
    • E2 (涉及 A_{-K}):利用条件2(范数保持)和条件3(AMM 1),结合投影矩阵 P 的秩为K这一事实,通过柯西-施瓦茨不等式进行放缩。
    • E3 和 E4 (交叉项):这是证明中最技巧性的部分。需要利用 E3 的列在 A 的列空间中的性质,引入伪逆 C^+,并巧妙地应用矩阵的柯西-施瓦茨不等式。最终,误差可以转化为涉及 Π V_KΠ A_{-K}^T 的弗罗贝尼乌斯范数,这正是条件4(AMM 2)所控制的对象。
  4. 合并:将四个误差项 bound 相加,即可得到总的误差在 O(ε) * Trace(Y C Y) 以内,从而完成证明。

详细的代数推导包含大量符号运算,但核心思想是利用给定的四个条件,分别控制来自最优近似部分 (A_K)、残差部分 (A_{-K) 以及两者交叉相互作用所产生的误差。

总结

本节课中我们一起学习了投影成本保持草图。我们首先了解了它的定义:一种能保持任何秩-K投影所产生误差成本的草图。接着,我们看到了它如何支持“草图求解”范式,用于低秩近似问题,并且可以扩展到约束设置,特别是K均值聚类问题。然后,我们介绍了确保PCP草图存在的关键定理,该定理要求草图矩阵同时满足子空间嵌入JL矩性质(后者可推出范数保持和近似矩阵乘法)。最后,我们概述了定理的证明思路,展示了如何通过分解和逐项控制来达成目标。PCP草图是线性代数草图化中一个基础而强大的工具,它将随机投影的理论保证与重要的机器学习优化问题紧密联系了起来。

019:P19 K-means, 压缩感知, RIP, 基追踪

在本节课中,我们将要学习如何将素描技术应用于K-means聚类问题,并开始探讨压缩感知的核心概念,包括受限等距性质(RIP)和基追踪算法。

K-means 与 PCP 素描

上一节我们介绍了线性代数中的素描技术。本节中我们来看看如何将这些技术应用于K-means聚类问题。

K-means可以表示为一个带约束的低秩近似问题。我们之前介绍了投影成本保持(PCP)素描的概念。PCP素描本质上就是子空间嵌入和JL矩阵,用于投影问题,从而为近似求解带约束的最小二乘问题提供一种“素描-求解”方法。

这是一个约束低秩近似的特例,意味着我们可以通过PCP素描将维度降至O(k/ε²)来求解。值得注意的是,这个界限完全不依赖于数据点的数量N,而只依赖于聚类数K。如果K是常数,那么这是一个很好的结果。

实际上,有更新的理论表明,可以进一步将素描维度降至O(log k/ε²)。这个分析更为复杂,但核心思想是利用JL性质。

一个常数因子近似证明

以下定理展示了如何使用仅O(log k)维的素描来获得K-means的一个常数因子近似解。

定理:给定数据矩阵A、聚类数k和参数ε。设Π是一个将数据行从d维映射到m维的素描矩阵,满足特定条件(稍后详述)。令P*为原始问题的最优聚类投影矩阵。令P̃为在素描后问题上找到的γ近似最优解,即满足:
||(I - P̃)AΠᵀ||_F² ≤ γ * min_P ||(I - P)AΠᵀ||_F²
那么,只要Π满足所需条件,P̃将是原始问题的一个(9 + O(ε))γ近似解,即:
||(I - P̃)A||_F² ≤ (9 + O(ε))γ * ||(I - P*)A||_F²
其中,所需的素描维度m可以是O(log k / ε²)。

证明思路

  1. 分解:将数据矩阵A分解为最优聚类部分和残差部分:A = P*A + (I - P*)A。记B = P*AB̄ = (I - P*)A。注意,B只有k个不同的行(即k个质心)。
  2. 应用三角不等式:目标量||(I - P̃)A||_F可以分解为||(I - P̃)B||_F + ||(I - P̃)B̄||_F
  3. 处理残差部分:由于投影不增加Frobenius范数,||(I - P̃)B̄||_F ≤ ||B̄||_F,而||B̄||_F就是最优K-means成本的平方根。
  4. 处理聚类部分:关键观察是B只有k个不同的行。因此,如果我们要求Π能保持这k个点之间的所有成对距离(即JL性质),那么只需要m = O(log k / ε²)维。这样,||(I - P̃)B||_F就可以用其素描版本的范数来近似。
  5. 结合条件:我们还需要Π满足JL矩性质,以保持固定矩阵B̄的范数。
  6. 代入与化简:利用P̃是素描后问题的γ近似解这一事实,以及P*在素描后问题上成本不劣于其原始成本(经过Π保持),最终通过一系列三角不等式和范数比较,可以得到所述界限。

这个证明的精妙之处在于利用了最优聚类矩阵B只有k个不同行这一事实,从而将对n个点的JL要求降低为对k个点的JL要求。


压缩感知简介

现在,让我们从K-means转向压缩感知。压缩感知的核心思想是,如果一个信号在某个已知基下是近似稀疏的,我们可以用少量的线性测量来近似恢复它。

具体来说,假设我们有一个高维向量x ∈ ℝⁿ,它是近似k稀疏的(即存在一个k稀疏向量z,使得x - z的某种范数很小)。我们希望找到一个恢复算法,给定线性测量值Πx(其中Π是m × n的测量矩阵,m << n),能输出一个估计值x̃,使得x - x̃很小。

精确稀疏恢复与 Vandermonde 矩阵

如果x是严格k稀疏的,一个经典方法是使用Vandermonde矩阵作为Π。只要取m ≥ 2k个测量值,并选择不同的采样点,就能保证任何2k列组成的子矩阵都是满秩的。恢复时,可以遍历所有可能的k列支撑集,求解线性系统,找到那个能产生k稀疏解的唯一支撑集。虽然解码是指数时间的,但在有限域上存在多项式时间算法(如伴随解码),在实数域上也有Prony方法。

近似稀疏恢复与 Count-Min/Count Sketch

在实际中,信号往往是近似稀疏而非严格稀疏。我们之前学过的Count-Min Sketch和Count Sketch算法实际上已经提供了非均匀的压缩感知方案。

  • Count-Min Sketch:通过O(k log n)次测量,可以保证估计值x̃满足||x - x̃||_∞ ≤ (1/k) * ||x_{tail}||_1。如果进一步将x̃投影到其最大的O(k)个分量上,还能得到||x - x̃||_1 ≤ (1+ε) * ||x_{tail}||_1的保证。
  • Count Sketch:通过O(k log n)次测量,可以保证||x - x̃||_∞ ≤ (1/√k) * ||x_{tail}||_2,同样可以通过投影得到L2范数下的保证。

这些方案的优点是,如果x确实是k稀疏的,则恢复误差为零。

均匀 vs 非均匀保证

在压缩感知文献中,有两种常见的保证类型:

  • 非均匀(For-each)保证:对于一个固定的信号x,测量矩阵Π和恢复算法以高概率(≥ 1-δ)成功恢复。Count-Min/Count Sketch属于此类。
  • 均匀(For-all)保证:存在一个固定的测量矩阵Π和恢复算法,使得对于所有可能的信号x,都能以高概率成功恢复。这是一种更强的保证。

接下来,我们将探讨如何获得均匀保证。


受限等距性质 (RIP) 与基追踪 (Basis Pursuit)

为了获得均匀保证,我们需要引入一个关键概念:受限等距性质。

定义:一个矩阵Π满足(k, ε)-RIP,如果对于所有k稀疏向量z,都有:
(1 - ε) ||z||_2² ≤ ||Πz||_2² ≤ (1 + ε) ||z||_2²
这意味着Π近似地保持所有k稀疏向量的L2范数。等价地,Π是所有由k个标准基张成的子空间的(1±ε)子空间嵌入。

如何获得RIP矩阵?

  1. 通过JL引理构造:随机高斯矩阵、随机符号矩阵等满足分布JL引理,可以证明它们以高概率满足RIP,所需行数m = O(k log(n/k) / ε²)。这是最常用的方法。
  2. 通过不相干矩阵构造:如果矩阵Π的列都是单位范数,且任意两列的内积绝对值(相干性)最多为α,那么当α ≤ ε/(k-1)时,Π满足(k, ε)-RIP。使用Reed-Solomon码可以构造这样的矩阵,但需要m = O(k²/ε²)行,效率较低。
  3. 显式构造:这是一个重要的开放问题。目前最好的显式构造需要m = O(k^{2 - ο(1)})行,并且构造非常复杂,涉及解析数论。

基追踪算法

给定RIP测量矩阵Π和测量值y = Πx,一个自然的恢复算法是求解以下L0最小化问题:
min ||z||_0 subject to Πz = y
其中||·||_0表示非零元素个数。这能找到最稀疏的相容解。然而,这个问题是NP难的。

一个实用的松弛方法是求解L1最小化问题,称为基追踪
min ||z||_1 subject to Πz = y
这是一个线性规划问题,可以在多项式时间内求解。

RIP保证基追踪成功

定理 (Candes, 2008):如果测量矩阵Π满足(2k, ε)-RIP,且ε < √2 - 1 ≈ 0.414,那么对于任意信号x,基追踪的解x̃满足以下误差界:
||x - x̃||_2 ≤ O(1/√k) * ||x_{tail}||_1
其中x_{tail}是x去掉其最大的k个分量(按绝对值)后的部分。

证明思路(简述)

  1. 定义h = x̃ - x为恢复误差。由于Πx̃ = Πx,所以Πh = 0,即h在Π的零空间中。
  2. 将坐标集按x绝对值大小分成块:T0是x中最大的k个坐标的集合,T1是h(在T0补集上)中最大的k个坐标的集合,T2是下一个k个坐标的集合,依此类推。
  3. 目标是通过RIP性质来界定h的范数。关键步骤是证明一个引理:∑_{j≥2} ||h_{T_j}||_2 ≤ (2/√k) * ||x_{tail}||_1 + ||h_{T_0 ∪ T_1}||_2。这个证明使用了“壳层”技巧,类似于证明Cramer-Rao界时的方法。
  4. 利用Πh = 0,将Π h_{T_0 ∪ T_1}表示为-∑_{j≥2} Π h_{T_j}
  5. 计算||Π h_{T_0 ∪ T_1}||_2²,并利用RIP性质将其与||h_{T_0 ∪ T_1}||_2²关联起来。同时,将点积<Π h_{T_0 ∪ T_1}, Π h_{T_j}>用原始向量点积(为零)加上RIP引入的误差来界定。
  6. 将第3步的引理代入,最终得到一个关于||h_{T_0 ∪ T_1}||_2的不等式,通过代数整理即可得到定理中的误差界。常数√2 - 1保证了整理过程中分母为正。

这个定理表明,只要测量矩阵满足适当的RIP条件,基追踪这个多项式时间算法就能为所有信号提供强大的、均匀的恢复保证。


总结

本节课中我们一起学习了:

  1. K-means的素描:利用PCP素描和最优聚类矩阵只有k个不同行的特性,可以将数据投影到O(log k/ε²)维并仍获得常数因子近似解。
  2. 压缩感知基础:从精确稀疏恢复(Vandermonde矩阵)到近似稀疏恢复(Count-Min/Count Sketch),区分了非均匀和均匀恢复保证。
  3. RIP与基追踪:引入了受限等距性质作为获得均匀保证的关键工具。证明了只要测量矩阵满足(2k, ε)-RIP(ε足够小),基追踪算法就能稳定地恢复近似稀疏信号,误差界由信号尾部能量的L1范数决定。

这为我们处理高维稀疏信号提供了一套强有力的理论工具和实用算法框架。

020:P20 迭代硬阈值、扩展图与RIP1

在本节课中,我们将继续讨论压缩感知,特别是使用RIP矩阵的“对所有”压缩感知。我们将介绍一种比求解线性规划更快的解码算法——迭代硬阈值法,并分析其性能保证。最后,我们还将简要介绍另一种基于扩展图的压缩感知方法。

回顾与引入

上一节我们讨论了如何使用RIP矩阵进行压缩感知。我们提到,如果矩阵满足RIP性质,那么通过求解基追踪(Basis Pursuit)线性规划问题,我们可以从测量值中恢复出原始稀疏向量 x,其误差由 x 的尾部 L1 范数控制。

然而,求解线性规划需要 O(poly(m, n)) 的时间,对于高维数据(例如百万像素的图像)来说可能太慢。因此,我们需要寻找更快的解码算法。

本节我们将重点介绍一种迭代算法——迭代硬阈值法,它能以更快的速度实现近似恢复。

迭代硬阈值法

迭代硬阈值法是一种快速的压缩感知解码算法。其核心思想是通过迭代更新,逐步逼近原始稀疏信号 x

以下是算法的伪代码描述:

输入:测量值 y = Φx + e, 测量矩阵 Φ, 稀疏度 k, 迭代次数 T
初始化:x_0 = 0
for t = 0 to T-1:
    a_{t+1} = x_t + Φ^T (y - Φ x_t)
    x_{t+1} = H_k(a_{t+1})  # H_k 是硬阈值算子,保留 a_{t+1} 中幅值最大的 k 个元素,其余置零
输出:x_T

算法直观理解:理想情况下,我们希望 Φ^T Φ 接近单位矩阵。那么,a_{t+1} ≈ x_t + (x - x_t) = x。因此,每次迭代都在向真实信号 x 靠近。硬阈值操作 H_k 是为了保证每次迭代的中间结果 x_t 都是 k-稀疏的,这样我们才能利用矩阵 Φ 在稀疏向量上的RIP性质进行分析。

定理与性能保证

以下是迭代硬阈值法的主要定理(由Blumensath, Davies 等人于2009年提出):

定理:假设测量矩阵 Φ 满足 (ε, 3k)-RIP 性质,且 ε < 1/(4√2)。那么,对于所有 T ≥ 1,算法输出的 x_T 满足:

||x_T - x||_2 ≤ (1/2)^T * ||x||_2 + O( (1/√k) * ||x - x_k||_1 + ||e||_2 )

其中 x_kx 的最佳 k-项近似,e 是测量噪声。

定理含义

  1. 指数收敛:误差项 (1/2)^T * ||x||_2 表明,随着迭代次数 T 增加,算法会指数级收敛到真实信号(在初始猜测为0的前提下)。
  2. 近似稀疏误差:项 O( (1/√k) * ||x - x_k||_1 ) 与基追踪的误差界类似,处理信号并非严格 k-稀疏的情况。
  3. 噪声鲁棒性:项 O( ||e||_2 ) 表明算法对测量噪声是稳健的。

上一节我们介绍了算法和定理,本节我们来看看证明的核心思路。

证明思路简述

证明的关键是定义残差 r_t = x - x_t,并推导其范数如何随着迭代衰减。

  1. 简化问题:通过将信号的尾部 x - x_k 归入测量噪声 e,可以假设 x 本身是 k-稀疏的。这是因为RIP矩阵不会过度放大任何向量的范数。
  2. 关键步骤:分析残差 r_{t+1} 的范数。通过代数展开和三角不等式,可以将其表示为前一轮残差 r_t 和噪声 e 的函数。
  3. 利用RIP性质:由于每次迭代的 x_t 都是 k-稀疏的,且我们考虑了其支撑集与真实信号支撑集的并集(大小至多 2k3k),因此可以应用矩阵 Φ 的RIP性质来约束某些算子范数。
  4. 得到递归不等式:最终可以得到形如 ||r_{t+1}||_2 ≤ (1/2) * ||r_t||_2 + C * ||e||_2 的递归关系。
  5. 求解递归式:通过展开递归式,即可得到定理中所示的误差界,其中包含指数衰减项和常数倍的噪声项。

另一种方法:基于扩展图的压缩感知

除了使用满足RIP2(保持L2范数)的随机矩阵(如通过JL引理构造),还存在一种完全不同的压缩感知方法,它基于组合对象——扩展图。

扩展图:一个二分扩展图 G=(U, V, E),其中 |U|=n(对应信号维度),|V|=m(对应测量数),每个左节点 u∈U 的度数为 d。它是 (k, ε, d)-扩展器,如果对于任意大小至多为 k 的左节点子集 S,其邻居集 Γ(S) 的大小至少为 (1-ε) * d * |S|。这意味着小集合的邻居几乎最大化地扩展。

方法与结论:如果将测量矩阵 Φ 设为该二分扩展图的邻接矩阵(并进行适当的归一化),那么它满足一种称为 RIP1 的性质(即它近似保持稀疏向量的 L1 范数)。可以证明,使用这样的矩阵,并通过特定的迭代算法或甚至基追踪,可以从测量值中恢复信号,并得到 L1 范数下的误差保证:

||x - x~||_1 ≤ C * ||x - x_k||_1

其中 C 是一个常数。这与之前基于RIP2的 L2/L1 误差保证不同。

扩展图的直观:扩展图可以看作是“没有稀疏割”的图。任何一个小节点集合,其连向集合外部的边(“边界”)数量都相对于集合内部的边数量很大。这保证了信息在图中能快速混合,也是其能用于压缩感知的深层原因。

总结

本节课我们一起学习了压缩感知中更快的解码算法。

  1. 我们介绍了迭代硬阈值法,分析了其在RIP矩阵下的性能,并概述了其误差呈指数衰减的证明思路。
  2. 我们简要了解了另一种基于扩展图的压缩感知框架,它使用满足RIP1性质的稀疏二值矩阵,并能提供 L1 范数下的恢复保证。

这两种方法提供了在测量数和解码时间之间进行权衡的不同工具,扩展了压缩感知在实际应用中的潜力。

021:高斯过程的上确界、Dudley不等式、通用链与Gordon定理

概述

在本节课中,我们将学习Gordon的“穿越网格”定理及其在降维(如Johnson-Lindenstrauss引理)中的应用。我们将探讨如何通过高斯过程的均值宽度来统一理解多种降维场景,并介绍Dudley不等式和通用链等关键工具。

高斯过程均值宽度与Gordon定理

上一节我们回顾了课程中已涉及的几种欧几里得降维方法。本节中,我们来看看一个能统一这些方法的通用定理。

Gordon定理(1988年)指出:设 Π 是一个 M × N 的随机矩阵,其元素独立同分布于 N(0, 1/M)。那么,对于单位球面 S^(N-1) 的任意子集 T,以下不等式在期望意义上成立:

sup_{x ∈ T} | ||Πx||² - 1 | ≤ ε

只要矩阵 Π 的行数 M 满足:

M ≥ C * (w(T)² / ε²)

其中,C 是一个常数,w(T) 是集合 T高斯均值宽度

高斯均值宽度 w(T) 的定义如下:设 G 是一个 N 维标准高斯随机向量(各分量独立,服从 N(0,1)),则:

w(T) = E_G [ sup_{x ∈ T} G·x ]

这个定理比之前提到的所有降维结果都更强。它可以直接推出Johnson-Lindenstrauss引理、子空间嵌入所需的最优行数,以及限制等距性质矩阵的构造。

如何估计高斯均值宽度

为了应用Gordon定理,我们需要理解如何估计给定集合 T 的高斯均值宽度 w(T)。以下是几种主要方法。

方法一:简单联合界

最直接的方法是使用联合界。由于对于每个固定的 x ∈ TG·x 是一个方差为 ||x||² = 1 的高斯随机变量。因此,sup_{x ∈ T} G·x 的期望最多是 √(log |T|)。具体推导如下:

w(T) = E[ sup_{x ∈ T} G·x ] ≤ λ + |T| * e^{-Θ(λ²)}

通过选择 λ ≈ √(log |T|),我们得到 w(T) = O(√(log |T|))

将这个界代入Gordon定理,我们得到 M ≥ C * (log |T| / ε²),这正是Johnson-Lindenstrauss引理的形式。因此,仅用联合界就足以从Gordon定理推出JL引理。

方法二:利用ε-网

当集合 T 中的点彼此非常接近时,联合界是浪费的。我们可以利用ε-网来获得更紧的界。

T'TL2 范数下的一个 ε-网。那么,我们可以将 w(T) 分解:

w(T) ≤ w(T') + ε * E[ ||G|| ]

由于 E[ ||G|| ] = O(√N),且 w(T') = O(√(log |T'|)),我们得到:

w(T) ≤ O( √(log N(T, L2, ε)) + ε√N )

其中 N(T, L2, ε)Tε-覆盖数。我们可以选择 ε 来优化这个上界。

方法三:Dudley不等式

我们可以将ε-网的思想进一步推广为一系列不同精度的网,从而得到Dudley积分不等式:

w(T) ≤ C * ∫_0^∞ √(log N(T, L2, u)) du

Dudley不等式通常已经能给出相当好的结果。例如,对于一个 d 维子空间 E 与单位球的交集 T,其覆盖数约为 (C/u)^d。代入Dudley积分可得 w(T) = O(√d),从而从Gordon定理推出子空间嵌入所需行数为 O(d/ε²)

然而,Dudley不等式有时不是最优的。例如,对于 L1 单位球,Dudley不等式给出的界是 O(√(n) log n),而实际的高斯均值宽度是 Θ(√(log n)),存在一个 log n 因子的差距。

方法四:通用链与Majorizing Measures定理

为了获得最优的、适用于任何集合 T 的界,我们需要更复杂的工具——通用链和Majorizing Measures定理。

首先定义容许序列:一个序列 {T_r}_{r≥0} 称为容许的,如果满足:

  1. T_0 ⊂ T_1 ⊂ ... ⊂ T
  2. |T_0| = 1
  3. |T_r| ≤ 2{2r}

定义 γ₂ 泛函:

γ₂(T, L2) = inf_{ {T_r} 容许 } sup_{x ∈ T} ∑_{r≥1} 2^{r/2} * dist_{L2}(x, T_r)

Talagrand的Majorizing Measures定理指出,对于任何 T,其高斯均值宽度满足:

c * γ₂(T, L2) ≤ w(T) ≤ C * γ₂(T, L2)

即,高斯均值宽度与 γ₂ 泛函在常数倍内等价。这个定理总是给出最优的界。

Gordon定理的证明思路:KMR引理

Gordon定理的完整证明涉及较多内容。这里我们概述一个由Krahmer、Mendelson和Rauhut在2014年提出的证明思路,该证明适用于次高斯分布(包括伯努利±1分布)。

KMR引理陈述如下:设矩阵 A 的条目为独立同分布的 ±1 随机变量。那么有:

E_σ [ sup_{A ∈ 𝒜} | ||Aσ||² - E[||Aσ||²] | ] ≤ C * ( γ₂(𝒜, op)² + γ₂(𝒜, op) * rad_F(𝒜) + rad_F(𝒜) * rad_op(𝒜) )

其中:

  • γ₂(𝒜, op) 是关于算子范数的 γ₂ 泛函。
  • rad_F(𝒜) 是集合 𝒜 的Frobenius范数半径。
  • rad_op(𝒜) 是集合 𝒜 的算子范数半径。

如何用KMR引理证明Gordon定理(对于±1条目)?

关键观察是,对于向量 x,我们可以构造一个特定的矩阵 A_x,使得随机投影的范数平方可以表示为:

||Πx||² = ||A_x σ||²

其中 σ 是将 Π 的所有条目连接起来的长向量。令 𝒜 = {A_x : x ∈ T},那么Gordon定理中要控制的量就等同于KMR引理左边的形式。

通过计算可以发现:

  • γ₂(𝒜, op)γ₂(T, L2) 成正比。
  • rad_F(𝒜)rad_op(𝒜) 都是 O(1)

将这些关系代入KMR引理的界,并利用Majorizing Measures定理(γ₂(T, L2) ≈ w(T)),我们最终得到:

M ≥ C * (w(T)² / ε²)

这正是Gordon定理的结论。

KMR引理本身的证明使用了鞅差序列解耦技巧Kinchin不等式以及对容许序列的复杂处理。核心思想是将目标表达式通过一个 telescoping sum(伸缩和)进行分解,然后逐项控制其期望。

总结

本节课我们一起学习了Gordon定理及其在降维中的强大应用。我们了解到,高斯均值宽度 w(T) 是衡量集合复杂度的关键几何量。为了估计它,我们介绍了从简单的联合界,到ε-网和Dudley积分不等式,再到最优的通用链与Majorizing Measures定理的一系列工具。最后,我们概述了通过KMR引理证明Gordon定理的思路,该证明揭示了随机投影保距性背后深刻的概率与几何原理。

022:KMR定理证明收尾与BPTree算法

在本节课中,我们将完成KMR定理的证明,并学习如何利用Dudley不等式等工具,为L2重击者问题设计一个内存效率更高的算法——BPTree。

KMR定理证明收尾

上一节我们介绍了KMR定理的证明框架,并定义了相关符号。本节中,我们来看看如何完成证明的最后部分。

我们定义了随机变量集合的期望上确界E,并已将其分解为几项之和。我们之前已经处理了第一项,现在需要处理中间项和右侧项。

以下是证明的核心步骤:

  1. 我们定义了一个随机变量Z,它是所有a和所有r的|X_ra|之和的上确界。
  2. 利用非负随机变量期望是其尾部概率的积分这一性质,我们通过变量替换引入了γ₂。
  3. 通过应用Khinchine不等式并对所有r和所有T_r中的元素进行联合界,我们得到了一个收敛的和式。
  4. 通过柯西-施瓦茨不等式和三角不等式,我们将E的界与E的平方根联系起来,最终得到一个关于√E的二次不等式。
  5. 解这个二次不等式,取其较大的根,平方后即可得到E的最终上界:E = O(γ₂² + γ₂ * d_F + d_F * d_op)

这个证明技巧被称为“平方根技巧”,在概率论中应用广泛。

随机过程上确界的其他应用

证明了KMR定理后,我们来看看随机过程上确界分析的其他应用。

以下是两个重要的应用方向:

  1. 证明SH满足RIP:其中S是采样矩阵,H是哈达玛矩阵。这可以用于分析快速JL变换。
  2. 为插入流中的L2重击者问题设计更低内存的算法:例如BPTree算法,它使用O(k log k)字的内存,而Count Sketch需要O(k log n)字。目前尚不清楚O(k)是否可能。

接下来,我们将重点介绍BPTree算法。

BPTree算法:归约到超重问题

BPTree算法的核心思想是将L2重击者问题归约到一个更简单的问题——超重问题。

超重问题定义:一个项i被称为超重的,如果其频率的平方x_i²大于其他所有项频率平方和的1000倍(或某个大常数倍)。

归约过程如下:

  1. 我们创建q = O(k)个独立的超重问题求解器副本D_1, ..., D_q
  2. 使用一个哈希函数h将全域[n]映射到{1,...,q}
  3. 当在流中看到项i时,将其发送到对应的求解器D_h(i)进行处理。
  4. 对于任意一个重击者i,在它被哈希到的桶中,由于其他重击者碰撞进来的期望数量很少(k/q),且来自非重击者的噪声平方和期望也很小,因此i在其桶中极有可能是一个超重项。
  5. 为了以高概率找出所有k个重击者,我们将上述结构重复r = Θ(log k)次。最终,如果一个项在超过r/2个副本中被报告,则将其输出。

这个归约过程保证了输出列表大小为O(k),且不会漏掉任何真正的重击者。总内存开销为O(k log k * S),其中S是解决单个超重问题所需的内存。

BPTree算法:以常数内存解决超重问题

现在,我们来看如何用常数内存S解决超重问题。这是算法的核心创新。

假设存在一个超重项h。我们的目标是逐位学习其二进制表示h_0, h_1, ..., h_{log n}

核心引理:设y_0, y_1, ..., y_T是一个仅插入流的频率向量演化过程。设σ_1,..., σ_n是四向独立的±1随机变量。那么,期望E[sup_t |σ · y_t|] = O(||y_T||_2)。这个引理是Levy极大不等式在四向独立和任意插入流上的推广。

利用引理逐位学习

  1. 为了学习第0位,我们初始化两个计数器B_0B_1,并生成四向独立的随机符号σ
  2. 当流中出现项i时,检查其最低位。若为0,则向B_0σ_i;若为1,则向B_1σ_i
  3. 设Δ是超重项h的最终频率x_h。我们等待,直到某个计数器的绝对值超过Δ/10
  4. 根据引理,代表错误项的计数器(如B_0)的噪声始终很小(O(||y_T||_2)),而||y_T||_2远小于Δ。
  5. 同时,当h出现足够多次(约Δ/10次)后,代表正确项的计数器(如B_1)的值将接近σ_h * (Δ/10),其幅度会超过阈值。
  6. 因此,首先超过阈值的计数器所对应的位,极有可能就是h的第0位。

学习所有位

  1. 学习第0位会消耗h的一小部分出现次数(如10%)。
  2. 在学习下一位时,我们只处理那些前缀与已学位匹配的项。由于使用了随机排列,剩余噪声项的数量会以因子2的指数级减少。
  3. 因此,h相对于剩余噪声的“超重”程度指数级增加,学习后续位所需消耗的h出现次数指数级减少。
  4. 总消耗次数是一个收敛的几何级数,可以在流结束前学完h的所有位。

处理未知的Δ

  1. 我们并行运行多个算法副本,每个副本假设不同的Δ猜测值(2的幂次)。
  2. 同时,运行一个AMS草图来持续估计当前流的L2范数。
  3. 当AMS草图指示当前L2范数远超某个猜测值时,就丢弃对应的副本,并启动一个具有更高猜测值的新副本。
  4. 利用引理的类似版本,可以证明AMS草图在所有时间步上的最大误差有界,因此无需为联合界付出额外内存开销。

引理的证明思路

  1. 将问题转化为求一个Rademacher过程的上确界期望。
  2. 应用Dudley熵积分,其上界依赖于集合在L2范数下的覆盖数。
  3. 对于插入流产生的向量序列,可以证明其ε-覆盖数至多为O(1/ε²)
  4. 将此覆盖数界代入Dudley积分,得到一个收敛的级数,从而证明期望上确界为O(||y_T||_2)
  5. 关键点在于,四向独立性足以保证Rademacher过程的尾部概率衰减足够快,使得Dudley积分收敛。

本节课中,我们一起学习了如何完成KMR定理的证明,并深入探讨了BPTree算法。该算法通过将L2重击者问题归约到超重问题,并利用关于随机过程上确界的新颖引理,实现了O(k log k)字的内存消耗,是数据流算法设计中的一个优美范例。

023:通过采样DFT实现RIP,稀疏傅里叶变换

在本节课中,我们将学习稀疏傅里叶变换。核心思想是,当我们知道一个信号的傅里叶变换是稀疏或近似稀疏时,能否设计出比标准快速傅里叶变换更快的算法来近似计算它。

概述

我们被给予对一个n维信号 x 的查询访问权限。我们希望计算其离散傅里叶变换 。标准FFT算法可以在 O(n log n) 时间内完成此计算。然而,在许多应用中,例如图像和音频压缩,信号的傅里叶变换是近似稀疏的。如果 是近似k稀疏的,我们能否在接近 k 的时间(例如 k polylog(n))内近似计算它?本节课将探讨这个可能性。

首先,我们将展示,即使只查看信号 x 的极少部分条目,理论上也包含足够的信息来恢复稀疏的傅里叶变换。我们将通过证明对逆傅里叶变换矩阵进行随机行采样能满足限制等距性质来实现这一点。这为从次线性数量的时域样本中进行压缩感知恢复提供了概念证明。

然而,基于RIP的恢复算法(如基追踪)在计算上并不高效。因此,我们将转向一种专门针对稀疏傅里叶变换设计的、真正次线性时间的算法。我们将重点讨论 k=1 的情况,展示如何通过巧妙的采样和比较,以 O(log n log log n) 的查询复杂度和运行时间,鲁棒地恢复出主导的傅里叶系数及其位置。

从概念到可行性:采样与RIP

上一节我们提出了能否通过次线性查询近似稀疏傅里叶变换的问题。本节中,我们首先来证明这在信息论上是可行的:通过随机采样时域信号,我们获得的测量值包含足够的信息。

我们考虑逆离散傅里叶变换矩阵 F⁻¹。定义采样矩阵 S 为一个 n×n 的对角矩阵,其对角线元素 δ_j 是独立的伯努利随机变量,取1的概率为 m/n,取0的概率为 1 - m/n。这相当于从时域信号 x 中随机选取 m 个样本。

核心主张:缩放后的矩阵 (1/√m) S F⁻¹ 以高概率满足 k 阶限制等距性质,只要 m = O(k polylog(n))

这意味着,给定测量值 y = (1/√m) S x(即随机时域样本),我们可以通过解决一个优化问题(如基追踪)来近似恢复 k 稀疏的 ,并保证如下形式的误差界:
‖x̂ - x̃‖₂ ≤ C * ‖x̂_{tail(k)}‖₂
其中 x̂_{tail(k)} 中除最大k个分量外的部分。

证明思路简述
证明的关键在于界定以下期望值:
E_δ [ sup_{|T|=k} ‖ I_T - (1/m) Σ_{j} δ_j x_j_T x_j_T^ ‖ ]*
其中 x_jF⁻¹ 的第j行。通过对称化、引入Rademacher随机变量 σ,并利用算子范数的定义,我们可以将问题转化为界定一个Rademacher过程的期望上确界。

应用 Dudley 积分不等式,该上确界可由覆盖数的积分来界定。这里的距离度量是一种依赖于采样集 Ω(即 δ_j=1 的索引j的集合)的“星范数”。通过结合体积论证(对小半径)和 Maurey 经验方法(对大半径)来界定覆盖数,并利用一个关键的平方根技巧来处理自引用项,最终可以证明,当 m = O(k polylog(n)) 时,RIP 常数可以足够小。

这个结果非常重要,因为它从理论上证实了从 O(k polylog(n)) 个时域样本中恢复稀疏频谱是可能的。然而,正如之前指出的,基于RIP的恢复算法通常需要至少 O(n) 的计算时间,这并不比FFT快。因此,我们需要寻找真正次线性的算法。

精确单稀疏情况下的次线性算法

在深入噪声鲁棒性算法之前,我们先看一个最简单的情况: 是严格单稀疏的,即只有一个非零傅里叶系数。这展示了次线性时间算法的基本可能性。

假设 x̂_u ≠ 0,且对于所有 v ≠ u,有 x̂_v = 0。根据逆DFT公式,时域信号为:
x_j = (1/n) Σ_{a=0}^{n-1} x̂_a ω^{-a j} = (1/n) x̂_u ω^{-u j}
其中 ω = e^{-2πi / n}

算法

  1. 查询 x_0。根据公式,x_0 = (1/n) x̂_u
  2. 查询 x_1。根据公式,x_1 = (1/n) x̂_u ω^{-u}
  3. 计算 x̂_u = n * x_0
  4. 计算 ω^{-u} = x_1 / x_0。由于 ω 是本原单位根,我们可以从中解出频率索引 u

因此,仅通过查看两个时域条目,我们就可以精确计算出单稀疏的傅里叶变换。然而,这个方案对噪声非常敏感。如果 只是近似稀疏(即存在小能量的尾部),该算法将失效。

鲁棒的次线性时间算法(k=1)

上一节我们看到了无噪声下的简单方案。本节中,我们来看看如何设计一个能够容忍近似稀疏性(即存在噪声尾部)的算法。我们将专注于 k=1 的情况,目标是找到一个估计 ,使得:
‖x̃ - x̂‖₂ ≤ O( ‖x̂_{tail(1)}‖₂ )
其中 x̂_{tail(1)} 是去除最大分量后的向量。

算法的核心思想是逐比特确定主导频率 u 的二进制表示。假设 u 是一个 log₂ n 比特的数字。

基本无噪声情况下的比特提取
首先,我们如何获取 u 的最低有效位(奇偶性)?

  1. 计算 A = x_0 + x_{n/2}B = x_0 - x_{n/2}
  2. 根据DFT性质,可以推导出:
    • 如果 u 是偶数,则 x_{n/2} = (1/n) x̂_u,因此 A ≈ 2*(1/n)x̂_uB ≈ 0
    • 如果 u 是奇数,则 x_{n/2} = -(1/n) x̂_u,因此 A ≈ 0B ≈ 2*(1/n)x̂_u
  3. 比较 |A||B|。较大的那个指示了 u 的奇偶性。

在得知最低有效位后,我们可以通过一个简单的技巧(在时域乘以一个相位因子 ω^j,相当于在频域进行循环移位),“归零”已确定的低位,使得 u 在新的等效问题中变成2的倍数。然后,我们可以通过比较 x_0 + x_{n/4}x_0 - x_{n/4} 来提取下一位比特。以此类推,重复 log₂ n 次即可确定所有比特。

引入噪声与随机化
当存在噪声时,问题在于噪声可能会集中在某次比较所使用的特定时域点(如 x_0x_{n/2})上,从而导致比特判断错误。

解决方案是随机化采样位置,以平均掉噪声的影响。
为了判断第 i 个比特:

  1. 不再固定使用 x_0x_{n/2^i},而是随机均匀选取一个偏移量 r
  2. 计算 A_r = x_r + x_{r + n/2^i}B_r = x_r - x_{r + n/2^i}(索引模 n)。
  3. 比较 |A_r||B_r|
  4. 独立重复此过程 O(log log n) 次,取多数结果作为该比特的估计。

原理:由于逆傅里叶变换是正交的,频域噪声在时域是均匀扩散的(具有相同的总能量)。随机选择 r 使得在 x_rx_{r+n/2^i} 处遇到大噪声的概率很小。通过多次重复和多数表决,可以以高概率正确判断每个比特。

算法复杂度

  • 查询复杂度:需要判断 log₂ n 个比特,每个比特重复 O(log log n) 次,每次查询2个点。总查询次数为 O(log n log log n)
  • 时间复杂度和:与查询复杂度同阶,为 O(log n log log n)
  • 恢复值:在确定频率索引 u 后,可以通过对少数几个 x_j 采样并解一个小线性方程组来估计 x̂_u 的值。

这种方法实现了真正的次线性时间,并且对近似稀疏性具有鲁棒性。

总结

本节课中我们一起学习了稀疏傅里叶变换的基本概念和算法思路。

我们首先探讨了从次线性时域样本中恢复稀疏频谱的理论可能性,通过证明随机采样逆DFT矩阵满足RIP性质,确立了信息论上的可行性。这为压缩感知在此场景的应用奠定了基础。

接着,我们针对最简单的单稀疏精确情况,展示了一个仅需两次查询的极简算法,揭示了次线性时间算法的潜力。

最后,也是最重要的,我们详细介绍了一种鲁棒的、次线性时间的算法,用于处理近似单稀疏的情况。该算法通过随机化采样位置和逐比特频率恢复的策略,以 O(log n log log n) 的查询和计算时间,能够高概率地恢复出主导频率及其系数,并满足 L₂ 误差界。这为在诸如压缩等实际应用中快速计算近似稀疏的傅里叶变换提供了有希望的途径。

对于更大的 k,算法思想可以推广,但需要更复杂的技术来识别和分离多个频率分量,这超出了本节课的范围。

谢谢,教授。

024:核心集 (Coresets)

在本节课中,我们将要学习核心集 (Coresets) 这一概念。核心集是一种数据压缩技术,它通过保留一个加权的数据子集,使得在该子集上进行优化计算的结果,能够近似等于在整个原始数据集上进行优化的结果。这为设计低内存的流式算法提供了一种通用框架。


核心集的定义

核心集分为弱核心集和强核心集,它们的共同点是:对于一个给定的优化问题,核心集是一个加权的输入数据子集。

  • 弱核心集:在该加权子集上找到的(近似)最优解,对于原始数据集也是(近似)最优的。
  • 强核心集:对于任意一个解,该解在加权子集上的目标函数值,与在原始数据集上的目标函数值近似相等。

我们可以通过一个 K-means 聚类的例子来直观理解。假设我们有一些数据点(黑色),最优的聚类中心是红色点。一个核心集(黄色点)是一个加权子集。弱核心集意味着,如果我们只在黄色点上运行 K-means 算法并找到聚类中心(例如蓝色点),那么这些中心对于原始的黑色点集也是近似最优的。


核心集与流式算法

核心集的构造可以系统地转化为低内存的流式算法。其基本思想是使用一种合并方法。

假设我们有一个流式数据点序列。我们可以构建一棵概念上的完美二叉树:

  1. 将到达的数据点暂存为微小的核心集(例如单个点本身)。
  2. 当暂存的数据达到一定数量时,我们运行核心集构造算法,将它们合并为一个新的、更小的核心集。
  3. 重复此过程,像合并二叉树节点一样,逐层向上构建核心集。

在这个过程中,我们最多同时存储 O(log n) 个核心集。如果原始核心集的大小为 S,那么流式算法的总空间复杂度就是 O(S log n)。核心集的近似质量 β 在合并过程中可能会累积,但通过精心设计(例如使 β ≈ 1 + ε/log n),最终可以控制总体误差。


示例一:K-Median 核心集构造

本节我们来看一个针对 K-Median 问题的核心集构造方法。K-Median 的目标是找到 k 个中心点,最小化所有数据点到其最近中心点的距离之和。

公式

给定数据集 X,最小化:Σ_{p ∈ X} d(p, C)
其中 C 是大小不超过 k 的中心点集合,d(p, C) = min_{c ∈ C} d(p, c)

我们假设存在一个离线的 B 近似算法(B 为常数),它可以在线性空间内运行。以下是构造方法:

  1. n 个数据点分成 n/B 块,每块大小为 B
  2. 在每个数据块上独立运行 B 近似算法,得到该块的 k 个近似中心点。
  3. 将所有块的中心点收集起来,形成核心集。

核心集的大小为 (n/B) * k。我们需要在内存中同时存储每个块的数据(O(B))和所有中心点(O((n/B)*k))。通过设置 B = √(nk) 来平衡,我们得到总空间和核心集大小为 O(√(nk))

可以证明,对此核心集使用 B 近似算法得到的解,对于原始数据集的近似比约为 4B(B+1)。虽然这个核心集大小仍依赖于 √n,不如已知的最佳结果(例如 poly(k, d, 1/ε)),但它是一个在任意度量空间中都有效的非平凡构造。


示例二:最小包围球核心集构造

本节我们看看 最小包围球 问题的核心集构造。该问题是在欧几里得空间中,找到能包围所有点的最小半径的球。

公式

给定点集 X,最小化:max_{p ∈ X} || p - o ||_2
其中 o 是球的中心。

我们利用 网 (Net) 的概念来构造核心集。具体步骤如下:

  1. 在单位球面上构造一个 α-网 V = {v1, v2, ..., vT}。这意味着对于球面上任何方向 x,网中都有一个点 vi,使得它们之间的夹角不超过 α。这样一个网的大小为 T = O(1/α)^{d-1},其中 d 是空间维度。
  2. 对于网中的每一个方向 vi,我们存储原始点集中在该方向上投影最大的点,即 argmax_{p ∈ X} (p · vi)
  3. 这些存储的点就构成了我们的核心集,其大小为 T

分析
设核心集的最小包围球半径为 R_c,其中心为 O(可通过平移假设为原点)。对于任何原始点 p,设其方向为 w。在 α-网中存在一个方向 v',使得 ∠(w, v') ≤ α。由于 p 不在核心集中,那么核心集中必然存在另一个点 p' 满足 p' · v' ≥ p · v'。因为 p' 在半径为 R_c 的球内,所以 p' · v' ≤ R_c。通过三角关系和点积计算,可以推导出 ||p|| ≤ R_c / cos(α)。当 α 很小时,cos(α) ≈ 1 - Θ(α^2)。为了使 ||p|| ≤ R_c * (1+ε),我们需要 α^2 = Θ(ε),即 α = Θ(√ε)。因此,核心集的大小为 O(1/ε^{(d-1)/2})

这个构造的优点是核心集大小与原始数据点数量 n 无关,只与维度 d 和精度 ε 有关。


总结

本节课我们一起学习了核心集的概念及其在流式算法中的应用。我们了解到:

  1. 核心集是一种数据摘要技术,能在压缩数据的同时保留解决优化问题所需的关键信息。
  2. 弱核心集保证了解的质量,而强核心集进一步保证了目标函数值的近似保持。
  3. 通过树形合并方法,可以将核心集构造转化为高效的流式算法。
  4. 我们分析了 K-Median最小包围球 两个问题的核心集构造实例,看到了不同方法在核心集大小和近似质量上的权衡。

核心集是处理大规模数据,特别是在内存受限的流式场景下,一个非常有力的工具。

posted @ 2026-03-29 09:26  布客飞龙II  阅读(2)  评论(0)    收藏  举报