哈佛-CS229r-大数据算法笔记-全-

哈佛 CS229r 大数据算法笔记（全）

1：课程介绍与近似计数算法 🚀

在本节课中，我们将学习大数据算法课程的基本框架，并深入探讨第一个核心算法——莫里斯近似计数算法。我们将从课程概述开始，了解大数据算法所面临的挑战和主要研究方向，然后详细分析一个经典的流式算法，该算法能以极小的空间代价近似地统计事件数量。

课程概述与物流信息 📋

欢迎来到CS 229R《大数据算法》课程。课程网站包含了所有必要信息。请通过课程指定的邮箱地址联系我和助教Yak Yaroswab，以确保邮件能被正确分类和处理。

课程讲座将被录制，视频链接会发布在课程网站上。课程要求包括完成习题集、参与一个最终项目并进行演示。此外，每位学生都需要至少负责一次讲座笔记的记录工作。

什么是大数据算法？ 🤔

本课程关注处理大规模数据的算法。所谓“大数据”，指的是数据规模如此之大，以至于传统的算法分析方法（例如在本科算法课程中学到的）变得不再适用。

例如，数据可能大到无法完全装入内存，因此算法只能对数据进行单次扫描；或者数据虽然能存储在计算机上，但无法放入内存，必须使用磁盘，而磁盘访问速度远慢于内存。这些都是大数据算法需要解决的问题。

课程主题概览 📚

以下是本课程将涵盖的几个核心主题：

1. 素描与流式计算

素描：数据 X 的素描 C(X) 是关于某个函数 F 的一种数据压缩。其目标是仅通过访问压缩后的素描 C(X)（有时是 C(X) 和 C(Y)），就能计算或近似原始数据 X 的函数值 F(X)（或 F(X, Y)）。
流式计算：指在数据 X 动态更新（例如不断有数据到达）时，在线维护其素描 C(X) 的算法。

2. 降维
数据通常是高维向量。降维技术能够将数据映射到低维空间，使得在低维数据上解决计算问题后，可以将解“提升”回原始高维空间。由于算法运行时间或所需资源通常依赖于输入维度，降维可以显著加速算法，例如聚类或高维最近邻搜索。

3. 大规模机器学习
我们将从算法角度探讨一些大规模机器学习问题，例如如何加速求解过程。一个典型例子是大规模线性回归。

线性回归公式：
假设我们有观测数据 (Z_i, b_i)，满足 b_i ≈ F(Z_i) + 噪声。在线性回归中，F 是线性函数，即 F(Z) = X · Z。目标是找到参数向量 X。最小二乘回归通过最小化残差平方和来求解：
X_ls = argmin_X ||Z * X - b||_2
其中，矩阵 Z 的行是 Z_i，向量 b 的元素是 b_i。其闭式解为（当 Z 列满秩时）：
X_ls = (Z^T Z)^{-1} Z^T b
当数据量极大时，直接计算这个解（涉及大矩阵乘法）非常昂贵。本课程将介绍能更快获得良好近似解的技术。

4. 压缩感知
其核心思想是利用高维信号（如图像）在某个基（如小波基）下的稀疏性，通过少量线性测量来廉价地获取或压缩信号。这使得从少量测量中重建信号成为可能，应用包括JPEG图像压缩和加速MRI扫描。

5. 外部存储模型
在传统算法分析中，我们衡量算法效率的标准是执行指令的数量。而在外部存储模型中，我们主要衡量磁盘I/O次数。这是因为内存访问和磁盘随机寻道之间存在巨大的速度差异（可达六个数量级）。

该模型将磁盘视为由大小为 B 的块组成的无限磁带，内存则能容纳 M 个数据项，分为 M/B 个页。当要访问的数据不在内存中时，需要从磁盘读取相应块到内存（可能还需将内存中的某页写回磁盘），这计为一次I/O操作。算法的目标是最小化I/O次数。B树就是为该模型设计的数据结构。

此外，时间允许的话，我们可能还会探讨其他模型，如MapReduce。

第一个流式算法：近似计数 🧮

现在，我们开始探讨第一个流式算法：近似计数问题。该问题由莫里斯在20世纪70年代提出，被认为是流式算法的开创性工作之一。

问题定义：
有一系列事件发生。每当一个事件到达，我们需要更新计数。目标是设计一个算法，在任意时刻当被询问“发生了多少事件”时，能输出一个估计值 N_tilde，满足：
P(|N_tilde - N| > εN) ≤ δ
其中，N 是真实事件数，ε 和 δ 是预先给定的误差参数和失败概率。算法的核心目标是最小化空间使用。

精确计数的空间下限：
对于精确计数，只需维护一个计数器，每次事件到达时加1。如果事件数 N 的上限已知，存储计数器需要 O(log N) 比特。这是最优的，否则由鸽巢原理可知必然出错。

近似计数的空间下限：
对于近似计数（例如，答案允许有常数倍误差），直觉上我们只需要知道 N 落在哪个2的幂次区间内。区间数量为 O(log N)，因此区分它们只需要 O(log log N) 比特。这给出了一个空间下界。

莫里斯算法

莫里斯算法达到了 O(log log N) 的空间复杂度。

算法步骤：

初始化变量 X = 0。
对每个到达的事件，以概率 1 / (2^X) 将 X 增加1。
当需要输出估计值时，返回 N_tilde = 2^X - 1。

算法分析：
我们设 X_N 为处理完 N 个事件后 X 的值。

引理1：期望的正确性
E[2^{X_N}] = N + 1
证明（归纳法）：

基础情况：N=0 时，X=0，2^0 = 1 = 0+1，成立。
归纳步骤：假设对 N 成立，即 E[2^{X_N}] = N+1。考虑第 N+1 个事件。当 X_N = j 时，X_{N+1} 有两种可能：
- 以概率 1 - 1/(2^j) 保持为 j，此时 2^{X_{N+1}} = 2^j。
- 以概率 1/(2^j) 增加为 j+1，此时 2^{X_{N+1}} = 2^{j+1}。
  因此，条件期望 E[2^{X_{N+1}} | X_N = j] = 2^j * (1 - 1/(2^j)) + 2^{j+1} * (1/(2^j)) = 2^j + 1。
  对 j 求期望，E[2^{X_{N+1}}] = E[2^{X_N}] + 1 = (N+1) + 1 = N+2。归纳成立。

引理2：二阶矩分析
通过类似的归纳法可以证明（详细步骤见课程笔记）：
E[2^{2 X_N}] = (3/2) N^2 + (3/2) N + 1
由此可以推导出估计值的方差：
E[(N_tilde - N)^2] = E[(2^{X_N} - 1 - N)^2] ≤ (1/2) N^2

应用切比雪夫不等式：
根据切比雪夫不等式：
P(|N_tilde - N| > εN) ≤ Var(N_tilde) / (ε^2 N^2) ≤ (1/2) / ε^2
这意味着失败概率有一个常数上界（例如，当 ε=1 时，失败概率不超过 1/2）。

改进算法：Morris+ 和 Morris++

为了得到任意小的误差 ε 和失败概率 δ，我们可以对基础算法进行增强。

Morris+ 算法：

并行运行 s 个独立的莫里斯算法实例。
输出这 s 个实例估计值的平均值。
由于方差可加，求平均后估计值的方差减小为原来的 1/s。因此，失败概率变为 O(1/(s ε^2))。要使得失败概率小于 δ，需要设置 s = O(1/(ε^2 δ))。空间复杂度为 O(s * log log N)。

Morris++ 算法：

运行 t 个独立的 Morris+ 实例，每个实例配置其失败概率为常数（例如 1/3）。
输出这 t 个 Morris+ 实例估计值的中位数。
中位数估计器失败，意味着超过一半的 Morris+ 实例失败。每个 Morris+ 实例失败概率为常数（如 1/3）。利用切尔诺夫不等式可以证明，当 t = O(log(1/δ)) 时，整体失败概率可降至 δ。因此，总空间复杂度为 O( (log(1/δ))/ε^2 * log log N )。

这种“先取平均，再取中位数”的方法是一种通用技巧，可以将一个较弱的估计器增强为高概率成功的强估计器。

总结 ✨

本节课中，我们一起学习了大数据算法课程的总体框架，明确了处理海量数据时在内存、I/O等方面面临的独特挑战。我们深入探讨了第一个流式算法——莫里斯近似计数算法。该算法巧妙地利用随机性，仅用 O(log log N) 比特的空间就能对事件数量进行近似统计。通过分析其期望和方差，并运用切比雪夫不等式，我们理解了其理论保证。进一步，我们看到了如何通过“并行求平均”和“取中位数”的增强技巧（Morris+, Morris++），将算法推广到满足任意精度 ε 和置信度 δ 的要求。这为我们后续学习更复杂的素描与流式算法奠定了坚实的基础。

2：流式算法与不同元素问题 🧮

在本节课中，我们将继续学习流式算法，并重点探讨一个经典问题——不同元素问题。我们将从理想化的算法开始，逐步将其转化为实际可实现的版本，并分析其空间复杂度。

回顾与引言

上一讲我们介绍了近似计数算法，即莫里斯算法。该算法提供了一个无偏估计器，并通过取多个估计器的平均值（Morris+）以及这些平均值的中位数（Morris++）来获得高概率下的良好近似。这种“先平均，再取中位数”的思想在流式算法设计中非常常见。

本节课我们将应用类似的思想，来研究另一个极端问题：不同元素问题。

不同元素问题（F0问题）🔍

问题定义

我们观察一个由整数组成的流，每个更新是一个介于1到n之间的数字。算法的目标是输出流中不同整数的数量。

我们用 T 表示真实的答案。与之前一样，我们的目标是最小化空间消耗。

两种简单算法

以下是两种直接的解决方案：

位数组法：维护一个长度为 n 的位数组。当看到整数 i 时，将第 i 位设为1。最后，统计数组中1的个数。空间复杂度为 O(n) 比特。
存储列表法：直接存储所有出现过的数字。如果流长度为 M，每个数字需要 log n 比特存储，总空间为 O(M log n) 比特。

因此，精确解的空间复杂度至少为 min(n, M log n)。然而，我们将允许算法输出一个随机化的近似解。

近似目标

我们希望算法输出一个估计值 T̃，使得以下不等式以高概率成立：
Pr[ |T - T̃| > εT ] ≤ δ

其中 ε 是近似参数，δ 是失败概率。

理想化算法：Flajolet-Martin 算法

首先，我们介绍一个由 Flajolet 和 Martin 在1983年提出的理想化算法。该算法是后续实际算法的基础。

算法步骤

选取哈希函数：随机选取一个哈希函数 h，将输入域 [1, n] 映射到连续区间 [0, 1]。
维护最小值：在内存中仅维护一个数字 X，它是所有看到的 h(i) 的最小值。
输出估计：处理完流后，输出 (1/X) - 1。

算法原理与期望

X 实际上是 T 个独立均匀分布于 [0,1] 的随机变量的最小值。直觉上，T 个随机点均匀分布在 [0,1] 区间，第一个点的期望位置大约是 1/(T+1)。

证明期望：
E[X] = ∫₀¹ (1-λ)^T dλ = 1/(T+1)

因此，1/X 的期望大约是 T+1，我们的估计器 (1/X)-1 的期望大约是 T。

方差分析

为了应用切比雪夫不等式，我们还需要知道 X 的方差。经过计算可得：
Var(X) ≤ 1/(T+1)^2

这意味着方差随着 T 增大而减小，是一个良好的性质。

从理想到现实：改进算法

上述理想化算法存在两个问题：

存储一个完全随机的哈希函数 h 需要 O(n) 空间。
算法操作的是实数，而计算机只有有限精度。

接下来，我们通过“平均”和“取中位数”的技巧来改进算法，并最终解决哈希函数的存储问题。

FM+ 算法：通过平均降低方差

我们独立运行 Q = 1/(ε²η) 个上述理想化算法实例。设第 i 个实例维护的最小值为 X_i。

FM+ 算法输出：(1 / ( (1/Q) * Σ X_i )) - 1

分析：

期望：由于线性性，平均值 (1/Q) Σ X_i 的期望仍然是 1/(T+1)。
方差：独立随机变量和的方差等于方差之和。平均值后的方差变为原来的 1/Q，即 Var_avg ≤ (1/(T+1)^2) / Q。

应用切比雪夫不等式，可以证明估计值以至少 1-η 的概率满足 (1-ε)T ≤ T̃ ≤ (1+ε)T。

FM++ 算法：通过中位数提升成功概率

FM+ 算法以常数概率（例如 η=1/3）给出好估计。为了将成功概率提升到 1-δ，我们采用取中位数的方法。

FM++ 算法步骤：

独立实例化 S = O(log(1/δ)) 个 FM+ 算法（每个设置 η=1/3）。
输出这 S 个估计值的中位数。

分析：
定义指示变量 Y_j，当第 j 个 FM+ 实例失败（估计误差过大）时为1，否则为0。E[Y_j] ≤ 1/3。
中位数失败意味着超过一半的实例失败，即 Σ Y_j > S/2。这要求 Σ Y_j 显著偏离其期望 ≤ S/3。利用切尔诺夫界可以证明，当 S = O(log(1/δ)) 时，失败概率被限制在 δ 以内。

最终空间复杂度（忽略哈希函数存储）：O( (1/ε²) * log(1/δ) )

实现真正的有限空间：哈希函数与层次化算法

现在，我们解决哈希函数的存储问题，并给出一个实际可实现的算法。

有限独立哈希函数族

我们无法存储完全随机的函数。替代方案是使用 k 阶独立哈希函数族。一个关键的例子是低阶多项式哈希：

构造：对于素数 p，所有次数小于 k 的多项式构成一个哈希函数族。
性质：该族是 k 阶独立的。存储这样一个函数只需要存储 k 个系数，即 O(k log p) 比特。当 k=2（两两独立）时，仅需 O(log p) 比特。

一个常数因子近似算法

我们首先设计一个空间为 O(log n) 比特的算法，它能以高概率给出 T 的常数因子近似。

算法思路：

选取一个两两独立的哈希函数 h，将 [1, n] 映射到 [1, n]（假设 n 是2的幂）。
维护一个变量 X，它是流中所有元素 h(i) 的最低有效非零位索引的最大值（例如，h(i)=...001000，则最低有效非零位索引是3）。
输出 2^X。

直观理解：
对于一个特定的 i，h(i) 的最低有效非零位索引为 j 的概率约为 1/2^(j+1)。因此，期望有 T/2^(j+1) 个元素的哈希值最低有效非零位索引为 j。当 j ≈ log₂ T 时，期望数量约为常数。算法记录的 X 大致就在 log₂ T 附近，因此 2^X 是 T 的常数因子近似。通过方差分析（利用两两独立性）和切尔诺夫/马尔可夫不等式，可以严格证明这一点。

提升到 (1+ε) 近似：层次化算法

为了得到更精确的 (1+ε) 近似，我们采用层次化（分桶）的思想。

算法框架：

并行运行多个“基础估计器”，每个负责一个“层次”。
每个层次 j 对应一个两两独立的哈希函数 g_j，它将元素哈希到两个桶之一（模拟最低有效位为 j 的判断）。
对于层次 j，我们运行一个简易求解器，它只精确记录前 C/ε² 个映射到该桶的不同元素。如果不同元素数量超过此阈值，求解器只输出阈值。
同时，运行前述的常数因子近似算法，以确定一个层次 j*，使得期望映射到该桶的元素数量 E ≈ 1/ε²。
最终，输出 2^(j*+1) * (层次 j* 上简易求解器的输出)。

为何有效：

在选定的层次 j*，期望有 Θ(1/ε²) 个不同元素被映射过来。由于两两独立性，实际数量将高度集中在期望值附近（切比雪夫不等式），相对误差约为 ε。
简易求解器在该层次上能精确计数（因为数量级约为 1/ε²，未超过其容量）。
常数因子近似算法确保了我们能找到这样一个合适的层次 j*。

最终空间复杂度：

常数因子近似算法：O(log n) 比特。
哈希函数族：O(log n) 比特。
简易求解器：共 O(log n) 个，每个使用 O((1/ε²) log n) 比特。
总计：O((1/ε²) log² n) 比特。

最优性：已知存在算法可以达到 O((1/ε²) + log n) 比特的空间，并且这已被证明是最优的。

总结与下节预告

本节课我们一起深入学习了不同元素问题（F0问题）的流式算法解决方案：

我们从 Flajolet-Martin 理想化算法 出发，理解了其基于最小哈希值的核心思想。
通过 FM+（平均）和 FM++（中位数）技术，我们将算法转化为具有高成功概率的实用随机近似算法。
我们探讨了 有限独立哈希函数族（特别是两两独立族）如何用少量空间模拟完全随机函数，从而解决了理想算法中哈希函数存储开销大的问题。
最后，我们介绍了一个 层次化算法，通过并行运行多个不同精度的估计器并结合一个常数因子近似器作为“指导”，最终实现了空间复杂度为 O((1/ε²) log² n) 的非理想化、可实现的 (1+ε) 近似算法。

在下一讲中，我们将探讨该问题的下界，即证明如果不允许随机化或近似，任何算法都无法在小于 min(n, M log n) 的空间内解决此问题。我们还将简要介绍另一个流式算法问题——巨大偏差问题。

3：下界与转盘模型 🧮

在本节课中，我们将要学习大数据流算法中的两个核心概念：下界证明与转盘模型。我们将首先探讨为什么对于“不同元素计数”问题，随机化和近似是必不可少的。接着，我们将转向一个更通用的流模型——转盘模型，并学习一种强大的算法设计技术：线性草图。

下界：为什么我们需要随机化和近似？ 🚫

上一节我们介绍了用于不同元素计数的Flajolet-Martin算法及其变种。该算法使用了随机化和近似。你可能会问，我们能否放松这些要求？例如，设计一个确定性的精确算法？或者一个确定性的近似算法？又或者一个随机化的精确算法？本节中，我们将通过下界证明来展示，对于“不同元素计数”问题，随机化和近似两者都是必需的。

确定性精确算法需要线性空间

核心概念：任何使用空间少于 n 比特的确定性精确算法都是不存在的。

证明思路：我们通过构造一个压缩方案来证明。假设存在一个空间为 S 比特的确定性精确算法 A。我们可以利用 A 来编码任意一个 n 比特的向量 x。

编码：根据向量 x 创建一个流。对于 x 中每一个为 1 的位 i，将元素 i 加入流中。然后运行算法 A 处理这个流，并将其最终的内存状态（S 比特）作为 x 的编码。
解码：给定编码（即 A 的内存状态），解码器可以通过“询问”算法 A 来逐位恢复 x。对于每一位 i，解码器将 i 追加到当前流中，然后查询 A 得到的不同元素计数。如果计数不变，则 x_i = 1；如果计数增加，则 x_i = 0。

由于这是一个从所有 2^n 个可能的 n 比特向量到 S 比特字符串的单射，根据鸽巢原理，必须有 S ≥ n。因此，任何确定性精确算法都需要至少 n 比特的线性空间。

确定性近似算法也需要线性空间

核心概念：即使允许近似（例如 (1 ± 0.1) 近似），任何确定性算法仍然需要 Ω(n) 比特的空间。

证明思路：我们同样通过编码论证，但需要构造一个特殊的向量集合 C。这个集合满足：

所有向量中 1 的数量（权重）相同，例如 n/100。
集合 C 的大小是指数级的（例如 2^{cn}）。
集合中任意两个不同向量的交集很小（例如最多 n/1000）。

如果我们有一个空间为 S 的确定性 (1±0.1) 近似算法 A，我们可以用它来编码集合 C 中的向量。

编码：与之前类似，用向量 x_S ∈ C 创建流，并存储 A 处理后的内存状态。
解码：解码器遍历集合 C 中的每一个候选向量 x_S‘。它将 x_S‘ 对应的流追加到当前状态，然后查询 A。
- 如果 x_S‘ 是正确的原始向量，那么追加后不同元素数量不变（仍为 n/100），A 的输出将接近 n/100。
- 如果 x_S‘ 是错误的，由于两个向量交集很小，追加后不同元素数量会接近 n/50，A 的输出将接近 n/50。
  解码器可以通过 A 的输出区分这两种情况，从而找到正确的原始向量。

这同样定义了一个从大小为 |C| = 2^{cn} 的集合到 S 比特字符串的单射，因此 S ≥ cn，即空间是线性的。

集合 C 的存在性：可以通过概率方法证明。随机选择大量大小为 n/100 的子集，利用切尔诺夫边界论证，可以证明存在一个满足上述条件的指数级大集合。

随机化精确算法也需要线性空间

这个问题将作为习题（PSet 1），其证明思路与上述类似，同样需要利用那个特殊的向量集合 C。

总结：对于不同元素计数问题，要获得亚线性空间（即 o(n) 空间）的算法，我们必须同时接受随机化和近似。这是该问题内在的复杂性所决定的。

转盘模型与线性草图 🌀

现在，让我们转向另一种更通用的流处理模型。到目前为止，我们看到的算法都属于“仅插入”模型。但是，如果流中允许删除操作呢？

什么是转盘模型？

在转盘模型中，我们维护一个初始为零的 n 维向量 x。流是一系列更新操作 (i, Δ)，表示 x_i ← x_i + Δ，其中 Δ 可以是正数（插入）或负数（删除）。我们的目标是计算或近似向量 x 的某个函数 f(x)，同时使用远小于存储整个向量 x（需要 O(n) 空间）的内存。

例子：当所有 Δ = 1 时，f(x) = F0（非零坐标的数量）就是带删除的不同元素计数问题。

核心算法技术：线性草图

设计转盘模型算法的主要（通常是唯一）技术是线性草图。

核心概念：算法在内存中维护一个“草图”向量 y，它是原始向量 x 的一个线性变换：y = Πx。其中 Π 是一个 m × n 的矩阵（m << n）。

处理更新：当收到更新 (i, Δ) 时，新的向量 x’ = x + Δ * e_i（e_i 是第 i 个标准基向量）。因此，新的草图 y’ = Πx’ = Πx + Δ * Π_i，其中 Π_i 是矩阵 Π 的第 i 列。这意味着更新草图 y 只需要加上 Δ 乘以第 i 列即可。关键在于，我们不需要显式存储巨大的矩阵 Π，而只需能够高效生成其任意一列（例如，通过哈希函数指定）。

矩估计问题 📊

现在，我们在线性草图的框架下看一个具体问题：矩估计。我们希望估计向量 x 的 L_p 范数的 p 次方，即 F_p = Σ_i |x_i|^p。当 p=0 时（定义 0^0=0），F_0 就是不同元素计数。

这个问题有一个有趣的复杂度相变：

当 0 ≤ p ≤ 2 时，存在空间复杂度为 O(ε^{-2} log n) 的 (1±ε) 近似算法。
当 p > 2 时，任何常数因子近似算法都需要 Ω(n^{1 - 2/p}) 的空间，并且存在匹配此下界的算法。

`p=2` 时的算法：AMS草图

我们首先看 p=2 的情况，即估计 F_2 = Σ_i x_i^2（L_2 范数的平方）。算法由Alon、Matias和Szegedy在1996年提出。

算法描述：

选择一个 m × n 的矩阵 Π，其每一行 σ_i 是一个从 {1, ..., n} 到 {-1, +1} 的4-wise独立哈希函数。我们可以用 O(log n) 比特存储整个 Π。
维护草图 y = Πx，其中 y_i = σ_i · x（点积）。
当需要估计时，输出 (1/m) * Σ_{i=1}^m y_i^2。

工作原理：

期望：E[y_i^2] = F_2。因为 E[σ_{i,j}^2] = 1，且对于 j ≠ j‘，由于2-wise独立性，E[σ_{i,j} σ_{i,j‘}] = 0。
方差：利用4-wise独立性，可以计算得到 Var[y_i^2] ≤ 2 * F_2^2。
通过设置 m = O(1/ε^2) 并取平均值，再利用切比雪夫不等式，我们可以得到一个 (1±ε) 近似估计。为了进一步提高成功概率，可以取多个估计值的中位数。

这个算法是线性草图的一个完美例子：更新高效，且空间复杂度为 O(ε^{-2} log n)。

扩展到 `0 < p ≤ 2`：Indyk算法 🧪

对于更一般的 p（0 < p ≤ 2），我们可以使用Indyk提出的算法。该算法的核心依赖于一类称为 p-稳定分布 的概率分布。

p-稳定分布的定义：一个分布 D 是 p-稳定的，如果对于任意向量 x ∈ R^n 和独立同分布于 D 的随机变量 Z_1, ..., Z_n，随机变量 Σ_i x_i Z_i 的分布与 (Σ_i |x_i|^p)^{1/p} * Z 相同，其中 Z 也服从分布 D。换句话说，线性组合的分布是 x 的 L_p 范数乘以一个同分布的随机变量。

例子：

p=2：高斯分布是2-稳定的。
p=1：柯西分布是1-稳定的。

一个关键定理（Levy, 1920s）指出：p-稳定分布当且仅当 0 < p ≤ 2 时存在。对于 p > 2，不存在这样的分布。这从算法角度解释了为什么 p > 2 时问题的复杂度会剧增。

理想化算法（假设无限精度和真正的 p-稳定分布）：

设矩阵 Π 的每一行 Z_i 是一个独立同分布的 p-稳定随机向量。
维护草图 y = Πx，其中 y_i = Z_i · x。
估计时，输出 median(|y_1|, |y_2|, ..., |y_m|)。这里我们假设所使用的 p-稳定分布 D_p 满足 P_{Z~D_p}(|Z| ∈ [1-ε, 1+ε]) 是一个常数（可通过缩放实现）。

分析思路（下节课详述）：可以证明，每个 |y_i| 是 ||x||_p 的一个常数因子近似估计器，并且具有常数概率。通过取 m = O(1/ε^2) 个这样的估计值并计算它们的中位数，我们可以得到一个 (1±ε) 近似估计。在实际实现中，我们需要处理有限精度问题，并用有限独立哈希函数来模拟 p-稳定分布的采样。

总结 🎯

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

下界证明：通过巧妙的编码论证，我们证明了对于不同元素计数（F_0）问题，任何使用亚线性空间的算法都必须同时是随机化和近似的。确定性精确、确定性近似、随机化精确的算法都需要线性空间。
转盘模型：一个允许插入和删除操作的通用流模型。
线性草图：解决转盘模型问题的核心技术，通过维护一个原始向量的低维线性投影 y = Πx 来工作。
矩估计 F_p：
- 介绍了 p=2 时高效的AMS草图算法。
- 概述了 0 < p ≤ 2 时基于 p-稳定分布的Indyk算法框架，并指出了 p > 2 时复杂度会发生相变的理论原因（p-稳定分布不存在）。

下一节课，我们将深入分析Indyk算法，并探讨 p > 2 时的情况，然后继续研究更多的流计算问题。

4：Lp范数估计（下）与Nisan伪随机生成器 🧮

在本节课中，我们将继续学习用于估计数据流中向量Lp范数的算法。我们将首先完成对Indyk算法的分析，然后探讨如何利用Nisan的伪随机生成器来减少算法所需的随机比特数。最后，我们将简要介绍一种适用于p>2情况的Lp范数估计算法。

回顾：Indyk算法与p稳定分布

上一节我们介绍了Indyk算法，它利用p稳定分布来估计向量的Lp范数。让我们快速回顾一下核心概念。

一个分布D_p被称为p稳定分布，如果对于任意实数向量x = (x₁, x₂, ..., xₙ)和从D_p中独立抽取的随机变量Z₁, Z₂, ..., Zₙ，其线性组合的分布满足：

∑ᵢ Zᵢ * xᵢ  ~  ||x||_p * D_p

其中，||x||_p 是向量x的Lp范数。这种分布仅当 0 < p ≤ 2 时存在。

算法步骤如下：

构造一个大小为 m × n 的线性投影矩阵 Π，其中每个元素 Πᵢⱼ 独立地从p稳定分布D_p中抽取。
计算投影向量 y = Πx。
估计Lp范数为 median(|y₁|, |y₂|, ..., |yₘ|)。

其中，行数 m 的数量级为 O(1/ε²)，以获得 (1±ε) 的近似。

分析Indyk算法的正确性

本节中，我们来看看如何分析这个算法的正确性。核心思想是证明投影值 yᵢ 的中位数以高概率落在期望的区间内。

我们定义指示函数 I_a,b，当x在区间[a, b]内时值为1，否则为0。

考虑随机变量 I_[ -1-ε, 1+ε ]( yᵢ / ||x||_p )。由于 yᵢ / ||x||_p 服从 D_p 分布，且我们假设D_p的半数质量位于[-1,1]区间内，因此该指示函数的期望值为：

E[ I_[ -1-ε, 1+ε ]( yᵢ / ||x||_p ) ] = 1/2 + Θ(ε)

这意味着，在期望下，超过半数的 yᵢ 满足 |yᵢ| ≤ (1+ε)||x||_p。

同理，对于区间 [-1+ε, 1-ε]，其期望值为 1/2 - Θ(ε)，意味着少于半数的 yᵢ 满足 |yᵢ| ≤ (1-ε)||x||_p。

如果这两个事件同时发生，那么 |yᵢ| 的中位数必然落在 [(1-ε)||x||_p, (1+ε)||x||_p] 区间内，这正是我们想要的。

然而，这只是期望。为了以高概率保证这两个事件发生，我们需要使用切尔诺夫不等式，这要求我们控制随机变量的方差。

指示函数 I_[...]( yᵢ / ||x||_p ) 的方差最多为1。对于m个独立同分布的此类变量，其平均值的方差最多为 1/m。因此，通过设置 m = O(1/ε²)，我们可以应用切尔诺夫界，证明算法以至少2/3的概率成功。

随机性挑战与Nisan伪随机生成器 🎲

上述分析假设矩阵Π的每一行、每一列的元素都是完全独立的随机数。存储这样一个矩阵需要海量的随机比特，这在实践中是不可行的。

一个自然的想法是使用k-独立哈希函数来生成“看似随机”但相关性受控的随机数。理论上，对于p稳定分布，只要k足够大（例如 k = O(1/ε^p)），使用k-独立随机变量构造的矩阵Π仍然能使算法近似正确。

但是，为每一行单独存储一个k-独立哈希函数的种子，总空间开销仍然可能很大。这里的关键洞察是：行与行之间只需要满足两两独立。这允许我们使用一个全局的、空间效率更高的伪随机方案。

接下来，我们将介绍Nisan的伪随机生成器，它是一种专门用于“欺骗”空间受限计算模型的强大工具。

Nisan伪随机生成器简介

Nisan的定理为解决上述随机性问题提供了优雅的方案。它主要针对一种称为分支程序的计算模型。

分支程序可以视作一个网格状的有向图，它模拟了一个空间上限为S比特的算法：

算法从左边的初始节点开始。
算法每次读取S比特的输入。
根据当前状态和读入的S比特，算法转移到下一层的某个节点。
经过R轮后，算法到达最终状态（最后一层的某个节点）。

Nisan定理指出：存在一个伪随机生成器H，它仅使用 T = O(S * log(R/S)) 个真正的随机比特，就能生成一个长度为 R*S 的“伪随机”比特串。对于任何分支程序B和任何输出判断函数F，用H生成的串作为输入，与用真正均匀随机串作为输入，导致B输出特定结果的概率差异极小（不超过 2^{-S}）。

这意味着，对于一个空间为S的算法，外部只需要提供 O(S log R) 个真随机比特，就能模拟出需要 R*S 个真随机比特才能产生的行为，且算法几乎无法区分。

生成器H的构造基于一个完美二叉树：

树根是一个随机的S比特串。
树的每一层使用一个独立的两两独立哈希函数。
左子节点复制父节点的值。
右子节点是父节点值经过该层哈希函数变换的结果。
所有叶子节点的值拼接起来，就构成了输出的伪随机长比特串。

存储整个生成器（包括所有哈希函数描述和种子）仅需 O(S log R) 比特。

将Nisan生成器应用于Indyk算法

现在，我们来看看如何利用Nisan的生成器来为Indyk算法高效地生成矩阵Π。

我们构造一个特定的分支程序B_x（它依赖于输入向量x），该程序模拟Indyk算法的验证过程：

程序维护三个计数器：low, high, cur，初始化为0。
按顺序读入矩阵Π的每一个元素 Πᵢⱼ。
对于每个元素，更新 cur += Πᵢⱼ * xⱼ。
当处理完一行（即一个yᵢ）后，检查 cur：
- 如果 |cur| ≤ (1-ε)||x||_p，则 low++。
- 如果 |cur| ≤ (1+ε)||x||_p，则 high++。
处理完所有行后，如果 high > m/2 且 low < m/2，则程序输出“成功”。

这个程序的空间复杂度 S = O(log n + log m + 精度比特数) = O(log n)。输入长度 R*S 对应整个Π矩阵的所有元素。

我们的分析表明，如果Π是真正随机的，那么这个分支程序以高概率（≥2/3）输出“成功”。根据Nisan定理，如果我们使用Nisan生成器产生的伪随机比特来填充Π，分支程序输出“成功”的概率将与使用真随机时非常接近（差异小于 1/poly(n)）。

因此，最终的Indyk算法实现只需存储一个 O(log² n) 比特的Nisan生成器种子，即可按需动态生成整个矩阵Π的元素，将总空间复杂度控制为 O(ε⁻² log n + log² n)。

对于p > 2的情况

Indyk算法和p稳定分布仅适用于 p ≤ 2。对于 p > 2 的情况，需要不同的方法。

Andoni等人提出了一种基于指数分布和稀疏投影的算法。其核心思想是：

构造一个对角矩阵D，其对角线元素为 1 / (u_i)^{1/p}，其中 u_i 服从指数分布。
构造一个稀疏投影矩阵P，每列只有一个随机±1的非零元素。
维护的草图是 y = P * D * x。
用 ||y||_∞（即y中绝对值的最大值）来估计 ||x||_p。

该算法可以达到 Õ(n^{1-2/p}) 的空间复杂度，这被证明是估计p>2时Lp范数所需的空间下界。

简要分析思路：

首先，可以证明 z = Dx 的最大值 ||z||_∞ 本身就以常数概率给出 ||x||_p 的常数倍近似。
然后，证明稀疏投影P不会显著破坏这个最大值估计的性质。
由于P非常稀疏，最终只需要维护 O(n^{1-2/p}) 个计数器。

总结

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

完成了对Indyk算法的分析：利用p稳定分布和切尔诺夫界，证明了其中位数估计器能以 O(1/ε²) 的空间提供(1±ε)近似的Lp范数估计（p≤2）。
引入了Nisan伪随机生成器：这是一个强大的工具，可以用 O(S log R) 的真随机比特，为空间S的算法模拟 R*S 长的随机输入，从而极大减少算法对随机性的需求。
应用Nisan生成器优化Indyk算法：通过构造一个模拟算法验证过程的分支程序，我们可以使用Nisan生成器动态产生矩阵Π，将存储开销从巨大的完全随机矩阵降至 O(log² n)。
简要了解了p>2的范数估计方法：使用指数分布和稀疏投影，可以在 Õ(n^{1-2/p}) 空间内获得常数倍近似。

这些技术展示了在大数据算法中，如何巧妙地结合概率论、伪随机性和计算模型理论，设计出既高效又实用的解决方案。

5：P>2范数估计与点查询问题 🧮

在本节课中，我们将学习如何估计P>2的向量范数，并初步了解流式模型中的另一个重要问题——点查询（Point Query）及其相关的频繁项（Heavy Hitters）问题。

P>2范数估计回顾

上一节我们介绍了使用稳定分布（Stable Distributions）和Indyk算法来估计P≤2的向量范数。本节中，我们来看看当P>2时，如何设计一个高效的流式算法来估计向量的LP范数。

我们回顾一下Andoni等人提出的算法框架。该算法维护三个部分：一个随机符号矩阵 P，一个对角矩阵 D，以及我们的高维向量 x。

矩阵 P：每一列只有一个随机符号（±1），其余位置为0。
矩阵 D：是一个对角矩阵，对角线元素为 u_i^{-1/p}，其中 u_i 服从指数分布。
向量 Z：定义为 Z = Dx。

上节课我们证明了一个关键引理：Z 的无穷范数 ||Z||_∞ 能以高概率在常数因子内近似 x 的LP范数 ||x||_p。然而，直接存储 Z 需要 O(n) 空间，这在大数据场景下不可行。

因此，算法的核心是维护一个低维的素描（sketch）向量 Y，定义为：
Y = P * D * x
其中 Y 的维度 m 远小于 n。我们的目标是证明，通过输出 ||Y||_∞，我们仍然能得到 ||x||_p 的一个良好近似。

算法直观理解与游戏计划 🎯

Y 的构造可以直观理解为“哈希计数”过程：

高维向量 Z 的每个坐标 Z_i 被一个哈希函数 h 映射到 m 个桶（即 Y 的坐标）中的一个。
在放入桶之前，Z_i 会乘以一个随机符号 σ_i。
每个桶 Y_j 的值就是所有哈希到该桶的 σ_i * Z_i 之和。

这带来了两个挑战：

大值碰撞：如果多个“大”的 Z_i（其值接近 ||x||_p）被哈希到同一个桶，它们可能因随机符号而相互抵消，或者叠加得过大，从而破坏近似精度。
小值干扰：即使大值被完美哈希（即不碰撞），桶中也可能混入许多“小”的 Z_i，这些噪声可能干扰我们对桶中大值的观测。

我们的分析将围绕解决这两个挑战展开：

控制大值数量：我们将证明，值超过 T/(C log n)（其中 T = ||x||_p）的“重”坐标（Heavy indices）的期望数量是 (C log n)^p，即 poly(log n)。通过设置足够多的桶数 m，我们可以以高概率保证这些少量的重坐标被完美哈希（无碰撞）。
控制小值噪声：对于每个桶，来自“轻”坐标（Light indices）的贡献和是一个期望为0的随机变量。我们将使用伯恩斯坦不等式（Bernstein‘s Inequality）证明，以高概率，每个桶中的噪声幅度都很小（例如，小于 T/10）。

基于以上两点，我们可以得出结论：||Y||_∞ 能够以高概率在常数因子内近似 ||x||_p。

严格分析

定义与引理1：重坐标数量的期望

设 T = ||x||_p。对于任意 L > 0，定义重坐标集合 H 为满足 |Z_i| > T/L 的索引 i 的集合。

引理1：重坐标的期望数量 E[|H|] ≤ L^p。

证明概要：
定义指示变量 Q_i = 1 当且仅当 |Z_i| > T/L。
Z_i = x_i / u_i^{1/p}，其中 u_i 服从指数分布 Exp(1)。
因此，
Pr[Q_i = 1] = Pr[|Z_i| > T/L] = Pr[u_i < (|x_i|^p * L^p) / T^p] = 1 - e^{- (|x_i|^p * L^p) / T^p}。
利用不等式 1 - e^{-x} ≤ x (对于 x ≥ 0)，我们有：
Pr[Q_i = 1] ≤ (|x_i|^p * L^p) / T^p。
求和得：E[|H|] = Σ_i Pr[Q_i = 1] ≤ (L^p / T^p) * Σ_i |x_i|^p = L^p。
证毕。

取 L = C log n，则 E[|H|] = O((log n)^p)。根据马尔可夫不等式，|H| 实际为 poly(log n) 的概率很高。如果我们设置桶的数量 m 约为 (log n)^{2p}，根据生日悖论原理，这些重坐标被完美哈希（无碰撞）的概率就很高。

控制轻坐标的噪声

现在，我们关注单个桶 j 中轻坐标的贡献。设 L 为轻坐标集合（|Z_i| ≤ T/L）。桶 j 中的噪声为：
Noise_j = Σ_{i ∈ L, h(i)=j} σ_i * Z_i

我们将对固定的 D（即固定的 Z）应用伯恩斯坦不等式。对于每个轻坐标 i，定义随机变量 R_i = σ_i * Z_i * 1_{h(i)=j}。这些 R_i 在给定 Z 的条件下是独立的（因为哈希函数 h 和符号 σ 独立于 Z）。

期望： E[R_i] = 0。
界 K：由于是轻坐标，|R_i| ≤ |Z_i| ≤ T/L。取 K = T/L。
方差上界：可以证明，Var(Σ_i R_i) ≤ ||Z||_2^2 / m。

然而，||Z||_2^2 本身依赖于随机矩阵 D。我们需要证明它以高概率不大。计算其期望：
E[||Z||_2^2] = Σ_i x_i^2 * E[1 / u_i^{2/p}]。
由于 p > 2，积分 ∫_0^∞ e^{-x} / x^{2/p} dx 收敛到一个常数。因此，E[||Z||_2^2] = O(||x||_2^2)。

接着，利用霍尔德不等式（Hölder‘s Inequality）关联 L2 和 Lp 范数：
||x||_2^2 ≤ ||x||_p^2 * n^{1 - 2/p}。
因此，E[||Z||_2^2] = O(T^2 * n^{1 - 2/p})。

根据马尔可夫不等式，以高概率（例如99%），||Z||_2^2 = O(T^2 * n^{1 - 2/p})。我们选择桶数 m = Θ(n^{1 - 2/p} * log n)，使得 Var(Σ_i R_i) ≤ O(T^2 / log n)。

现在应用伯恩斯坦不等式。我们想证明 Pr[|Noise_j| > T/10] 非常小。将方差上界和 K 值代入伯恩斯坦不等式，并通过调整常数（即增大 m 的常数因子），我们可以使这个失败概率小于 1/n^3。

由于只有 m ≤ n 个桶，通过并集界限（Union Bound），所有桶的噪声都小于 T/10 的概率很高。

结论

结合以上分析：

重坐标数量少，且以高概率被完美哈希。
每个桶中的噪声以高概率很小。

因此，对于包含重坐标的桶，其值 Y_j 主要由该重坐标主导（± 一个小的噪声），从而 |Y_j| 是 T 的一个良好近似。对于不包含重坐标的桶，其值 Y_j 就是噪声本身，很小。所以，||Y||_∞ 就是某个重坐标值的良好近似，进而也是 ||x||_p 的常数因子近似。

空间复杂度为 O(m) = O(n^{1 - 2/p} * log n)。值得注意的是，当 p=2 时，空间为 polylog(n)；当 p>2 时，空间为 n^{1 - 2/p}。这暗示 p=2 可能是一个计算复杂度的临界点，我们将在后续关于下界的课程中深入探讨。

从范数估计到点查询与频繁项 🔍

现在，我们转向流式模型中的另一个实际问题。假设有一个高维向量 x（例如，索引是所有可能的搜索词，值是搜索次数），它在动态更新（流式数据）。我们关心两类查询：

点查询（Point Query）：给定索引 i，返回 x_i 的一个近似值，保证误差在 ±ε * ||x||_1 范围内。
Output ∈ [x_i - ε||x||_1, x_i + ε||x||_1]
频繁项/热点发现（Heavy Hitters）：返回一个列表 L，满足：
- 若 x_i ≥ ε||x||_1，则 i ∈ L。（包含所有真热点）
- 若 x_i < (ε/2)||x||_1，则 i ∉ L。（避免明显的假阳性）

点查询是更基础的操作。事实上，如果能解决点查询，就能解决频繁项发现问题（尽管时间复杂度可能较高）。方法如下：
运行一个误差参数为 ε/10 的点查询数据结构。当需要获取频繁项列表时，遍历所有可能的索引 i（从1到 n），并对每个 i 执行点查询。如果查询结果 ≥ 0.9ε||x||_1，则将 i 加入输出列表。这样可以保证捕获所有真正的频繁项，且列表中的项至少具有 0.8ε||x||_1 的频率。

确定性点查询算法简介

接下来，我们介绍一种基于不相干矩阵（Incoherent Matrix）的确定性方法来解决点查询问题。

定义（ε-不相干矩阵）：
一个 m × n 的矩阵 Φ 是 ε-不相干的，如果：

对任意列 i，||Φ_i||_2 = 1。（列是单位范数）
对任意不同的列 i 和 j，|⟨Φ_i, Φ_j⟩| ≤ ε。（任意两列近似正交）

不相干矩阵可以看作一个“好”的素描矩阵。

定义（ε-纠错码）：
一个参数为 (ε, T, Q, n) 的纠错码 C 是一个包含 n 个码字的集合，每个码字 c_i 是长度为 T、字母表大小为 Q 的向量。并且，对于任意两个不同的码字 c_i 和 c_j，它们的汉明距离（不同位置的个数）至少为 (1 - ε)T。

引理：一个 (ε, T, Q, n) 纠错码可以构造出一个 ε-不相干 矩阵 Φ，其维度为 m = Q * T。

构造方法：
矩阵 Φ 的每一列对应一个码字 c_i。我们将 Φ 的 m 行分成 T 个块，每块有 Q 行。对于第 i 列，在第 t 个块中，我们在对应于码字 c_i 的第 t 个符号的位置上放置一个 1，该块其他位置为 0。最后，将整个列乘以归一化因子 1/√T，使得每列的 L2 范数为1。由于码字间汉明距离大，任意两列的内积至多为 ε。

最终结论（下节课详细证明）：
一个 ε-不相干 矩阵可以被用作一个素描（sketch）来解决点查询问题。并且，存在随机的编码方案，其参数 T = O((1/ε) log n), Q = O(1/ε)，能以高概率构造出所需的 (ε, T, Q, n) 码。因此，我们可以先（随机地）构造这样一个码，进而得到一个确定性的不相干矩阵 Φ，用于解决点查询。

总结 📝

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

P>2范数估计：我们深入分析了Andoni等人的算法，通过维护一个低维素描 Y = P D x，并巧妙地控制“重”坐标的碰撞和“轻”坐标的噪声，实现了对 ||x||_p 的常数因子近似估计，所需空间为 O(n^{1 - 2/p} log n)。
点查询与频繁项问题：我们定义了流式模型中的点查询和频繁项发现问题，并指出点查询是更基础的操作。频繁项发现问题可以通过多次调用点查询来解决。
确定性算法基础：我们引入了不相干矩阵和纠错码的概念，并展示了如何利用纠错码构造不相干矩阵，为下节课讲解基于不相干矩阵的确定性点查询算法奠定了基础。

下节课，我们将完成点查询算法的分析，并探讨更高效的（随机化）频繁项发现算法。

6：点查询与频繁项检测 🎯

在本节课中，我们将要学习大数据流处理中的两个核心问题：点查询和频繁项检测。我们将探讨如何使用不相干矩阵和Count-Min Sketch这两种关键技术来解决这些问题，并比较它们的优缺点。

点查询问题 📍

上一节我们介绍了流处理的基本模型。本节中我们来看看点查询的具体定义。

点查询问题定义如下：

存在一个高维向量 x（维度为 n），在“旋转门”流模型下持续更新。
查询时，输入一个索引 i（1 ≤ i ≤ n），需要输出 xᵢ 的值，并允许一定的误差。
误差要求：输出值必须在 xᵢ ± ε‖x‖₁ 范围内。其中 ‖x‖₁ 是向量的 L1 范数（所有元素绝对值之和）。

这被称为 L1 点查询。

频繁项检测问题 🔍

与点查询紧密相关的问题是频繁项检测（或称“Heavy Hitters”问题）。

频繁项检测问题定义如下：

给定一个参数 ε，需要返回一个索引的列表 L。
完整性：所有“ε-频繁项”（即满足 xᵢ ≥ ε‖x‖₁ 的项）都必须包含在列表 L 中。
准确性：列表 L 中不应包含任何“非 (ε/2)-频繁项”（即满足 xᵢ < (ε/2)‖x‖₁ 的项）。

这意味着列表 L 的大小最多为 2/ε。

使用 ε-不相干矩阵求解 🧮

上一节我们定义了ε-不相干矩阵。本节中我们来看看如何利用它来同时解决上述两个问题。

一个矩阵 Π 是 ε-不相干的，如果其所有列向量都是单位范数，且任意两个不同列向量的内积绝对值不超过 ε。可以将其视为几乎正交的向量集合。

核心思路：使用 Π 作为草图（sketch）。

草图构建：

维护向量 y = Πx。

点查询估计：

当查询索引 i 时，输出估计值 x̃ᵢ = (Πᵀ y)ᵢ，即 Πᵀ 的第 i 行与 y 的点积。

原理分析：

因为 ΠᵀΠ 近似于单位矩阵。
具体推导：x̃ = Πᵀy = ΠᵀΠx。
因此，x̃ᵢ = xᵢ + Σ_{j≠i} (xⱼ * ⟨πᵢ, πⱼ⟩)。
由于 |⟨πᵢ, πⱼ⟩| ≤ ε，所以误差 |x̃ᵢ - xᵢ| ≤ ε * Σ_{j≠i} |xⱼ| ≤ ε‖x‖₁。

这提供了 L∞/L1 的误差保证：‖x̃ - x‖∞ ≤ ε‖x‖₁。

不相干矩阵的构造方法

以下是两种主要的构造方法：

1. 基于随机编码

思路：随机生成一个编码，其码字（codeword）之间的重合度很低。
参数：设置字母表大小 Q ≈ 2/ε，码字长度 T ≈ O(log n / ε²)。
方法：独立随机地为 n 个长度为 T 的向量选择每个位置的符号（1到Q）。
分析：利用切尔诺夫边界和并集界限，可以证明以高概率得到一个 ε-不相干矩阵。其草图大小 M = Q * T = O((log n)/ε²)。

2. 基于里德-所罗门码

思路：使用代数编码（多项式）构造。
参数：Q = T = Θ((log n)/(ε * (log log n + log(1/ε))))。
结果：草图大小 M = Q * T = Θ((log² n)/(ε² * (log log n + log(1/ε))²))。
比较：当 ε 不太小时，这种方法比随机编码更优；当 ε 极小时，两者渐近性能相同。

存在一个下界表明，任何不相干矩阵必须满足 M ≥ Ω((log n)/(ε² * log(1/ε)))。上述构造在某些参数范围内达到了这个下界。

从点查询到频繁项检测的黑盒转换

给定一个误差为 ε/4 的点查询数据结构，可以如下解决频繁项检测问题：

算法步骤：

实例化该点查询数据结构。
遍历所有索引 i（从 1 到 n）。
对每个 i 进行点查询，得到估计值 x̃ᵢ。
如果 x̃ᵢ ≥ (3ε/4)‖x‖₁，则将 i 加入输出列表 L。

正确性分析：

对于真正的 ε-频繁项，其真实值 xᵢ ≥ ε‖x‖₁。由于误差不超过 (ε/4)‖x‖₁，估计值 x̃ᵢ ≥ (3ε/4)‖x‖₁，因此会被包含。
对于非 (ε/2)-频繁项，其真实值 xᵢ < (ε/2)‖x‖₁。加上误差后，估计值 x̃ᵢ < (3ε/4)‖x‖₁，因此不会被包含。

缺点：需要遍历所有 n 个索引，查询时间慢。

随机化方法：Count-Min Sketch 🎲

上一节我们介绍了确定性的不相干矩阵方法。本节中我们来看看一种更节省空间的随机化方法——Count-Min Sketch。

Count-Min Sketch 基于哈希技术，在“严格旋转门”模型（所有 xᵢ ≥ 0）下工作。

数据结构：

创建一个 L 行、T 列的计数器网格 C。
选择 L 个相互独立的、两两独立的哈希函数：h₁, h₂, ..., hₗ，每个都将索引 {1...n} 映射到 {1...T}。

更新操作（将 xᵢ 增加 v）：

对于每一行 r（1 ≤ r ≤ L）：
- 计算哈希位置 b = hᵣ(i)。
- 将计数器 C[r][b] 增加 v。

点查询估计（查询 xᵢ）：

对于每一行 r，查看计数器值 C[r][hᵣ(i)]。
输出这些值中的最小值作为 xᵢ 的估计值。即：x̃ᵢ = min_{1≤r≤L} C[r][hᵣ(i)]。

参数设置与误差分析：

设置 T = ⌈2/ε⌉，L = ⌈log₂(1/δ)⌉。其中 δ 是失败概率。
对于某一行 r，计数器 C[r][hᵣ(i)] 包含 xᵢ 加上所有与 i 哈希到同一位置的其它项 xⱼ 之和。
由于哈希函数是两两独立的，其它任意项 j 与 i 碰撞的概率为 1/T。
因此，误差项的期望值 E[error] ≤ (Σ_{j≠i} xⱼ) / T ≤ (‖x‖₁) / T ≤ (ε/2) ‖x‖₁。
在严格旋转门模型下，所有 xⱼ ≥ 0，误差项非负，因此每个计数器给出的都是高估的值。
根据马尔可夫不等式，单行误差过大的概率：P(error > ε‖x‖₁) ≤ (E[error])/(ε‖x‖₁) ≤ 1/2。
由于 L 行是独立的，所有行同时误差过大（导致最小值误差过大）的概率 ≤ (1/2)^L = δ。

因此，Count-Min Sketch 以至少 1-δ 的概率保证：x̃ᵢ ≤ xᵢ + ε‖x‖₁。

空间复杂度：M = O((1/ε) * log(1/δ))。这比确定性方法的 O(1/ε²) 更优。

关于适应性查询的说明：此随机化保证适用于非适应性查询。如果查询基于之前查询的结果（适应性查询），则此保证可能失效，此时确定性方法更有优势。

基于树的快速频繁项检测算法 🌳

上一节提到黑盒转换的查询时间很慢。本节中我们利用 Count-Min Sketch 和树结构，实现更快的频繁项检测。

核心思想：

将 n 个索引视为完全二叉树的叶子节点。
为树的每一层维护一个独立的 Count-Min Sketch，用于跟踪该层上每个“节点”（即索引区间）的权重和。
在严格旋转门模型下，如果一个叶子是频繁项，那么从该叶子到根节点路径上的所有祖先节点，在其所在层也都是频繁项。

算法过程：

更新：当更新叶子 i 时，更新从叶子 i 到根节点路径上每一层对应的 Count-Min Sketch。
查询：
- 从根节点开始（根节点总是频繁项）。
- 检查当前节点的两个子节点（通过查询该层的 Count-Min Sketch）。
- 只递归进入那些估计值超过 (3ε/4)‖x‖₁ 的子节点。
- 最终到达叶子节点，将其加入输出列表。

正确性与效率：

每一层最多只有 O(1/ε) 个节点可能是频繁项。
因此，在整个树中，我们最多只会检查 O((1/ε) * log n) 个节点。
需要为每个 Count-Min Sketch 设置更小的失败概率 δ‘，通过并集界限保证整体算法成功的概率。

最终复杂度：

空间：O((1/ε) * log n * log(log n / (εδ)))。
查询时间：O((1/ε) * log n)。
更新时间：O(log n * log(1/δ))。

更强的保证：L1/L1 稀疏恢复 🏆

Count-Min Sketch 不仅能做点查询，还能提供更强的稀疏近似恢复保证。

问题定义：给定向量 x，希望找到一个 k-稀疏（即最多 k 个非零元）的向量 y，使得 ‖x - y‖₁ ≤ (1+ε) * ‖x_{-k}‖₁。其中 ‖x_{-k}‖₁ 是将 x 中绝对值最大的 k 个分量置零后的 L1 范数。

使用 Count-Min Sketch 的算法：

设置 Count-Min Sketch 参数：T = O(k/ε)， L = O(log(1/δ))。
对每个索引 i，通过 Count-Min Sketch 得到估计值 x̃ᵢ。
找出 x̃ 中绝对值最大的 k 个索引，构成集合 T。
输出向量 y，其中 yᵢ = x̃ᵢ（若 i ∈ T），否则 yᵢ = 0。

直觉解释：

当哈希桶数量 T 足够大时，一个特定索引 i 与“头部”（绝对值最大的 k 个元素）发生碰撞的概率很低。
因此，估计误差主要来自“尾部”元素与 i 的碰撞，其期望与 ‖x_{-k}‖₁ 成正比，而非整个 ‖x‖₁。

通过分析可以证明，得到的 k-稀疏向量 y 满足上述 L1/L1 近似保证。特别地，如果 x 本身就是 k-稀疏的，那么可以精确恢复 x。

总结 📝

本节课中我们一起学习了大数据流处理中的核心算法：

点查询和频繁项检测的问题定义。
使用 ε-不相干矩阵的确定性解法，其草图大小为 O((log n)/ε²)。
使用 Count-Min Sketch 的随机化解法，其草图大小为 O((log(1/δ))/ε)，效率更高，但保证适用于非适应性查询。
通过树形结构将 Count-Min Sketch 应用于频繁项检测，显著提升了查询速度。
Count-Min Sketch 还能提供更强的 L1/L1 稀疏恢复保证，用于压缩感知等任务。

这些工具为在大数据流中快速、节省空间地监控关键统计量提供了坚实的基础。

7：L0采样与图连通性算法 🧮➡️🌳

在本节课中，我们将学习如何利用L0采样这一强大的数据流处理原语，来解决图论中的连通性问题。我们将从L0采样的基本概念和构建方法开始，然后展示如何将其应用于动态图（支持边插入和删除）的连通分量计算中。

L0采样：基本概念与构建 🎯

上一节我们讨论了在数据流中寻找“重要项”（Heavy Hitters）的算法。本节中，我们来看看一个与之相关但目标不同的原语：Lp采样，特别是L0采样。

什么是Lp采样？

在数据流（尤其是回转式模型，支持插入和删除）中，我们维护一个向量 x。Lp采样的目标是：在流结束时，能够按照特定的概率分布随机抽取一个坐标索引 i。

具体来说，我们希望以正比于 |xi|^p 的概率抽取索引 i。即：
P(采样到索引 i) ∝ |xi|^p

什么是L0采样？

L0采样是Lp采样在 p=0 时的特例。在流计算文献中，“L0范数”通常指向量 x 的支撑集大小（即非零坐标的数量）。因此，L0采样的目标是：均匀地随机抽取一个 x 中的非零坐标。

更正式地说，我们希望算法能以接近均匀的概率输出一个非零坐标的索引，同时允许一定的失败概率。算法的保证如下：

如果 x 的支撑集非空，算法以高概率 (1-δ) 输出一个均匀随机的非零坐标索引。
算法也允许以概率 δ 输出“失败”。

如何构建L0采样器？🔨

我们将通过组合几个更简单的组件来构建L0采样器。

组件A：恢复严格1-稀疏向量

假设我们有一个向量 y，并且我们承诺它在流结束时恰好只有一个非零元素。我们的目标是构建一个数据结构，能够确定性地恢复这个非零元素的值和索引。

解决方案如下：
我们维护两个简单的线性测量（计数器）：

A = Σ yi （所有元素的和）
B = Σ i * yi （索引加权的和）

我们的草图矩阵有两行：第一行全是1，第二行是 [1, 2, 3, ..., n]。
当流结束时，如果 y 确实是1-稀疏的（即 yi 是唯一的非零项），那么：

A = yi
B = i * yi

因此，我们可以恢复：

索引 i = B / A
值 yi = A

这个方案是确定性的，不需要随机性。

组件B：检测非1-稀疏向量

现在，我们需要另一个组件。它的目标是：当向量 y 的支撑集大小不等于1时（即它是零向量或至少包含两个非零项），能够以高概率检测到这一情况并输出“失败”。

有多种方法可以实现这一点，以下是两种思路：

方法一（利用AMS草图）：

运行组件A的恢复算法，得到一个候选向量 y‘。
同时，我们维护一个AMS草图（用于估计L2范数）。
在查询时，我们计算 y‘ - y 的草图（由于草图是线性的，这可以在后处理中完成）。
检查 y‘ - y 的L2范数估计值是否为零。如果不为零，则说明 y 不是1-稀疏的。

方法二（利用随机哈希）：

选择一个随机哈希函数 h: [n] -> {0, 1}，将每个坐标随机映射到两个桶中。
为每个桶分别维护一个AMS草图，用于估计桶内元素的L2范数。
如果 y 的支撑集中至少有两个元素，那么至少有1/2的概率这两个元素会被哈希到不同的桶中。
因此，如果我们发现两个桶的L2范数都大于零，就可以断定支撑集大小至少为2。
通过重复此过程 O(log(1/δ)) 次并取“或”结果，可以将成功概率提升到 1-δ。

组件C：几何采样与最终整合

现在，我们有了能处理1-稀疏向量的组件A和能检测非1-稀疏向量的组件B。如何利用它们从任意向量 x 中进行L0采样呢？关键思想是几何采样。

几何采样步骤如下：

我们创建 log n + 1 个“虚拟流”，对应向量 y0, y1, ..., y_log n。
使用一个哈希函数 H: [n] -> {0, 1, ..., log n}，使得 P(H(i) = j) ≈ 1 / 2^(j+1)。这意味着每个坐标 i 以几何递减的概率被分配到不同层级的虚拟流中。
当流中对 xi 有更新时，我们将其（乘以相应的随机符号后）更新到对应的虚拟流 y_H(i) 中。
核心观察：设 x 的支撑集大小为 T。那么，对于虚拟流 y_j，其期望支撑集大小为 T / 2^(j+1)。因此，必然存在某个层级 j*，使得 y_j* 的期望支撑集大小在常数范围内（例如介于1和2之间）。
对于这个 j*，其支撑集大小恰好为1的概率是一个常数（例如，当采样率为 1/T 时，概率约为 1/e）。

最终算法流程：

我们并行运行所有 log n 个虚拟流，每个流都配备组件A和组件B的数据结构。
在查询时，我们检查所有虚拟流。对于每个流 yj：
- 先用组件B检测它是否不是1-稀疏的。如果是，则拒绝这个流。
- 如果组件B没有拒绝，则用组件A尝试恢复。如果恢复成功，我们就得到了一个候选样本。
由于层级 j* 以常数概率产生一个1-稀疏的向量，我们只需重复整个几何采样框架 O(log(1/δ)) 次，然后选择第一个成功的样本，即可将总体失败概率降至 δ。

空间复杂度：

每个虚拟流需要 O(log(1/δ)) 的计数器（用于组件B的AMS草图等）。
有 O(log n) 个虚拟流。
重复 O(log(1/δ)) 次以降低失败概率。
最终总空间为 O(log n * log²(1/δ)) 个机器字。通过更精细的构造（如使用 s-稀疏恢复和Prony方法），可以优化到 O(log n * log(1/δ))，这已被证明在常数 δ 下是最优的。

应用：动态图连通性 🌉

现在，我们来看如何将L0采样应用于图流算法中。考虑一个无向图 G=(V, E)，边以回转式流的形式出现（即边可以被插入和删除）。我们的目标是维护一个数据结构，使得在流结束后，能够回答任意两个顶点是否连通（即是否在同一个连通分量中）。

挑战与目标

在仅插入边的模型中，一个简单的 O(n) 空间算法是维护一个生成森林（最多 n-1 条边）。然而，对于回转式模型（有删除），这个问题更具挑战性。我们的目标是实现 O(n polylog n) 空间的算法，这被称为“半流式”算法。

图草图算法

核心思想是将图连通性问题转化为L0采样问题。我们为每个顶点 v 维护一个与所有潜在边相关的向量 f(v)。

定义向量 f(v)：

维度是所有可能的边 (a, b)（共 n choose 2 维）。
对于边 e = (a, b)：
- 如果 v 不是 a 或 b，则 f(v)_e = 0。
- 如果边 e 不存在，则 f(v)_e = 0。
- 如果边 e 存在，且 v 是 a 和 b 中较小的那个，则 f(v)_e = +1。
- 如果边 e 存在，且 v 是 a 和 b 中较大的那个，则 f(v)_e = -1。

关键性质：
对于一个顶点集合 S（代表当前的一个连通分量或超节点），定义 f(S) = Σ_{v in S} f(v)。

f(S) 的支撑集恰好是那些一端在 S 内，另一端在 S 外的边（即割边）。因为 S 内部的边对应的 +1 和 -1 会相互抵消。

连通分量算法

我们可以模拟一个经典的“契约-合并”算法：

初始化：每个顶点自成一个分量。
迭代：重复 O(log n) 轮。在每一轮中，对于当前的每个分量 S，找到一条离开 S 的边 (u, v)（其中 u in S, v not in S），然后将 S 与 v 所在的分量合并（通过收缩边 (u, v)）。
正确性：每一轮，非最大的连通分量数量至少减半。O(log n) 轮后，所有连通分量都被正确计算出来。

流式实现

我们无法显式存储所有 f(v)。取而代之的是：

为每个顶点 v 维护一个 L0采样器草图，作用于向量 f(v)。这需要 polylog n 空间。
当需要找到离开分量 S 的边时，我们计算 f(S) 的草图。由于草图是线性的，f(S) 的草图等于 S 中所有顶点草图的和。
然后，我们对这个合并后的草图进行 L0采样。采样结果将以高概率返回一条离开 S 的边（因为 f(S) 的支撑集就是这些边）。
找到边后，合并两个分量，这对应于将它们对应的草图相加，形成新分量的草图。
注意自适应性问题：一轮中采样结果会影响下一轮的图状态。为确保算法正确性，我们需要为每一轮迭代使用独立的L0采样器草图。因此，我们需要预先存储 O(log n) 个独立的草图副本，这使总空间变为 O(n polylog n)。

总结 📚

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

L0采样的定义：从数据流向量中均匀随机抽取一个非零坐标。
L0采样器的构建方法：通过组合确定性1-稀疏恢复、非1-稀疏检测以及几何采样技术来实现。
L0采样在图流算法中的强大应用：通过为每个顶点维护L0采样草图，我们能够以 O(n polylog n) 的空间解决动态图（支持边插入和删除）的连通性问题。

L0采样作为一个基础原语，其应用远不止于此，它还可用于解决最小割、顶点连通性、顶点覆盖乃至最大权匹配等众多图流问题。这体现了流算法中“草图”和“采样”技术的强大与优美。

8：流式模型中的动态规划 🧮

在本节课中，我们将学习如何在流式模型下处理动态规划问题。我们将重点关注两个密切相关的问题：最长递增子序列 和 距离单调性。我们将看到，尽管这两个问题在离线情况下本质相同，但在流式模型下，它们的空间复杂度却截然不同。

上一节我们介绍了流式模型中的一些上界算法，本节中我们来看看如何将经典的动态规划问题适配到流式场景。

问题定义

在以下两个问题中，我们都有一个输入序列，序列中的每个数字介于1到M之间。

最长递增子序列：目标是找到一个尽可能长的子序列，使得其中的数字是非递减的（即弱递增）。
距离单调性：目标是找到需要从序列中删除的最少元素数量，使得剩余序列是单调递增的。

这两个问题是紧密相关的。事实上，距离单调性的值等于序列总长度 n 减去最长递增子序列的长度。因为要得到单调序列，最优策略就是保留最长递增子序列，删除其他所有元素。

流式算法的挑战与成果

标准的动态规划解法需要线性空间来存储DP表，这在流式大数据场景下是不可接受的。因此，研究转向了亚线性空间的近似算法。

以下是该领域的一些关键成果：

LIS的确定性算法：可以在约 O(√n / ε) 字的空间内，得到LIS长度的 (1+ε) 近似。这个上界是紧的，任何确定性算法都需要这么多空间。
DTM的随机算法：可以在仅 O((1/ε) * log² n) 字的空间内，得到DTM的 (1+ε) 近似。这比LIS所需的空间高效得多。
开放性问题：是否存在高效的随机化算法，能在多对数空间内近似LIS？这仍然是一个未解决的问题。

本节课将重点介绍Saks和Seshadhri在2013年提出的，用于解决距离单调性问题的随机流式算法。

基础动态规划解法

首先，我们回顾一下离线情况下，使用动态规划解决距离单调性问题的思路。为了处理边界情况，我们在序列首尾添加哨兵元素（值分别为 -∞ 和 +∞，权重为无穷大，确保它们永远不会被删除）。

定义：

w_i：第 i 个元素的权重（在基础问题中，所有 w_i = 1，哨兵权重为∞）。
W_i：前 i 个元素的权重前缀和（W_i = Σ_{j=1}^{i} w_j）。
S_i：考虑子序列 x_0, x_1, ..., x_i，并且强制保留元素 x_i 时，需要删除的最小权重和（即子问题的DTM值）。

动态规划递推关系如下：
对于每个时间点 t（从1到 n+1）：

更新前缀和：W_t = W_{t-1} + w_t
计算 S_t：S_t = min_{i∈R, x_i ≤ x_t} ( S_i + (W_{t-1} - W_i) )
- 其中 R 是当前已处理的所有索引集合 {0, 1, ..., t-1}。
- 这个递推式的含义是：为了以 t 结尾，我们需要找到一个之前的元素 i（其值不大于 x_t），然后删除 i 和 t 之间的所有元素。
将 t 加入集合 R。

最终，S_{n+1} 就是整个序列的距离单调性值。这个算法需要 O(n) 空间来存储所有的 S_i 值。

流式算法：随机化“遗忘”策略

流式算法的核心思想是：我们无法记住所有过去的 S_i 值。Saks和Seshadhri的巧妙之处在于，在每一步之后，以一定的概率随机“遗忘” DP表中过去的部分条目。

我们维护一个动态的、规模较小的集合 R，它只包含我们“记住”的过去索引。算法流程修改如下：

对于每个时间点 t：

步骤A, B, C：与基础DP类似，但只在当前的小集合 R 上计算最小值，得到近似的 r_t（对应原来的 S_t）。
关键步骤D：对于当前 R 中的每个索引 i，我们以概率 (1 - p_{i,t}) 将其从 R 中移除。

概率 p_{i,t} 的定义：
我们定义 q_{i,t} = min( 1, (1+ε)/ε * ln(4t³/δ) * (w_i / W_{[i, t]}) )，其中 W_{[i, t]} 是区间 [i, t] 内元素的权重和。
那么，p_{i,t} = q_{i,t} / q_{i, t-1}（当 t > i），且 p_{i,i} = 1。

这个设计的精妙之处在于：

q_{i,t} 恰好是元素 i 在时间 t 之后仍然保留在集合 R 中的概率。
当所有权重 w_i = 1 时，q_{i,t} ≈ O(1/(t-i))。这意味着一个元素被记住的概率大致反比于它距今的时间长度。
因此，在任何时刻 t，我们期望记住的过去元素数量大约是 O(log t)，实现了多对数空间。

算法分析概要

我们需要证明两件事：

正确性：算法输出的 r_{n+1} 以高概率是真实DTM值 S_{n+1} 的 (1+ε) 近似。
空间界：集合 R 的大小以高概率保持为多对数级别。

正确性证明的核心思路：

固定一个最优的LIS序列 C。
定义“不安全”的时刻：如果在某个时刻 t，区间 [i, t] 内所有属于 C 的元素都已被我们从 R 中遗忘，则称 i 在 t 时刻是“不安全”的。
可以证明，算法的输出 r_{n+1} 最多不超过 S_{n+1} 加上所有“不安全”元素的权重和。
通过精心设计的概率 p_{i,t}，我们可以证明，任何一个元素成为“不安全”元素的概率非常低。再利用“危险区间”的概念进行聚合分析，最终可以证明，所有“不安全”元素的总权重以高概率不超过 ε * S_{n+1}。这就保证了近似比。

空间界证明：
在任意时刻 t，R 中元素的数量是多个独立伯努利指示变量之和。其期望值是 O(log t)。通过切尔诺夫界和所有时间步上的联合界，可以证明 R 的大小始终以高概率保持在多对数级别。

总结

本节课中我们一起学习了如何将动态规划问题应用于流式模型。我们重点探讨了距离单调性问题，并介绍了一个巧妙的随机算法。该算法通过在每一步以精心设计的概率“遗忘”过去的DP状态，将空间复杂度从线性降低到了多对数级别，同时保证了 (1+ε) 的近似精度。

关键要点在于：

LIS 和 DTM 在流式模型下的复杂度不同，后者存在更高效的空间近似算法。
随机化和近似是突破流式算法空间下界的强大工具。
通过分析“不安全”事件和利用概率工具，可以证明这种“健忘型”动态规划的有效性。

这个算法展示了流式计算中一种重要的设计范式：为了获得亚线性空间，有时需要主动且智能地丢弃信息。

9：流算法空间下界证明 🧠

在本节课中，我们将学习如何证明流算法的空间下界。我们将通过通信复杂性这一工具，将流算法的空间需求与通信问题的难度联系起来，从而证明几个经典流问题所需的最小空间。

通信复杂性简介

上一节我们提到了通信复杂性是证明流算法空间下界的关键工具。本节中，我们来看看它的具体定义。

通信复杂性研究的是两个或多个参与者（例如 Alice 和 Bob）为了共同计算某个函数 F(x, y) 所需交换的最小信息量（比特数）。

参与者：Alice 持有私有输入 x ∈ X，Bob 持有私有输入 y ∈ Y。
目标：共同计算函数 F(x, y) 的值（通常输出为 0 或 1）。
协议：参与者事先约定好一套通信规则（协议）。根据规则和各自的输入，他们交换一系列消息。
成本：通信成本是所交换消息的总比特数。目标是设计成本最低的正确协议。

总是存在两种平凡的协议：

Alice 直接将整个 x 发送给 Bob，成本为 log|X| 比特。
Bob 直接将整个 y 发送给 Alice，成本为 log|Y| 比特。

研究的核心在于，是否存在成本显著低于这些平凡协议的非平凡协议。

通信复杂性与流算法下界的关联

理解了通信复杂性的基本概念后，我们来看看它如何帮助我们证明流算法的空间下界。

其核心思想是归约。如果我们知道某个通信问题 F 很难（即需要高通信成本），并且可以证明“如果能用较小的空间 S 解决某个流问题，就能用较小的通信成本解决 F”，那么流问题的空间 S 就不可能太小。

具体连接方式如下：
假设存在一个空间复杂度为 S 比特的单遍流算法 A。考虑一个通信问题，其中 Alice 的输入构成流的前半部分，Bob 的输入构成流的后半部分。为了计算关于整个流的函数 F，他们可以执行以下协议：

Alice 在她的输入（前半部分流）上运行流算法 A。
Alice 将算法 A 的完整内存状态（S 比特）发送给 Bob。
Bob 从他持有的输入（后半部分流）开始，继续在接收到的内存状态上运行算法 A。
Bob 查询算法 A 的最终输出，得到答案。

这个协议的总通信量正好是 S 比特。因此，任何解决该通信问题的协议都至少需要 S 比特，所以 S 必须至少等于该通信问题的（单向）通信复杂性下界。

确定性精确 F0 的线性下界

为了熟悉这个证明框架，我们首先用它来重新表述一个我们已经知道的结果：精确计算不同元素数量（F0）需要线性空间。

我们通过从相等函数归约来证明。

相等函数：EQ(x, y)，其中 x 和 y 都是 n 比特字符串。当 x = y 时输出 1，否则输出 0。
已知结论：EQ 的单向确定性通信复杂性为 Ω(n) 比特。直观上，Alice 必须发送几乎整个 x，Bob 才能判断是否相等。

归约证明：
假设存在一个空间为 S 比特的确定性流算法 A，可以精确计算 F0。

Alice 将她的 n 比特输入 x 视为一个流：如果 x_i = 1，则将元素 i 加入流；如果 x_i = 0，则不加入。
Alice 在此流上运行算法 A，并将最终的内存内容（S 比特）发送给 Bob。
Bob 现在可以恢复出 Alice 的整个输入 x：对于每个 i 从 1 到 n，Bob 将 i 追加到流中，并再次查询 F0。如果 F0 计数增加，说明 i 之前不在流中（即 x_i = 0）；如果未增加，说明 i 已在流中（即 x_i = 1）。
恢复 x 后，Bob 可以将其与自己的输入 y 比较，计算出 EQ(x, y)。

上述过程构成了一个解决 EQ 问题的、通信成本为 S 比特的单向协议。由于 EQ 需要 Ω(n) 比特通信，因此必然有 S = Ω(n)。这就证明了精确 F0 需要线性空间。

通信复杂性的不同度量

在讨论随机化流算法下界之前，我们需要区分通信复杂性的几种不同度量，它们对应不同的随机性模型。

以下是几个关键定义：

D(F)：计算函数 F 的最优确定性协议的最小成本（最坏情况输入）。
R_δ^{pub}(F)：计算函数 F 的最优随机化协议的最小成本（最坏情况输入），允许的失败概率为 δ，且随机性是公共的（Alice 和 Bob 都能看到同一个随机字符串）。
R_δ^{priv}(F)：与上类似，但随机性是私有的（Alice 和 Bob 各自拥有独立的随机字符串）。
D_δ^μ(F)：当输入 (x, y) 服从分布 μ 时，计算 F 的最优确定性协议的最小成本，允许的失败概率为 δ。

这些度量之间存在以下关系：
D(F) ≥ R_δ^{priv}(F) ≥ R_δ^{pub}(F) ≥ D_δ^μ(F)（对于任意分布 μ）

对于证明随机化流算法的下界，私有随机性模型最贴合实际，因为流算法中的哈希函数等随机性需要存储，在通信中必须显式传递。然而，证明私有模型的下界通常更难。因此，我们通常的策略是：

构造一个“困难”的输入分布 μ。
证明在该分布下，即使有公共随机性，任何协议的通信成本也很高，即 R_δ^{pub}(F) 很大。
由于 R_δ^{priv}(F) ≥ R_δ^{pub}(F)，这自动意味着私有模型的下界也很高。
通过归约，得到流算法的空间下界。

索引问题及其下界

现在，我们引入一个在流算法下界证明中非常核心的通信问题：索引问题。

索引问题：Alice 持有一个 n 比特字符串 x ∈ {0,1}^n。Bob 持有一个索引 j ∈ [n]。他们需要输出 x_j（x 的第 j 位）。
关键结论：索引问题的单向公共随机性通信复杂性下界为 Ω(n)。更精确地说，对于失败概率 δ，有 R_δ^{pub}(Index) ≥ n(1 - H_2(δ))，其中 H_2(δ) 是二元熵函数 -δ log δ - (1-δ)log(1-δ)。当 δ 为常数（如 1/3）时，这也是 Ω(n)。

直观理解：在单向通信中，Alice 不知道 Bob 关心哪一位 j。为了确保 Bob 能以高概率正确猜出任意 j 对应的 x_j，Alice 几乎必须发送整个字符串 x 的信息。

证明概要（利用信息论）：

设 Π 为最优公共随机性协议产生的消息（随机变量）。
协议长度 |Π| 至少等于其熵 H(Π)。
H(Π) ≥ I(X; Π)，即 Π 与整个输入 X 的互信息。
使用链式法则：I(X; Π) = Σ_{i=1}^n I(X_i; Π | X_{<i})。
对于均匀随机的 X，X_i 与 X_{<i} 独立，因此 I(X_i; Π | X_{<i}) = H(X_i) - H(X_i | Π, X_{<i})。
H(X_i) = 1（因为均匀分布）。
由于协议能以至少 1-δ 的概率正确输出 X_i（对任意 i），根据法诺不等式，条件熵 H(X_i | Π) ≤ H_2(δ)。进一步地，H(X_i | Π, X_{<i}) ≤ H(X_i | Π) ≤ H_2(δ)。
因此，每个 i 对应的互信息至少为 1 - H_2(δ)。
总和为 n(1 - H_2(δ))，即 |Π| = Ω(n)。

应用于中位数计算

我们首先展示如何将索引问题归约到精确中位数计算，从而证明后者需要线性空间。

归约构造：
假设存在一个空间为 S 的随机化流算法，能以至少 2/3 的概率精确计算流的中位数。

Alice 持有 x ∈ {0,1}^n，Bob 持有索引 j ∈ [n]。
Alice 构造流：对于 i = 1 to n，她将值 (2i + x_i) 插入流中。例如，若 x_i=0，插入 2i；若 x_i=1，插入 2i+1。这样，Alice 的 n 个比特被编码为 n 个不同的、成对出现的数字。
Alice 运行中位数算法，并将内存状态（S 比特）发送给 Bob。
Bob 扩展流：Bob 的目标是让第 j 个元素（即 2j + x_j）成为整个扩展流的中位数。他通过添加大量“小”值或“大”值来实现：
- 他添加 (j-1) 个值为 0 的元素（所有值都小于 2）。
- 他添加 (n-j) 个值为 (2n+2) 的元素（所有值都大于 2n+1）。
现在，扩展后的流共有 n + (j-1) + (n-j) = 2n-1 个元素。其中位数是第 n 小的元素。由于添加的 0 和 (2n+2) 分别位于排序序列的两端，原始流中的第 j 个元素恰好成为总体的第 n 个元素，即中位数。
Bob 继续运行中位数算法，得到结果。如果结果是 2j，则 x_j = 0；如果结果是 2j+1，则 x_j = 1。

这样，我们就用 S 比特的通信解决了索引问题。由于索引问题需要 Ω(n) 比特的通信，因此 S = Ω(n)。这证明了精确中位数计算需要线性空间。

应用于近似 F0

接下来，我们证明近似计算不同元素数量（F0）到 (1±ε) 因子需要 Ω(ε^{-2}) 空间。这通过从间隙汉明距离问题归约来实现。

间隙汉明距离问题：Alice 和 Bob 分别持有 t 比特字符串 x 和 y。他们被承诺以下两种情况之一：
- HAM(x, y) ≥ t/2 + √t
- HAM(x, y) ≤ t/2 - √t
  他们的目标是区分属于哪种情况。
已知结论：间隙汉明距离的单向通信复杂性是 Ω(t)。并且，可以通过索引问题归约来证明此结论。

关键等式：对于 0/1 向量，不同元素数量满足 2 * F0(x, y) = |x| + |y| + HAM(x, y)，其中 |x| 表示 x 中 1 的个数（权重）。

归约证明：

设 t = Θ(1/ε^2)。假设存在一个空间为 S 的流算法，能以 2/3 概率将 F0 近似到 (1±ε/2) 因子。
Alice 和 Bob 获得间隙汉明距离问题的输入 x 和 y（各 t 比特）。
Alice 构造流：将她的输入 x 视为流：如果 x_i=1，则插入元素 i。
Alice 运行 F0 算法，并将内存状态（S 比特）连同她向量的权重 |x|（需 log t 比特）一起发送给 Bob。
Bob 计算：
- Bob 知道 |y|（他自己的输入）。
- Bob 将他的输入 y 作为流后半部分继续运行算法，得到 F0 的估计值 est，满足 |2*F0_true - 2*est| ≤ ε t（因为 (1±ε/2) 近似意味着 2*F0 的误差在 ±ε t 内）。
- 根据关键等式，汉明距离 HAM(x, y) = 2*F0_true - |x| - |y|。
- Bob 可以计算汉明距离的估计：est_HAM = 2*est - |x| - |y|。
- 这个估计的误差满足 |est_HAM - HAM_true| ≤ ε t。
由于间隙汉明距离问题的两个情况相差 2√t ≈ 2ε t（因为 t = Θ(1/ε^2)），只要估计误差小于 ε t，Bob 就能以高概率区分这两种情况，从而解决间隙汉明距离问题。

因此，解决 F0 近似问题的算法空间 S，必须至少与间隙汉明距离问题的通信下界 Ω(t) = Ω(1/ε^2) 一样大。这就证明了 Ω(1/ε^2) 的空间下界。

总结与展望

本节课中，我们一起学习了如何利用通信复杂性这一强大工具来证明流算法的空间下界。

我们首先建立了通信复杂性与流算法空间下界之间的通用归约框架。然后，我们回顾了如何用该框架证明精确 F0 的线性下界。接着，我们介绍了索引问题及其信息论下界证明，并展示了如何将其归约到精确中位数计算，从而得到后者的线性下界。最后，我们通过从间隙汉明距离问题归约，证明了 (1±ε) 近似 F0 需要 Ω(1/ε^2) 空间。

这些技巧是流算法理论中的核心工具。索引问题及其变体（如多玩家下的不相交问题）可用于证明更多问题的下界，例如 p > 2 时的 F_p 矩估计下界。掌握这种归约思想，对于理解流算法能力的根本限制至关重要。

10：流算法下界与降维技术入门 🚀

在本节课中，我们将完成对流算法空间下界的讨论，并开启一个全新的主题：降维技术。我们将首先回顾如何通过通信复杂性证明流算法的空间下界，然后介绍降维的基本概念及其重要性，特别是著名的Johnson-Lindenstrauss引理。

回顾：通过通信复杂性证明流算法下界 📉

上一节我们介绍了如何通过通信复杂性来证明流算法的空间下界。其核心思路是：首先证明某些通信问题是困难的，然后将这些困难问题归约到我们的流问题上。

以下是几个关键归约的例子：

精确估计不同元素数量（F0）：可以通过从“相等性”问题归约，得到确定性的线性空间下界。
随机近似F0：可以通过从“索引”问题归约到“间隙汉明距离”问题，再归约到随机近似F0，从而得到Ω(1/ε²)的空间下界。
随机中位数查找：同样可以通过归约得到下界。

本节中，我们将完成蓝色部分标注内容的证明，然后结束对流算法的讨论。

流算法下界的进一步证明 🔍

1. 随机精确F0的线性空间下界

断言：随机化的、失败概率为常数的精确F0算法需要Ω(n)的空间。

证明：我们通过从“索引”问题归约来证明。

索引问题：Alice有一个n位字符串x，Bob有一个索引j。目标是让Bob知道x[j]是0还是1。
归约构造：
1. Alice根据她的输入x构造流：如果x[i]=1，则将i放入流中。
2. Alice在流上运行假定的精确F0算法，将算法内存状态连同x中1的个数（即支撑集大小）发送给Bob。
3. Bob将元素j追加到流中，然后查询F0值。
4. 如果F0值增加，则说明j原本不在流中，即x[j]=0；否则x[j]=1。
通信成本：Alice发送了S（算法空间） + log n（支撑集大小）比特。
下界：已知索引问题需要Ω(n)比特的通信。因此，S + log n = Ω(n)，这意味着S = Ω(n)。

这个证明简洁地展示了“索引”问题作为下界工具的威力。

2. 随机近似F0的对数空间下界

定理：随机近似F0算法也需要Ω(log n)的空间。

证明思路：

首先，一个已知定理指出：任何通信问题的随机私有硬币通信复杂度，至少是其确定性通信复杂度的对数。即：R(f) = Ω(log D(f))。
回顾之前证明确定性近似F0下界时，我们构造了一个大小为2^Ω(n)的子集族C，其中每个子集大小固定，且任意两个不同子集的交集很小。
在这个子集族C上，确定性“相等性”测试需要Ω(n)通信（根据鸽巢原理）。
根据上述定理，在C上的随机“相等性”测试则需要Ω(log n)通信。
最后，我们可以将C上的随机“相等性”测试归约到随机近似F0问题（方法类似：Alice和Bob分别插入其集合对应的元素到流中，通过查询F0值是否近似翻倍来判断集合是否相等）。
因此，随机近似F0算法的空间复杂度也必须为Ω(log n)。

结合之前得到的Ω(1/ε²)下界，随机近似F0的完整空间下界是Ω(1/ε² + log n)，这与已知的最优算法匹配。

3. 高阶矩Fp (p>2)的空间下界

对于p>2，近似计算Fp（Lp范数的p次方）需要Ω(n^(1-2/p))的空间，这同样是最优的（忽略对数因子）。

证明思路（基于T玩家不相交问题）：

T玩家不相交问题：有T个玩家，每人有一个n位字符串（视为集合）。承诺两种情况之一：要么所有集合两两不相交；要么存在一个公共元素k，使得每对集合的交集恰好是{k}。目标是区分这两种情况。
下界：该问题的随机通信复杂度下界为Ω(n / T)。这意味着某个玩家必须发送Ω(n / T²)比特。
归约到Fp：
1. 玩家们依次将各自集合中的元素作为流处理，并传递流算法的状态。
2. 在“全不相交”情况下，每个元素频率最多为1，因此Fp ≤ n。
3. 在“存在公共元素k”情况下，元素k的频率为T，因此Fp ≥ T^p。
4. 如果p>2，且T约等于n^(1/p)，那么这两种情况的Fp值相差一个大于2的因子。因此，一个足够好的近似Fp算法可以区分它们。
结论：每个玩家传递的空间S必须满足S = Ω(n / T²)。代入T = n^(1/p)，得到S = Ω(n^(1 - 2/p))。

关于L∞范数的注记：L∞范数（最大频率）本质上对应p = log n的情况。通过两玩家不相交问题的归约，可以证明其需要Ω(n)的空间，这与上述公式当p=log n时（n^(1 - 2/log n) ~ 常数）所暗示的线性下界精神一致。

新篇章：降维技术 📉

现在，我们告别流算法，开启新的章节——降维技术。

降维的动机

在许多算法问题中，输入是高维向量（例如，图像特征、文档向量化表示）。高维度会导致：

运行时间增加
存储成本上升
通信开销变大

降维的目标：在运行核心算法之前，对高维数据进行预处理，将其映射到低维空间，同时保留关键的结构信息，从而使得后续算法更高效。

需要保留的几何结构可能包括：

点对之间的距离
向量之间的角度
点集的体积
特定优化问题的最优解性质

本节我们将主要关注距离保持。

失真度与JL引理

定义（失真度）：给定两个度量空间(X, d_X)和(Y, d_Y)，以及映射f: X → Y。若存在常数C1, C2，使得对所有x, x‘ ∈ X，有：
C1 * d_X(x, x') ≤ d_Y(f(x), f(x')) ≤ C2 * d_X(x, x')
且 D = C2 / C1，则称映射f的失真度最多为D。

我们将专注于赋范向量空间，其中距离由范数定义：d(x, y) = ||x - y||。

Johnson-Lindenstrauss (JL) 引理：对于任意一组n个点（向量）{x1, ..., xN} ⊂ R^n（欧几里得空间L2），以及任意ε ∈ (0, 1/2)，存在一个线性映射π: R^n → R^m，其中m = O(ε^(-2) log N)，使得该映射到低维空间后，所有点对之间的欧氏距离保持(1±ε)因子不变。即，映射的失真度至多为(1+ε)。

重要性：JL引理表明，对于欧氏距离，我们可以将点集降至仅与点数量对数相关的维度，而与原始维度n无关！这为处理高维数据提供了强大工具。

为什么专注于L2范数？

研究表明，L2范数在降维中具有特殊性：

对于L1范数，要实现常数失真度降维，目标维度需要是n^Ω(1)，而非对数维度。
更一般地，Johnson和Naor的定理指出：如果一个赋范空间具有“任意n点集都可嵌入到O(log n)维子空间且失真度为常数”的性质（即JL性质），那么该空间本质上必须近似于L2空间。

因此，L2范数是实现高效降维的“最佳舞台”，这也是我们重点研究它的原因。

分布式JL引理

为了证明JL引理，我们通常先证明一个更简单的“分布式”版本。

分布式JL引理：对于任意维度n，任意ε, δ ∈ (0, 1/2)，存在一个分布在矩阵π上（π是m×n矩阵，m = O(ε^(-2) log(1/δ))），使得对于任意固定单位向量x (||x||=1)，有：
Pr[ | ||πx||² - 1 | > ε ] ≤ δ.
也就是说，随机选取的矩阵π以高概率(1-δ)近似保持单个向量的长度。

如何从分布式JL推出标准JL：

给定N个点，考虑所有(N choose 2)个归一化的差向量 (xi - xj)/||xi - xj||。
令分布式JL中的失败概率δ < 1/(N choose 2)。
根据联合界，随机矩阵π在所有差向量上同时保持长度的概率 > 0。
因此，存在一个确定的矩阵π满足要求，其维度m = O(ε^(-2) log N)。

与流算法的联系：你可能会觉得这个随机矩阵π很眼熟。它正是AMS草图等流算法中使用的随机符号矩阵（每个元素独立地以1/√m缩放为+1或-1）。在流算法中，我们通常用切比雪夫不等式分析，得到m = O(1/(ε²δ))。为了得到更好的对数依赖log(1/δ)，我们需要更强大的工具来分析这种随机二次型的尾部行为。

下节课预告 🎯

本节课我们一起学习了流算法下界的最终证明，并引入了降维技术的基本概念和JL引理。

下节课我们将：

证明分布式JL引理。
引入Hanson-Wright不等式，它是分析独立随机变量二次型尾部概率的强大工具。
利用Hanson-Wright不等式，我们可以更优雅地证明分布式JL，并探索更高效的JL构造方法。

我们将看到，大数据算法中的不同领域（如流算法、降维）在技术层面上是紧密相连的。

11：Johnson-Lindenstrauss引理与Hanson-Wright不等式 🎯

在本节课中，我们将学习Johnson-Lindenstrauss引理，这是一个在大数据降维中至关重要的工具。我们将看到如何通过一个简单的随机线性映射，将高维数据点嵌入到低维空间，同时几乎保持所有点对之间的距离。为了证明这一点，我们将深入探讨Hanson-Wright不等式，这是一个强大的概率工具，用于控制随机二次型的集中性。

概述 📋

Johnson-Lindenstrauss引理指出，对于任意一组高维点，存在一个到低维空间的线性映射，能够以高概率保持所有点对之间的距离。我们将首先形式化地定义这个引理，然后介绍一个更易于处理的“分布化”版本。证明的核心将依赖于Hanson-Wright不等式，该不等式为随机符号向量（Rademacher向量）的二次型提供了尾概率界。

Johnson-Lindenstrauss引理与分布化版本

设 X 是 N 维空间中的一个点集，其大小为 |X| = N。对于任意 ε ∈ (0, 1/2)，存在一个线性映射 π: R^N -> R^M，其中 M = O(ε^(-2) log N)，使得对于所有 x, y ∈ X，有：

(1 - ε) ||x - y|| ≤ ||π(x) - π(y)|| ≤ (1 + ε) ||x - y||

这里的 M 远小于 N，实现了降维。

然而，直接证明这个引理是复杂的。一个更实用的途径是证明其“分布化”版本。

分布化Johnson-Lindenstrauss引理：对于任意 ε, δ ∈ (0, 1/2)，存在一个随机矩阵 π 的概率分布（π ∈ R^{M×N}），使得对于任意单位范数向量 x（即 ||x|| = 1），有：

Pr[ | ||π(x)||^2 - 1 | > ε ] ≤ δ

其中 M = O(ε^(-2) log(1/δ))。

这个引理意味着映射 π 能以高概率保持向量 x 的范数。由于 π 是线性的，这自然就保持了所有向量间的距离。通过设置 δ = 1/N^2 并对所有可能的向量对取并集，可以从分布化版本推导出原始的JL引理。

上一节我们介绍了JL引理及其分布化形式的目标，本节中我们来看看如何构造这样的随机映射并证明其性质。

历史构造与简单随机映射

在Johnson和Lindenstrauss于1984年的原始论文中，他们通过随机旋转矩阵后投影到前 M 个坐标来构造 π。然而，随机旋转矩阵的条目不是独立的，实现起来有一定开销。

后续的研究表明，使用具有独立同分布条目的随机矩阵同样有效，这大大简化了实现。我们将分析以下这种简单的构造：

令 π 是一个 M×N 的矩阵，其每个条目 π_{ij} 独立地以 1/√M 的概率取 +1 或 -1。即：

π_{ij} = σ_{ij} / √M, 其中 σ_{ij} 是独立的Rademacher随机变量（等概率取 ±1）。

我们将证明，这个简单的分布满足分布化JL引理。

Hanson-Wright不等式：核心工具

证明分布化JL引理的关键是Hanson-Wright不等式。这个定理由Hanson和Wright于1971年证明，它控制了随机二次型的偏差。

Hanson-Wright不等式：设 A 是一个 n×n 的实对称矩阵。设 σ = (σ_1, ..., σ_n) 是独立的Rademacher随机变量。那么，对于任意 λ > 0，有：

Pr[ | σ^T A σ - E[σ^T A σ] | > λ ] ≤ C_1 * exp(-C_2 * min( λ^2 / ||A||_F^2, λ / ||A||_{op} ))

其中：

||A||_F 是 A 的Frobenius范数（所有条目平方和的平方根）。
||A||_{op} 是 A 的算子范数（最大奇异值，对于实对称矩阵即最大特征值的绝对值）。
C_1, C_2 是普适常数。

这个尾界是混合型的，包含了高斯尾（λ^2 / ||A||_F^2）和指数尾（λ / ||A||_{op}）。它等价于对所有 p ≥ 1 的矩有界：

( E[ | σ^T A σ - E[σ^T A σ] |^p ] )^(1/p) ≲ √p * ||A||_F + p * ||A||_{op}

从Hanson-Wright推导分布化JL

现在，我们利用Hanson-Wright不等式来证明我们构造的随机映射 π 满足分布化JL引理。

证明：
对于单位向量 x，考虑随机变量 Z = ||π(x)||^2。

Z = ||π(x)||^2 = ∑_{r=1}^{M} ( (1/√M) ∑_{i=1}^{N} σ_{ri} x_i )^2 = (1/M) ∑_{r,i,j} x_i x_j σ_{ri} σ_{rj}

我们可以将其重写为二次型。定义一个大向量 σ，它由所有 σ_{ri}（r=1..M, i=1..N）拼接而成。构造一个分块对角矩阵 A_x，其每个 M×M 块都是 (1/M) * x x^T。可以验证：

Z = σ^T A_x σ

并且 E[Z] = 1（因为 E[σ_{ri}σ_{rj}] = 1 当 i=j，否则为0，且 ||x||^2=1）。

我们需要 bound Pr[ |Z - 1| > ε ]。这正是Hanson-Wright不等式的形式，其中 A = A_x。

计算 A_x 的范数：

Frobenius范数：||A_x||_F^2 = 1/M。
算子范数：||A_x||_{op} = 1/M。

将这两个范数代入Hanson-Wright不等式，我们得到：

Pr[ |Z - 1| > ε ] ≤ C_1 * exp( -C_2 * min( ε^2 / (1/M), ε / (1/M) ) )
                 = C_1 * exp( -C_2 * min( M ε^2, M ε ) )

为了使这个概率小于 δ，我们只需要选择 M 使得：

M ≥ C * max( ε^(-2) log(1/δ), ε^(-1) log(1/δ) )

由于 ε < 1/2，ε^(-2) 项占主导，因此 M = O(ε^(-2) log(1/δ)) 足以满足要求。证毕。

可以看到，一旦有了Hanson-Wright不等式，分布化JL引理的证明变得非常简洁。接下来，我们将深入探讨Hanson-Wright不等式本身的证明。

证明Hanson-Wright不等式：所需工具

Hanson-Wright不等式的证明需要三个核心工具：

Khintchine不等式：用于控制随机线性形式的矩。
解耦（Decoupling）：将二次型问题约化为线性形式问题。
平方根技巧（Rudelson's Square Root Trick）：用于处理递归的矩 bound。

我们将逐一介绍并证明它们。

1. Khintchine不等式

Khintchine不等式给出了Rademacher线性组合的矩 bound。

定理（Khintchine）：对于任意向量 x ∈ R^n 和任意 p ≥ 1，存在常数 C_p（仅依赖于 p），使得：

( E[ | ∑_{i=1}^n σ_i x_i |^p ] )^(1/p) ≤ C_p * ||x||_2

特别地，对于 p ≥ 2，有 C_p ≲ √p。这意味着其尾概率服从高斯型衰减：Pr[ |∑ σ_i x_i| > λ ] ≲ exp(-c λ^2 / ||x||_2^2)。

证明思路：我们只需对偶数 p 证明 C_p ≲ √p 即可（通过取整技巧可推广到所有 p）。将期望展开为单项式求和。对于存活项（所有索引出现偶数次），其系数在Rademacher情形下为1，在高斯情形下（将 σ_i 替换为标准高斯 g_i）至少为1。因此，Rademacher的 p 阶矩被高斯情形的 p 阶矩所控制。而高斯线性组合 ∑ g_i x_i 本身是一个方差为 ||x||_2^2 的高斯变量，其 p 阶矩为 ≲ √p * ||x||_2。这就完成了证明。

2. 解耦（Decoupling）

解耦技术允许我们将一个二次型的矩 bound 用两个独立的随机向量序列的线性型的矩来 bound。

定理（Decoupling）：设 A 是 n×n 实矩阵（对角线元素可为任意值）。设 σ 和 σ‘ 是两组独立的Rademacher向量。那么，对于任意 p ≥ 1，有：

( E_σ[ | ∑_{i≠j} A_{ij} σ_i σ_j |^p ] )^(1/p) ≤ 4 * ( E_{σ, σ'}[ | ∑_{i,j} A_{ij} σ_i σ‘_j |^p ] )^(1/p)

注意右边是对所有 i, j 求和，且 σ 和 σ‘ 独立。

证明思路：引入独立的伯努利(1/2)随机变量 η_i。将原和式写为 4 * E_η[ ∑_{i≠j} A_{ij} σ_i σ_j η_i (1-η_j) ]。利用Jensen不等式将期望移到绝对值 p 次方外面。然后，存在一个固定的 η’ ∈ {0,1}^n 使得 bound 成立。令 S = {i: η‘_i = 1}，则和式变为 ∑_{i∈S, j∉S} A_{ij} σ_i σ_j。由于 {σ_j : j∉S} 是独立随机符号，我们可以将其替换为另一个独立的序列 σ‘_j 而不改变分布。最后，通过添加期望为零的项（如 i,j ∈ S 或 i,j ∉ S 的项，因为 E[σ‘_j]=0），我们可以将和式补全为 ∑_{i,j} A_{ij} σ_i σ‘_j，再次利用Jensen不等式完成证明。

3. 证明Hanson-Wright不等式

现在，我们结合Khintchine不等式和解耦技术来证明Hanson-Wright不等式。我们证明其矩形式：对于 p ≥ 2，

( E[ | σ^T A σ - E[σ^T A σ] |^p ] )^(1/p) ≲ √p * ||A||_F + p * ||A||_{op}

（对于 1 ≤ p < 2 的情形，可由 p=2 的情形推出）。

证明：
令 Z = σ^T A σ - E[σ^T A σ]。注意 E[σ^T A σ] = trace(A)，而去心二次型 Z 等价于 ∑_{i≠j} A_{ij} σ_i σ_j（当 A 对称时）。

应用解耦：
```
|| Z ||_p ≲ || ∑_{i,j} A_{ij} σ_i σ‘_j ||_p
```
这里 ||·||_p 表示 L_p 范数（即 (E|·|^p)^(1/p)），且右边的期望是对 σ 和 σ‘ 取的。
条件化与Khintchine：将 σ‘ 视为条件，对 σ 应用Khintchine不等式：
```
|| ∑_{i,j} A_{ij} σ_i σ‘_j ||_p ≤ √p * ( E_{σ‘}[ || A σ‘ ||_2^p ] )^(1/p)
```
因为对于固定的 σ‘，∑_i (∑_j A_{ij} σ‘_j) σ_i 是关于 σ 的线性形式，其 L_2 系数向量的范数正是 || A σ‘ ||_2。
矩的转换：注意 (E[ || A σ‘ ||_2^p ])^(1/p) = || || A σ‘ ||_2^2 ||_{p/2}^{1/2}。通过广义三角不等式，这小于等于：
```
( E[|| A σ‘ ||_2^2 ] )^{1/2} + || || A σ‘ ||_2^2 - E[|| A σ‘ ||_2^2] ||_{p/2}^{1/2}
```
第一项等于 ||A||_F（Frobenius范数）。第二项是另一个去心二次型的 L_{p/2} 范数的平方根，即 || σ‘^T (A^T A) σ‘ - E[...] ||_{p/2}^{1/2}。
递归与平方根技巧：我们对这个新的二次型再次应用解耦和Khintchine（步骤类似，但阶数变为 p/2）。经过推导，我们会得到一个关于变量 E = ||Z||_p / √p 的不等式，其形式类似于：
```
E^2 ≲ ||A||_F + p^{1/2} ||A||_{op}^{1/2} * E
```
这是一个关于 E 的二次不等式。解这个不等式（取较大的根），并利用 √(a^2+b^2) ≲ a + b，最终可得：
```
E ≲ ||A||_F + p^{1/2} ||A||_{op}
```
因此，
```
||Z||_p = √p * E ≲ √p * ||A||_F + p * ||A||_{op}
```
这正是我们想要证明的矩 bound。通过矩生成函数或直接代入尾概率公式，即可得到定理陈述中的指数尾 bound。

总结 🏁

在本节课中，我们一起学习了：

Johnson-Lindenstrauss引理：一个强大的降维工具，它保证可以通过随机线性映射将高维数据压缩到低维空间，同时近乎完美地保持距离信息。
分布化JL引理：一个更易于证明的版本，它要求随机映射能以高概率保持单个单位向量的范数。
简单随机映射构造：使用独立同分布的 ±1/√M 随机变量作为矩阵条目的映射，即可满足分布化JL引理。
Hanson-Wright不等式：证明的核心，它给出了Rademacher二次型集中性的精确控制。
关键证明工具：我们深入探讨了Khintchine不等式（控制线性形式）、解耦技术（将二次型转化为线性型）以及平方根技巧（处理递归矩 bound），并利用它们完成了Hanson-Wright不等式的证明。

通过这一系列内容，我们不仅理解了一个重要的大数据算法理论基础，也掌握了一套分析随机矩阵和二次型集中性的有力工具。这些工具在机器学习、数据科学和算法设计的多个领域都有广泛应用。

12：Johnson-Lindenstrauss引理的下界与超越最坏情况分析 🎯

在本节课中，我们将学习Johnson-Lindenstrauss引理的下界，并探讨如何超越最坏情况分析，以获得针对特定数据集的更优维度约简保证。

1. 下界介绍 📉

上一节我们证明了分布式的Johnson-Lindenstrauss引理。本节中，我们来看看其下界，即证明JL引理在某种意义上是最优的。

这些下界的形式如下：对于任意给定的点数N，总存在一个包含N个点的集合，使得任何嵌入到低维空间的方法都需要至少大约 min(N, 1/ε² log N) 的维度。这表明JL引理几乎是“最坏情况”下最优的。

2. 分布式JL引理的下界

首先，我们区分两个不同的问题：分布式JL引理和标准JL引理。文献中常统称为JL，但它们是不同的问题，其中一个已解决，另一个仍有开放性问题。

分布式JL引理的下界由J. Ram和Woodruff等人证明。其结论是：对于任意ε, δ ∈ (0, 1/2)，任何(ε, δ)-分布式JL方案所需的维度m必须满足：
m ≥ min(n, 1/ε² log(1/δ))

以下是证明思路的简要概述：

分布式JL保证：对于所有向量x，映射π满足 | ||πx||² - 1 | > ε 的概率至多为δ。
这意味着对于球面上的任意分布μ，有：
Pr_{x~μ}[Pr_π[| ||πx||² - 1 | > ε] ] ≤ δ
通过交换概率顺序，可以论证：如果π只有m行，那么除非m足够大（至少为 1/ε² log(1/δ)），否则它无法以高概率保持随机向量的范数。

3. 标准JL引理的下界（线性映射）

如果我们要求嵌入映射F必须是线性的（即F(x) = Πx），那么存在一个下界。

与Casper Green Larson合作的工作表明，对于线性映射，所需的维度m必须满足：
m ≥ min(n, 1/ε² log n)
这个下界是紧的，因为单位矩阵可以达到n维，而JL引理可以达到 1/ε² log n 维。

证明思路是使用概率方法：对于一个固定的线性映射Π，如果它的行数不足，那么它无法以高概率保持一个随机高斯向量的范数。通过考虑一组随机高斯向量，并对其所有可能的距离（包括到0点的距离）进行保持，可以论证，如果Π不能以良好概率保持其中一个向量的范数，那么它保持所有向量范数的概率就更小。最后，通过对所有可能的Π矩阵进行并集论证（需要处理无限多矩阵的问题，可通过构造有限网来解决），可以得出不存在满足条件的线性映射。

4. 标准JL引理的下界（任意映射）

如果我们允许嵌入映射F是任意（非线性）函数，情况则不同。Alon在2003年证明了一个下界：
m ≥ min(n, (1/ε² log n) / log(1/ε))
目前仍有一个开放性问题：能否移除分母中的 log(1/ε) 项，使下界与线性映射情况下的 1/ε² log n 匹配？

证明基于“不相干向量”的概念。考虑一个简单的点集：X = {0, e1, e2, ..., en}，即零向量和标准基向量。任何保持这些点距离的嵌入F，在经过平移使F(0)=0后，可以定义一个矩阵V，其列向量为 vi = F(ei)。根据距离保持性质，可以推导出这些vi具有单位范数，且彼此之间的点积为O(ε)。换句话说，我们得到了一组近似正交（不相干）的向量。

核心问题转化为：在低维空间中，不可能存在太多这样的不相干向量。定义矩阵 A = V^T V，其对角线元素为1，非对角线元素绝对值≤ ε。Alon的结论是，任何如此接近单位矩阵的实对称矩阵A，其秩r必须满足：
r ≥ min(n, (1/ε² log n) / log(1/ε))

为了证明这个秩的下界，我们首先证明一个较弱的版本（Alon-）：如果非对角线元素≤ 1/√n，则秩r = Ω(n)。然后，通过考虑矩阵A的逐元素k次幂 A^(k)，并选择合适的k使得 ε^k ≤ 1/√n。可以证明，如果原矩阵A的秩为r，则 A^(k) 的秩至多为 C(r+k-1, k)。结合弱版本结论，A^(k) 的秩必须为Ω(n)，从而推导出关于r的下界。

关于秩 rank(A^(k)) ≤ C(r+k-1, k) 的证明，利用了行空间基的展开和“星与条”的组合计数论证。

5. 超越最坏情况分析：Gordon定理 🌟

最坏情况的下界告诉我们存在难以处理的向量集。但在实际应用中，我们面对的数据集可能并非这种最坏情况。我们希望能获得依赖于数据集本身特性的、更优的维度约简保证。

这就是Gordon定理（或称“逃离网格定理”）的内容。它提供了针对特定集合T的维度约简保证。

定理陈述：假设投影矩阵Π的每个元素独立服从 ±1/√m 分布（或任何衰减速度不低于高斯分布的、均值为0、方差为1的分布）。令T是单位球面上的一个子集。如果投影后的维度m满足：
m ≥ (G(T)² + log(1/δ)) / ε²
那么，以至少 1-δ 的概率，对于所有 x ∈ T，有：
| ||Πx||² - 1 | ≤ ε

其中，G(T) 是集合T的高斯均值宽度，定义为：
G(T) = E_g [ sup_{x∈T} | g · x | ]
这里，g是一个n维标准高斯随机向量。

理解G(T)：

意义：G(T) 衡量了集合T的“复杂度”或“宽度”。集合越分散、越复杂，G(T) 越大。
与JL的关系：可以证明，对于任意有限集T，总有 G(T) ≤ O(√log |T|)。因此，Gordon定理的要求 m ≥ O( (log |T| + log(1/δ))/ε² ) 永远不会比JL引理差。
优势：如果集合T中的向量高度相关或聚集在一起，G(T) 可能远小于 √log |T|，从而允许使用更低的维度m。例如，如果T中所有向量都几乎相同，则 G(T) 接近一个常数。

6. 如何估计高斯均值宽度 G(T) 📏

为了应用Gordon定理，我们需要估计或界定 G(T)。以下是几种方法，从简单到精细：

简单界：对于任何有限集T，有 G(T) ≤ O(√log |T|)。证明使用了积分尾部概率和 Gaussian 的尾部界。

通过ε-网：如果能为T构造一个大小为N‘的ε-网T’，那么 G(T) ≤ G(T') + ε√n。当向量聚集时，可以用很小的N‘构造网，从而得到更好的界。
Dudley不等式（链式积分）：这是更精细的界。
G(T) ≤ C * ∫_0^∞ √log N(T, u) du
其中 N(T, u) 是T在L2距离下、半径为u的最优覆盖网的大小。这个积分上界通常比简单界更紧。
最优刻画（Talagrand）：存在一个称为 γ₂(T) 的泛函（通过最优分层序列定义），满足：
c * γ₂(T) ≤ G(T) ≤ C * γ₂(T)
即，高斯均值宽度在常数因子内等于 γ₂(T)。这给出了估计 G(T) 的最优方式。

7. Gordon定理与JL引理的等价性

最近的研究表明，分布式JL引理与Gordon定理在本质上是等价的。

Gordon定理显然蕴含了JL引理（因为 G(T) ≤ √log |T|）。
反过来，也可以证明：如果存在一个满足特定参数的分布式JL方案，那么它可以用来推导出Gordon定理的结论。基本思想是，在由Dudley不等式或 γ₂ 泛函定义的不同尺度（scales）上，对数据集T的近似网（nets）应用JL引理，最终可以保证整个集合T的距离被保持。

总结 📚

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

JL引理的下界：证明了在最坏情况下，JL引理给出的目标维度 O((log N)/ε²) 几乎是紧的。我们分别探讨了分布式JL、线性映射JL和任意映射JL的下界。
超越最坏情况分析：引入了Gordon定理，它将维度约简所需的维度与数据集本身的复杂度 G(T)（高斯均值宽度）联系起来。对于结构简单、向量相关性高的数据集，我们可以获得比标准JL引理更好的维度约简保证。
高斯均值宽度的估计：介绍了多种估计 G(T) 的方法，从简单界到最优的 γ₂ 泛函刻画。
定理间的联系：了解到Gordon定理与JL引理在本质上是等价的，这深化了我们对维度约简理论的理解。

通过本讲内容，我们认识到维度约简不仅存在普遍的最坏情况限制，还可以根据数据的具体特性进行优化，这为实际应用提供了更强大的理论工具。

13：超越最坏情况分析与快速JL变换 🚀

在本节课中，我们将学习如何对Johnson-Lindenstrauss引理进行超越最坏情况的分析，并探索如何构建更快速的JL变换矩阵，以加速降维过程。

回顾：高斯过程与Gordon定理

上一节我们介绍了JL引理的下界，并转向了超越最坏情况的分析。我们以Gordon定理的陈述作为结尾。在深入Gordon定理之前，我们定义了高斯过程上确界的一些概念。

设T是R^N中的一个有界子集。我们定义了关于距离函数d的T的γ₂泛函，其形式如下：

γ₂(T, d) ≈ inf sup Σ 2^{r/2} * d(T, T_r)

其中，T_r是T的一个大小为2^{2r}的网。这个泛函与高斯平均宽度G(T)密切相关，G(T)定义为：

G(T) = E[ sup_{x∈T} g·x ]

其中g是一个随机高斯向量。定理指出，γ₂(T)和G(T)在常数因子内是相等的。

Gordon定理是JL引理的一个强化。它指出，如果T是单位球面上的一个子集，那么存在一个具有高斯条目的矩阵Π，其行数依赖于(G(T)/ε)²，使得以高概率同时保留T中的每个向量。

从分布JL到Gordon定理

Gordon定理可以通过分布JL来证明。关键在于满足一个多尺度的分布JL条件。

以下是证明的核心思路。我们首先定义一个引理，该引理列出了几个确定性条件。如果这些条件成立，那么矩阵Π就能保留T中的所有向量。

条件一与条件二（范数保留）：对于所有R，以及所有在特定集合（如T_R、T_{R-1}及其差集）中的向量v，要求Πv的范数接近v的范数。这本质上是标准的JL陈述。
条件三（内积保留）：对于所有在特定集合中的向量u和v，要求Πu和Πv的内积接近u和v的内积。这可以通过将内积表示为向量和与差的范数平方，并应用JL引理来保证。
条件四（算子范数界）：要求矩阵Π的算子范数有界。这可以通过在R^N的一个网上应用JL引理并取并集来证明。

如果矩阵Π的分布能满足所有这些不同尺度上的分布JL条件，那么通过联合界，所有这些事件同时以高概率成立，从而Gordon定理得证。有趣的是，对于亚高斯矩阵（如随机符号矩阵），只需满足基尺度（R=0）的分布JL，就能自动满足所有其他尺度的条件，这直接推导出了Gordon定理。

实现更快的降维：稀疏与结构化JL

到目前为止，我们主要关注如何将目标维度M降至最小。然而，降维过程本身（即计算Πx）也可能成为瓶颈，特别是当Π是稠密矩阵时。本节我们来看看如何加速这一步。

有两种主要思路：

使Π稀疏：如果Π是稀疏的，那么矩阵-向量乘法会更快。
使Π具有快速乘法结构：即使Π不稀疏，特定的结构（如基于哈希）也能实现快速计算。

我们将首先关注稀疏化的方法。

稀疏JL变换的历史与构造

早期的尝试是让Π的条目独立且稀疏。例如，Achlioptas提出了一种分布，其中每个条目以1/6的概率为+1，1/6的概率为-1，2/3的概率为0。这实现了常数因子的加速，但后续研究表明，在条目独立的条件下，加速程度存在基本限制。

为了突破这个限制，需要放弃条目独立性。Kane和Nelson提出了一种构造，其中每一列恰好有S个非零元，每个非零元是一个随机符号除以√S。另一种等价的构造是CountSketch矩阵：将行分成S个块，每个块大小为M/S；对于每一列，在每个块中随机选择一个位置放入一个随机符号。

分析这种稀疏JL变换可以使用Hanson-Wright不等式。关键步骤是将||Πx||² - 1写成一个二次型σᵀB_xσ，其中σ是随机符号向量，B_x是一个依赖于随机稀疏模式Δ的矩阵。然后，通过控制B_x的算子范数和Frobenius范数，并应用矩方法和马尔可夫不等式，可以证明其满足分布JL性质，所需行数M为O((1/ε²) log(1/δ))，每列非零元数S为O((1/ε) log(1/δ))。

总结

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

如何利用高斯平均宽度和Gordon定理对JL变换进行超越最坏情况的、与数据实例相关的分析。
分布JL性质如何蕴含Gordon定理。
为了加速降维计算，可以构造稀疏的JL变换矩阵（如基于CountSketch），在保持目标维度M近乎最优的同时，显著减少计算量。其分析工具仍然是Hanson-Wright不等式和矩方法。

下一讲，我们将继续探讨具有快速乘法结构的JL变换，并可能触及当前的一些开放问题。

14：稀疏JL与快速JL 🚀

在本节课中，我们将完成对稀疏JL（Johnson-Lindenstrauss）嵌入矩阵的证明分析，然后探讨快速JL变换，并介绍JL在高维近似最近邻搜索中的一个应用。

完成稀疏JL证明 🔍

上一节我们介绍了稀疏JL矩阵 Π。现在，我们来完成其证明。

Π 是一个缩放因子为 1/√s 的矩阵。每一列恰好有 s 个非零元素，这些非零元素的位置是随机选择的，并带有随机符号，其余位置为零。各列之间相互独立。

我们关注的关键量是 ||Πx||² - 1 的尾部概率。这可以写成一个二次型 σᵀ B(x) σ，其中 B(x) 是一个块对角矩阵。通过 Hanson-Wright 不等式，我们可以先对随机符号 σ 取期望，得到 √p 倍的 Frobenius 范数加上 p 倍的算子范数。

上一节我们已经论证过，无论 δ（决定非零元素位置的随机变量）如何，B(x) 的算子范数始终不超过 1/s。因此，分析的重点在于处理 Frobenius 范数。

Frobenius 范数的平方可以展开为：
||B(x)||_F² = (1/s²) * Σ_{i≠j} x_i² x_j² * Q_{ij}
其中，Q_{ij} 是一个随机变量，表示坐标 i 和 j 在哈希过程中发生碰撞的次数。

为了完成证明，我们需要一个关键引理。

引理：对于参数 s, m, p 的特定设置，Q_{ij} 的 p 范数至多为 p。

假设该引理成立，我们可以继续分析。利用 p 范数的性质和三角不等式，我们可以推导出：
|| ||Πx||² - 1 ||_p ≤ O(√p/m + p/s)

通过马尔可夫不等式，最终可以证明，对于设定的参数，||Πx||² 偏离 1 超过 ε 的概率至多为 δ。

现在，我们来证明这个关于 Q_{ij} 的引理。

引理证明：
Q_{ij} 是坐标 i 和 j 的碰撞次数。固定列 i，它有 s 个非零位置。列 j 独立于列 i，但列 j 内部的 s 个非零位置选择并非完全独立（因为是无放回抽样）。我们定义指示变量 Y_t，表示在列 i 的第 t 个非零位置上，列 j 是否也非零。那么 Q_{ij} = Σ_{t=1}^{s} Y_t。

虽然 Y_t 之间不独立，但我们可以论证其任意阶矩都被独立情况下的矩所控制。在独立情况下，每个 Y_t 为 1 的概率是 s/m。而在我们的无放回抽样模型中，给定前 t-1 次碰撞后，第 t 次碰撞的条件概率为 (s-t+1)/(m-t+1) ≤ s/m。因此，任何单项式 E[Π Y_{i_q}^{d_q}] 的期望值都小于等于独立情况下的 (s/m)^L。由于所有单项式非负，整个 E[(Σ Y_t)^p] 的期望也被独立情况所控制。

对于独立的指示变量和，我们可以使用切尔诺夫或伯恩斯坦不等式来得到其矩的界，从而证明引理中 Q_{ij} 的 p 范数界。

为什么需要恰好 s 个非零？
如果每列非零个数的期望是 s 而非恰好 s，那么在展开 ||Πx||² - 1 时，对角线项不会完全抵消，这将破坏我们整个二次型 B(x) 的构造，使得分析失效。

快速JL变换 ⚡

稀疏JL的嵌入时间与向量 x 的支撑集大小乘以 s 成正比。当 x 稠密时，我们希望有接近线性的嵌入时间。快速JL变换（FJLT）利用结构化矩阵实现了这一点。

FJLT 矩阵 Π 被构造为三个矩阵的乘积：Π = P * H * D。

以下是各组成部分：

D：一个 n×n 的对角矩阵，对角线元素为随机 ±1 符号。
H：一个 n×n 的正交矩阵，需要满足两个性质：
1. 最大条目幅度为 O(1/√n)。
2. 支持快速矩阵向量乘法（例如，O(n log n) 时间）。
  离散傅里叶变换（DFT）或哈达玛变换（Hadamard Transform）都满足这些条件。
P：一个 n×n 的“采样”矩阵。为了简化分析，我们可以将其视为一个随机对角矩阵，其中每个对角线元素以概率 m/n 为 1（并乘以缩放因子 √(n/m)），否则为 0。这相当于从 H 的结果中随机选择约 m 行。

设计直觉：
核心思想源于海森堡不确定性原理的启发：一个向量和它的傅里叶变换不能同时在少数坐标上高度集中。

直接对 x 采样（Px）方差可能很高，特别是当 x 本身集中在少数坐标时。
先应用 H 可以“铺开”向量，使 Hx 的坐标更均匀，从而降低采样方差。
然而，如果 x 本身已经很均匀，H 可能反而使其集中。随机符号矩阵 D 的作用就是防止这种最坏情况。可以证明，对于任何单位范数向量 x，HDx 的每个坐标的幅度以高概率不超过 O(√(log(n/δ)/n))，从而保证了良好的“铺开”性质。

通过这种预处理，再应用采样矩阵 P，并结合伯恩斯坦不等式等工具进行分析，可以证明 m = O(ε^{-2} log(1/δ) log(n/δ)) 足以以高概率保证JL性质。这比最优的 O(ε^{-2} log(1/δ)) 多了一个 log(n/δ) 因子。

为了获得精确的JL界，可以在 PHD 之后再乘一个具有正确参数的JL矩阵 G（如高斯矩阵或稀疏JL矩阵）。这样，总运行时间为应用 D（线性时间）、H（n log n 时间）、P（线性时间）以及将 G 应用于降维后向量的时间。

后续的研究通过不同的分析框架（例如，先论证 PH 具有限制等距性等性质，再分析 D 的作用）改进了这个结果，消除了对 δ 的对数平方依赖，对于需要 δ 极小的应用场景尤为重要。

应用：高维近似最近邻搜索 🎯

JL引理的一个直接应用是解决高维空间中的近似最近邻搜索问题。

问题定义：
给定一个静态的 n 个点集合 P ⊂ R^d（例如，图像特征向量），需要构建一个数据结构。对于查询向量 q ∈ R^d，返回一个点 p ∈ P，使得 ||p - q|| 至多是真正最近邻距离的 c 倍（c > 1 为近似比）。

维数灾难：
在低维（如2维）中，可以通过构建维诺图（Voronoi Diagram）和点定位在近线性空间和对数查询时间内解决精确最近邻问题。但在高维 d 下，维诺图本身的复杂度可能高达 n^{⌈d/2⌉}，使其不实用。

利用JL的解决方案：
一种思路是将近似最近邻搜索约简为“cR-近邻”问题：给定半径 R，如果存在点距离 q 在 R 以内，则返回一个距离在 cR 以内的点；否则报告失败。通过二分搜索 R，可以将近似最近邻问题转化为 cR-近邻问题。

以下是一个基于JL和网格的简单 cR-近邻数据结构（以 R=1, c=1+ε 为例）：

预处理：
- 选取一个JL变换 Φ: R^d → R^k，其中 k = O(ε^{-2} log n)。
- 对数据库中的所有点 p_i 应用 Φ，得到降维后的点集 Φ(p_i)。
- 在 R^k 空间中，施加一个边长为 ε/√k 的均匀网格。任何网格单元的直径（最大L2距离）不超过 ε。
- 对于每个降维后的点 Φ(p_i)，考虑以其为中心、半径为1的L2球。找出所有与该球相交的网格单元。
- 将所有这样的（网格单元，点索引i）对存储在一个哈希表中。
查询：
- 对于查询 q，首先计算其降维表示 Φ(q)。
- 确定 Φ(q) 落入哪个网格单元。
- 在哈希表中查找该网格单元。如果找到，则返回关联的任意一个点索引 i 对应的原始点 p_i。

正确性：
如果存在某个数据库点 p* 满足 ||Φ(p*) - Φ(q)|| ≤ 1（由JL性质，这近似意味着 ||p* - q|| 较小），那么 Φ(q) 所在的网格单元必定与以 Φ(p*) 为中心、半径为1的球相交，因此会被存储在哈希表中。查询返回的点 p 满足 ||Φ(p) - Φ(q)|| ≤ 1 + ε（因为网格单元直径至多 ε），再次由JL性质，这意味着 ||p - q|| ≈ (1+ε) * ||p* - q||。

空间复杂度：
关键在于每个半径为 1+ε 的球（包含所有与之相交的网格单元）所覆盖的网格单元数量。这个数量由球的体积与网格单元体积的比值决定，约为 O(1/ε)^k。代入 k = O(ε^{-2} log n)，总空间约为 O(dn) + n * (1/ε)^{O(ε^{-2} log n)} = O(dn) + n^{O(1/ε^2)}。

虽然这个空间复杂度对于理论分析是多项式，但指数中 1/ε^2 使其在实际中可能很大。因此，更流行的解决方案是使用局部敏感哈希系列方法，它们能以接近线性的空间和次线性的查询时间解决这个问题。

总结 📝

本节课我们一起学习了：

完成了稀疏JL嵌入的证明，核心在于分析哈希碰撞次数 Q_{ij} 的矩，并利用其被独立情形所控制的性质。
介绍了快速JL变换，它通过 PHD 的结构化矩阵乘积，在 O(n log n) 时间内实现嵌入，特别适用于稠密向量。
探讨了JL在高维近似最近邻搜索中的一个应用，展示了如何通过JL降维和网格划分来构建一个简单的数据结构，并分析了其原理和复杂度。

这些内容体现了降维技术在应对高维数据挑战中的核心作用，并为后续学习更高效的近似最近邻算法（如局部敏感哈希）奠定了基础。

15：近似矩阵乘法 🧮

在本节课中，我们将学习如何高效地近似计算两个大矩阵的乘积。当矩阵的行数 N 远大于列数 D 和 P 时，精确计算 A^T B 的复杂度为 O(NDP)。我们将探讨两种近似方法，它们能以更快的速度得到一个误差可控的结果。

问题定义与精确计算

我们有两个矩阵：A 是 N x D 矩阵，B 是 N x P 矩阵。目标是计算矩阵乘积 A^T B。

直接计算 A^T B 的标准算法时间复杂度为 O(NDP)。我们也可以利用快速矩阵乘法技术。例如，可以将输出矩阵划分为 D x D 的块，每个块的计算可以利用快速方阵乘法算法，其时间复杂度为 O(D^ω)，其中 ω 是矩阵乘法的指数（目前已知 ω < 2.373）。

然而，本节课我们关注的是近似计算。我们愿意接受一个矩阵 C，使得 A^T B - C 的某种范数（例如Frobenius范数）很小，但计算速度更快。

基于采样的方法（Frobenius范数误差）

第一种方法是基于行采样。注意到 A^T B 可以表示为外积的和：A^T B = Σ_{k=1}^{N} a_k b_k^T，其中 a_k 和 b_k 分别是矩阵 A 和 B 的第 k 行。

采样方法的核心思想是：我们并不对所有 N 个外积求和，而是根据某种概率分布采样 M 个行对，并对采样的外积进行缩放，以得到一个无偏估计量。

以下是具体步骤：

定义采样概率：对于第 k 行，其被采样的概率 p_k 与 ||a_k|| * ||b_k|| 成正比。
独立采样 M 次。对于第 t 次采样，若选中第 k_t 行，则构造估计项 Z_t = (1/(M * p_{k_t})) * a_{k_t} b_{k_t}^T。
最终的估计矩阵为 C = Σ_{t=1}^{M} Z_t。

可以证明，C 是 A^T B 的无偏估计。通过马尔可夫不等式分析，可以得出以下结论：当采样次数 M ≥ 1/(ε²δ) 时，以至少 1-δ 的概率满足 ||A^T B - C||_F ≤ ε * ||A||_F * ||B||_F。计算复杂度从 O(NDP) 降为 O(MDP)，而 M 是一个与 ε 和 δ 相关的常数。

改进运行时间对失败概率 δ 的依赖

上述方法的运行时间依赖于 1/δ。我们希望将其改进为 log(1/δ)。一个技巧是运行 R = O(log(1/δ)) 个独立的上述估计器，每个的失败概率设为常数（例如1/10）。根据切尔诺夫界，其中至少有2/3的估计器是“好的”（即误差很小）。

接下来的挑战是如何从这 R 个候选矩阵中识别出一个好的，而无需计算真实的 A^T B。方法如下：

对于每个候选矩阵 C_i，计算它与所有其他 C_j 的Frobenius距离。
如果某个 C_i 与超过半数的其他 C_j 距离都很小（≤ (ε/2)*||A||_F||B||_F），则返回这个 C_i。

分析表明，一个好的 C_i 一定会通过这个测试，而一个坏的 C_i 则无法欺骗我们，因为它要与超过半数的其他矩阵接近，而这些矩阵中必然包含好的，通过三角不等式可以推出这个坏的 C_i 本身误差也很小。这样就以 O(log(1/δ)) 的代价获得了高成功概率。

基于JL引理的方法（Frobenius范数误差）

第二种方法利用满足JL矩性质的随机投影矩阵 Π，这种方法对输入是** oblivious **（无需知晓数据内容）的。

一个分布 D 在 M x N 矩阵上满足 (ε, δ, p)-JL矩性质，如果对于任意单位范数向量 x ∈ R^N，从 D 中随机抽取的 Π 满足：E[ | ||Πx||² - 1 |^p ] ≤ ε^p * δ。

许多熟悉的构造具有此性质，例如：

AMS草图（p=2，需要 M = O(1/(ε²δ))）。
带符号的Count Sketch/Thorup-Zhang草图（p=2，需要 M = O(1/(ε²δ))）。
高斯或亚高斯投影（p=O(log(1/δ))，需要 M = O(log(1/δ)/ε²)）。

主要定理是：如果 Π 来自一个满足 (ε, δ, p)-JL矩性质的分布（p≥2），那么对于任意矩阵 A, B，以至少 1-δ 的概率有：||A^T B - (ΠA)^T (ΠB)||_F ≤ 3ε * ||A||_F * ||B||_F。

证明思路是通过马尔可夫不等式，将问题转化为控制 ||A^T B - (ΠA)^T (ΠB)||_F 的 p 阶矩。通过巧妙的分解和JL矩性质的定义，可以完成这个上界的估计。这种方法的好处是投影矩阵 Π 与输入 A, B 无关，适用于数据流等动态场景。

引入更强的近似：谱范数误差与子空间嵌入

以上两种方法都保证了Frobenius范数下的误差。接下来我们引入一种更强的近似保证：谱范数（算子范数）误差。

我们希望找到 Π，使得 || (ΠA)^T (ΠB) - A^T B ||_op ≤ ε * ||A||_op * ||B||_op。通过将矩阵 A 和 B 并排拼接成一个新矩阵 M = [A B]，可以证明上述问题等价于要求 Π 同时保留下述所有向量的范数：

对于所有向量 x，满足 ||Π(Ax)||² = (1 ± ε) ||Ax||²。

这比谱范数误差的条件更强。满足这个性质的 Π 被称为矩阵 A 列空间的子空间嵌入。子空间嵌入在加速最小二乘回归、主成分分析等许多线性代数问题中至关重要。

下一讲我们将探讨实现子空间嵌入的两种方法：一种基于更复杂的采样（杠杆得分采样），另一种基于JL类型的投影。需要注意的是，要达到这种更强的保证，所需的样本或投影维度 M 将依赖于矩阵 A 的秩或列数 D，而不再是一个常数。

总结

本节课我们一起学习了大数据下的近似矩阵乘法。

我们首先介绍了基于行采样的方法，可以在Frobenius范数下快速得到近似结果，并通过“中位数技巧”将运行时间对失败概率的依赖从 1/δ 改进到 log(1/δ)。
接着，我们介绍了基于JL矩性质的随机投影方法，同样获得Frobenius范数误差保证，并且具有对输入数据 oblivious 的优点。
最后，我们引出了更强大的谱范数误差保证和子空间嵌入的概念，这将是下一讲的重点，并为众多大规模线性代数问题提供加速基础。

16：子空间嵌入 🧮

在本节课中，我们将要学习子空间嵌入的概念。这是一种强大的降维技术，能够高效地近似高维矩阵的运算，并应用于诸如最小二乘回归等问题。我们将探讨其定义、存在性、应用以及两种高效的构造方法：采样法和随机投影法。

子空间嵌入的定义 📖

上一节我们介绍了近似矩阵乘法的概念，本节中我们来看看一个更强大的工具：子空间嵌入。

给定一个 RN 中的线性子空间 E，以及一个 M x N 的矩阵 π。如果对于所有属于 E 的向量 x，矩阵 π 都能保持其范数，即满足以下条件，则称 π 是 E 的一个 ε-子空间嵌入：

对于所有 x ∈ E，有：
(1 - ε) ||x||² ≤ ||πx||² ≤ (1 + ε) ||x||²

这可以看作是 Johnson-Lindenstrauss 引理在特定子空间上的应用。对我们而言，子空间 E 通常是某个矩阵 A 的列空间。这意味着对于所有向量 z（它可以表示为 Ax），其范数都能被 π 近似保持。

子空间嵌入的存在性与维度 🎯

一个自然的问题是：我们需要多少行（即 M 多大）才能构造一个子空间嵌入？

以下是关于子空间嵌入维度的两个关键事实：

不存在性：对于任何 ε < 1，如果 M < D（其中 D 是子空间 E 的维度），则不存在 ε-子空间嵌入。
存在性：总存在一个完美的（误差为0的）子空间嵌入，其行数 M 恰好等于 D。

证明概要：

对于（1），如果 M < D，则矩阵 π 的秩小于 D。其核空间与子空间 E 的交集中必然存在非零向量，该向量被 π 映射为零，无法保持范数。
对于（2），我们可以通过奇异值分解（SVD）来构造。对于任意矩阵 A，其 SVD 为 A = U Σ Vᵀ，其中 U 的列构成了 A 列空间的一组标准正交基。此时，令 π = Uᵀ，即可得到一个具有 D 行且零误差的完美子空间嵌入。

然而，通过 SVD 计算 Uᵀ 需要 O(ND²) 的时间，这在大数据场景下可能过于昂贵。因此，我们的目标是寻找能快速计算和应用的子空间嵌入。

子空间嵌入的应用：最小二乘回归 📉

在了解了子空间嵌入的基本概念后，我们来看看它的一个核心应用：加速最小二乘回归。

我们面临的问题如下：
给定一个高大的 N x D 矩阵 A（N 很大）和一个向量 b ∈ Rᴺ，我们希望找到：

x_ls = argmin over x ∈ Rᴰ ||Ax - b||²

传统的解法（如计算 SVD 或 AᵀA）需要 O(ND²) 的时间。

使用子空间嵌入进行加速：
假设我们有一个对于 A 的列空间以及向量 b 所张成子空间的 ε-子空间嵌入 π（一个 M x N 的矩阵，且 M ≈ D）。我们可以转而求解一个维度小得多的问题：

x_tilde = argmin over x ∈ Rᴰ ||πAx - πb||²

由于 πA 只是一个 M x D 的矩阵，求解这个新的回归问题只需要 O(MD²) ≈ O(D³) 的时间。

关键结论：
可以证明，这样得到的近似解 x_tilde 的质量是有保证的：

||A x_tilde - b||² ≤ (1+ε)/(1-ε) * ||A x_ls - b||²

也就是说，近似解的目标函数值在最坏情况下也不会比最优解差太多。

现在，整个算法的耗时取决于三部分：

找到子空间嵌入 π 的时间。
计算 πA 和 πb 的时间。
求解降维后回归问题的时间（O(D³)）。

我们的目标是让前两步也非常高效。

构造子空间嵌入的方法 🔧

接下来，我们将探讨两种构造子空间嵌入的方法，它们分别对应于我们上节课介绍的两种近似矩阵乘法技术。

方法一：基于杠杆得分的采样法

这种方法的思想是对矩阵 A 的行进行非均匀采样。

杠杆得分：
首先，我们定义矩阵 A 第 i 行的杠杆得分 lᵢ 为：

lᵢ = aᵢᵀ (AᵀA)⁻¹ aᵢ

其中 aᵢ 是 A 的第 i 行（视为列向量）。杠杆得分具有重要的几何意义：它等于将 A 的列空间正交化后，第 i 行范数的平方。所有行的杠杆得分之和等于矩阵 A 的秩 D。

采样过程：
我们独立地采样矩阵的行。采样第 i 行的概率 pᵢ 被设置为与 lᵢ 成正比（但需进行截断和缩放，具体公式涉及 log D 和 1/ε² 因子）。被采样的行会被一个权重因子 1/√pᵢ 进行缩放，以确保无偏性。

理论保证：
通过精心设置采样概率，可以证明，如果我们采样大约 O(D log D / ε²) 行，那么所得到的采样矩阵 π 以高概率是一个 ε-子空间嵌入。

挑战：
计算精确的杠杆得分本身需要 A 的 SVD，这又回到了我们最初想要避免的昂贵计算。不过，存在方法可以快速近似杠杆得分（例如，使用接下来要介绍的随机投影法进行预处理），这可以作为习题进行探索。

方法二：基于随机投影的JL法

这种方法利用随机投影来构造子空间嵌入，其灵感来源于 Johnson-Lindenstrauss 引理。

核心思想：
如果我们使用一个随机高斯矩阵或随机符号矩阵作为 π，那么根据 Gordon 定理的分析，只需要 M = O(D/ε²) 行，就能以高概率保证 π 是任意一个 D 维子空间的 ε-子空间嵌入。

优势与问题：

优势：M 很小（约 D/ε²），且构造 π 完全不需要查看数据 A，速度极快。
问题：计算 πA 是稠密矩阵乘法，需要 O(MND) ≈ O(ND²/ε²) 时间，这可能比原问题更慢。

解决方案：快速 oblivious 子空间嵌入
我们需要一个既能快速生成，又能快速与 A 相乘的 π。这引出了 oblivious 子空间嵌入 的概念。

定义：
一个 OSED 是一个矩阵分布 D，使得对于任意具有标准正交列的 N x D 矩阵 U，从 D 中随机抽取的矩阵 π 满足：

Pr[ || (πU)ᵀ(πU) - I || > ε ] ≤ δ

关键在于，这个分布 D 不依赖于具体的 U（或 A）。

研究表明，我们可以使用稀疏的JL变换或快速JL变换来实例化这样的分布 D。这些变换矩阵 π 非常稀疏或具有快速算法结构（如傅里叶变换），使得计算 πA 的时间可以接近 O(nnz(A))，即矩阵 A 中非零元素的数量，这对于稀疏矩阵是极大的加速。

总结与前瞻 🚀

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

子空间嵌入的定义：它是一种能保持某个子空间内所有向量范数的降维变换。
子空间嵌入的存在性与最小维度要求（至少为子空间维度 D）。
子空间嵌入在最小二乘回归中的应用，可将计算复杂度从 O(ND²) 降至近 O(D³)，前提是能快速获得嵌入矩阵。
构造子空间嵌入的两种主要方法：
- 基于杠杆得分的采样法：需要采样约 O(D log D / ε²) 行，但计算精确杠杆得分成本高。
- 基于随机投影的JL法：需要约 O(D/ε²) 行，使用稠密随机矩阵时乘法成本高，但可通过稀疏或快速的 oblivious 子空间嵌入（如稀疏JL）来解决，实现近乎线性的计算时间。

在下节课中，我们将更深入地探讨 oblivious 子空间嵌入，并介绍使用子空间嵌入加速回归的其他技巧（例如，仅使用常数精度嵌入配合迭代优化）。我们还将看到子空间嵌入在其他问题（如主成分分析PCA）中的应用。

17：获取与使用子空间嵌入 🧮

在本节课中，我们将继续探讨无意识子空间嵌入。我们将学习获取它们的方法，以及如何利用子空间嵌入进行回归分析。此外，我们还将简要介绍低秩近似。

获取OSE的方法 🔍

上一讲我们提到了一种方法，实际上还有多种其他方法可以获取无意识子空间嵌入。以下是五种主要方法：

网论证
使用非交换Khintchine不等式
矩方法
使用具有Frobenius误差的近似矩阵乘法
链化方法

接下来，我们将逐一简要介绍这些方法。

1. 网论证

网论证的基本思想是：对于RN中的任意d维子空间E，存在一个单位范数向量子集T，其大小为O(1)^d。如果投影Π能将T中所有向量的范数保持到(1 ± O(ε))以内，那么它就能将整个子空间E的范数保持到(1 ± ε)以内。

这意味着，分布式的Johnson-Lindenstrauss引理可以自动推导出无意识子空间嵌入。通过设置失败概率δ为1/|T|，并利用联合界，我们可以确保保留T中的所有向量，从而保留整个E。

公式表示：
分布JL（失败概率δ' = δ/|T|） ⇒ OSE（失败概率δ）

2. 非交换Khintchine不等式方法

此方法利用矩阵集中不等式。我们希望以下不等式以高概率成立：
[ |\Pi U^T \Pi U - I| > \epsilon ]
通过马尔可夫不等式，这最多是(1/ε^p)乘以该算子范数p次方的期望。通过选择p至少为log d，我们可以利用霍尔德不等式将Schatten p范数与算子范数联系起来。

核心技巧是对称化：引入独立同分布的副本Y_i，并利用行的独立性，将问题转化为可以应用非交换Khintchine不等式的形式。这种方法已被多篇论文使用，能够获得接近最优的行数m和稀疏度s。

近似结果示例（独立条目）：

m ≈ d * polylog(d/(δ)) / ε²
s ≈ polylog(d/(δ)) / ε

3. 矩方法

矩方法从基本原理出发，计算以下概率：
[ P(|\Pi U^T \Pi U - I| > \epsilon) \leq \frac{1}{\epsilon^p} E[|\Pi U^T \Pi U - I|^p] ]
对于偶数p，矩阵的p次幂的迹可以展开为所有环路上矩阵元素乘积的和。通过计算这个期望，并进行组合计数，可以得到界。这是经典随机矩阵理论的方法，虽然计算量大，但技术复杂性较低。

历史应用：
此方法被Wigner、Freddie Comossh等人用于分析随机矩阵的算子范数。对于稀疏JL，早期通过矩方法得到的结果在log因子数量上不如后来的方法优化。

4. 近似矩阵乘法（AMM）方法

此方法基于一个关键观察：如果矩阵的Frobenius范数小于ε，那么其算子范数也小于ε。
[ |\Pi U^T \Pi U - I|_{op} \leq |\Pi U^T \Pi U - I|_F ]
因此，我们希望：
[ P(|\Pi U^T \Pi U - I|_F > \epsilon) \leq \delta ]
注意到，ΠU^TΠU 接近 U^TU 正是具有Frobenius误差的近似矩阵乘法问题。我们已知，满足JL矩性质的分布都能实现AMM。通过设置ε' = ε/d，我们最终需要 m = O(d²/(ε²δ)) 行。

优势：
此方法允许我们使用任何满足JL矩性质的分布（例如，稀疏JL）来获得OSE，且分析相对直接。

5. 链化方法

链化是比朴素网论证更精细的技术。对于稀疏度s=1的情况，直接对指数级多的向量进行联合界要求m非常大。Clarkson和Woodruff观察到，子空间中的向量可以分解为“集中”向量和“分散”向量。分散向量即使s=1也有很好的尾界，而集中向量的数量很少。结合链化技术，他们得到了 m ≈ d² * polylog(d/ε) / ε² 且 s=1 的结果。

子空间嵌入在回归中的应用 📈

上一节我们介绍了获取OSE的各种方法，本节我们来看看如何利用OSE来解决回归问题。除了上节课提到的“素描后求解”方法，还有两种更高效的使用方式。

方法二：基于常数OSE的梯度下降加速

梯度下降是一种迭代优化方法。对于最小二乘问题，目标函数是 f(x) = |Ax - b|₂²。梯度下降的更新公式为：
[ x_{t+1} = x_t - \eta \nabla f(x_t) ]
其中，梯度 ∇f(x) = 2A^T(Ax - b)。

定理： 设κ(A)为A的条件数（最大奇异值/最小奇异值）。经过 L = O(κ(A) log(1/ε)) 次迭代后，误差 |A x_L - A x^*|₂² 可以降至初始误差的ε倍以内。

问题在于，κ(A)可能很大。解决方案是使用一个常数（如1/6）子空间嵌入Π。计算ΠA的SVD：ΠA = UΣV^T。令 R = VΣ^{-1}。由于Π是子空间嵌入，对于所有x，有：
[ |x|₂ ≈ |ΠARx|₂ ≈ |ARx|₂ ]
因此，AR是一个条件数为常数的良态矩阵。我们在AR上运行梯度下降，迭代次数仅为 O(log(1/ε))。

计算复杂度：

计算ΠA的SVD：O(m d²) ≈ O(d³) （若m≈d）
每次迭代计算AR和(AR)^T：O(nnz(A) + d²)
总复杂度：O(nnz(A) log(1/ε) + d³)

为了获得最终(1+ε)倍的近似解，可以先使用方法一（素描求解）获得一个常数近似解作为初始点x0，然后使用此方法进行精炼。

方法三：改进的素描求解分析

此方法表明，即使Π不是ε-OSE，只要满足以下两个稍弱的性质，素描求解方法仍然有效：

常数子空间嵌入（例如，因子为1/6）。
具有Frobenius误差的近似矩阵乘法，其中误差参数ε' = √(ε/d)。

证明思路：
设x是原问题最优解，x̃是素描问题最优解。目标证明 |A x̃ - b|₂² ≤ (1+O(ε)) |A x - b|₂²。
将误差分解为：
[ |A x̃ - b|₂² = |A(x̃ - x)|₂² + |A x - b|₂² ]
利用正交性，只需证明第一项 ≤ ε |A x* - b|₂²。
通过引入A的SVD和常数嵌入性质，可以将问题转化为评估项 |U^T Π^T Π w|₂²。这正好是近似矩阵乘法的形式，利用性质2可得其界为 (ε/d) * d * opt² = ε * opt²。

优势：
此分析允许使用更少的行数m。例如，对于稀疏JL（s=1），使用方法一需要 m = O(d²/ε²)，而使用方法三只需要 m = O(d² + d/ε)。

总结 🎯

本节课我们一起深入探讨了无意识子空间嵌入。

我们回顾了五种获取OSE的主要方法：网论证、非交换Khintchine不等式、矩方法、近似矩阵乘法以及链化方法，了解了它们的基本思路和相对优劣。
在回归应用方面，我们学习了三种利用OSE的策略：
1. 直接的ε-OSE用于“素描后求解”。
2. 常数OSE用于预处理矩阵，从而加速梯度下降，将迭代次数依赖从多项式级(1/ε)降低到对数级log(1/ε)。
3. 结合常数OSE和特定精度的AMM性质，为“素描后求解”提供了更宽松且高效的分析框架，有时能用更少的测量数达到目的。

这些工具和分析技巧是构建高效大数据算法的核心。在下节课中，我们将看到类似的证明思路如何应用于低秩近似问题，并随后进入压缩感知的主题。

18：大数据算法中的压缩感知 🎯

在本节课中，我们将学习如何利用子空间嵌入进行低秩矩阵近似，并初步介绍压缩感知的基本概念。我们将看到，这些主题虽然看似不同，但在核心思想上紧密相连。

低秩近似与子空间嵌入 📉

上一节我们介绍了利用子空间嵌入解决回归问题的三种方法。本节中，我们来看看如何将类似的思想应用于低秩矩阵近似问题。

低秩近似的动机

低秩近近似的基本思想是，我们有一个巨大的矩阵 A（例如，用户对电影的评分矩阵），其维度为 n × d。我们假设矩阵的行（例如，用户）实际上存在于一个低维空间中。这意味着，尽管矩阵有 d 个显式属性（例如，电影评分），但只有少数几个隐藏属性真正决定了所有数据。

用公式表示，我们希望找到一个秩至多为 k 的矩阵 B，使得 A ≈ B。更具体地说，我们希望找到矩阵 U（n × k）和 V（k × d），使得 A ≈ U V。

根据Eckart-Young定理，矩阵 A 的最佳秩 k 近似可以通过其奇异值分解（SVD）获得。设 A = U Σ V^T，则最佳近似 A_k = U_k Σ_k V_k^T，其中下标 k 表示仅保留前 k 个奇异值和对应的奇异向量。

然而，计算完整的SVD需要 O(n d^2) 时间，对于大规模矩阵来说代价高昂。我们的目标是利用随机化技术来加速这一过程。

基于子空间嵌入的快速算法

以下定理为我们提供了加速的可能性。

定理：定义矩阵 Ã_k = Proj_{AΠ^T, k} (A)，即先将矩阵 A 的列投影到 AΠ^T 的列空间，然后取该投影矩阵的最佳秩 k 近似。那么，只要投影矩阵 Π 满足以下两个条件，就有：
||A - Ã_k||_F ≤ (1 + O(ε)) * ||A - A_k||_F

Π 是矩阵 V_k（A 的右奇异向量矩阵的前 k 列）列空间的 (1/2)-子空间嵌入。
Π 满足关于矩阵 V_k^T 和 A - A_k 的 √(ε/k)-近似矩阵乘法保证（在Frobenius范数下）。

这个定理的证明思路与回归问题的“方法三”非常相似。核心在于，我们可以将低秩近似问题中的误差项，转化为一个类似于回归问题的形式，并利用之前为回归问题建立的分析框架进行误差控制。

算法效率分析

现在，我们来分析为什么基于此定理的算法比直接计算SVD更快。

首先，为了满足定理的条件，我们需要选择 Π 的维度 m。条件1要求 m = O(k)，条件2要求 m = O(k/ε)。因此，m 可以取为 O(k/ε)。

以下是算法的关键步骤及其时间复杂度：

计算 AΠ^T：利用快速子空间嵌入（如稀疏嵌入或基于Hadamard/FFT的方法），可以在接近 O(nnz(A)) 的时间内完成，其中 nnz(A) 是 A 的非零元数量。
计算 AΠ^T 的SVD：设 AΠ^T = U‘ Σ’ V‘^T。由于 AΠ^T 是 n × m 矩阵，计算其SVD的时间为 O(n m^2) = O(n k^2 / ε^2)。
近似求解回归问题：我们不需要精确计算 U‘^T A。相反，我们可以求解一个广义回归问题：min_{X, rank(X)≤k} ||U‘ X - A||_F。利用上一讲关于回归的快速算法，我们可以得到一个近似解 X̃，满足 ||U‘ X̃ - A||F^2 ≤ (1+ε) * min ||U‘ X - A||_F^2。这一步的时间复杂度可以控制在 O((nnz(A) + poly(k/ε)) * log(1/ε)) 量级。
输出分解：计算 X̃ 的SVD，得到 X̃ = U_x Σ_x V_x^T。最终输出的低秩近似分解为 (U‘ U_x) Σ_x V_x^T。计算 X̃ 的SVD需要 O(m d^2) = O((k/ε) d^2) 时间，但通常 d 很大，我们可以利用 X̃ 是通过回归得到的这一事实，有时可以避免显式计算其SVD。

总体而言，当 k/ε 远小于 d 时，该算法的时间复杂度 O(n k^2 / ε^2 + (k/ε) d^2) 优于直接计算SVD的 O(n d^2)。通过进一步的优化（如近似求解回归），甚至可以将依赖 d 的部分从 d^2 降低到 d * poly(k/ε)。

与K-Means聚类的联系

低秩近似的思想与其他问题密切相关，例如K-Means聚类。可以将数据点矩阵 A（每行是一个点）的K-Means目标函数，重新表述为寻找一个特定形式的秩 k 近似矩阵的问题。这种联系使得我们可以利用降维技术（如Johnson-Lindenstrauss变换）来加速K-Means算法。研究表明，将数据维度降至 O(k/ε^2)，然后求解降维后的K-Means问题，可以得到原始问题 (1+ε) 近似的解。

压缩感知初步 🔍

接下来，我们将转向一个新的主题——压缩感知，它将占据接下来几讲的内容。压缩感知的核心思想是，利用信号的“可压缩性”，以远少于信号维度（奈奎斯特采样定理要求）的测量次数来恢复信号。

基本问题设定

假设我们有一个高维信号向量 x ∈ R^n，它是“可压缩的”。可压缩性通常意味着 x 是近似稀疏的。我们固定一个参数 k，定义 x_{tail(k)} 为将 x 中幅度最大的 k 个分量（头部）置零后得到的向量。如果 ||x_{tail(k)}|| 很小，我们就说 x 是高度可压缩的。

压缩感知的目标是，通过少数 m 个线性测量值来近似恢复 x。测量过程表示为 y = Π x，其中 Π 是一个 m × n 的测量矩阵（m << n）。我们希望找到一个高效算法，能够从 y 中计算出一个估计值 x̃，使得恢复误差 ||x - x̃|| 以 x_{tail(k)} 的范数为上界。

为什么关注稀疏性？

信号在某个变换域（而非标准基）下稀疏的现象非常普遍。一个经典的例子是自然图像。

图像在像素域（标准基）下并不稀疏，但在某些变换基（如小波基）下是近似稀疏的。以小波变换为例，它通过递归地计算图像局部块的平均值和差值，将图像能量集中到少数系数上。在得到的变换系数中，大部分系数（对应于平滑区域）接近于零，只有少数系数（对应于边缘和纹理）具有较大值。这就是JPEG等图像压缩标准能够有效工作的基础。

在压缩感知框架下，如果 x = Ψ α，其中 Ψ 是稀疏变换基（如小波基），α 是稀疏或近似稀疏的系数向量，那么我们的测量就变为 y = Π x = (Π Ψ) α。此时，问题转化为从测量值 y 中恢复稀疏向量 α。

主要理论结果

一个里程碑式的结论由Candès、Romberg、Tao以及Donoho等人给出：

定理：存在 m × n 的测量矩阵 Π（例如，随机高斯矩阵或随机符号矩阵），其中 m = O(k log(n/k))，以及一个多项式时间恢复算法。该算法以 y = Π x 为输入，输出 x̃，并保证以下误差界：
||x - x̃||2 ≤ (C / √k) * ||x||_1
其中 C 是一个常数。

这个定理意味着，所需的测量次数 m 主要取决于信号的稀疏度 k，并带有一个对数因子 log(n/k)。如果信号是严格 k-稀疏的，则恢复是精确的。如果是近似稀疏的，则恢复误差由信号“尾部”的能量控制。

值得注意的是，如果事先知道信号是严格 k-稀疏的，并且不考虑噪声，那么理论上 2k 次测量就是必要且充分的，这可以通过诸如Prony方法等更古老的技巧来实现。然而，为了对近似稀疏性和测量噪声具有鲁棒性，引入额外的对数因子是必要的。

总结 📝

本节课中，我们一起学习了两个重要主题。

首先，我们探讨了如何将子空间嵌入技术应用于低秩矩阵近似问题。通过构造满足特定条件的随机投影矩阵，我们可以显著加速寻找最佳低秩近似的进程，其核心思想与之前解决的回归问题一脉相承。

其次，我们初步介绍了压缩感知的基本框架。它旨在利用信号在某个域中的稀疏性或可压缩性，以远低于信号维度的测量次数来恢复信号。我们了解了其问题设定、实际应用背景（如图像处理），以及一个关键的理论结果，该结果指出了所需测量次数与信号稀疏度之间的基本关系。

从下节课开始，我们将深入压缩感知的理论核心，证明上述主要定理，并探索更多相关的算法与概念。

19：压缩感知与恢复算法 🧠

在本节课中，我们将学习压缩感知的核心定理，并探讨如何通过一个称为“基追踪”的算法，从少量测量中高效地恢复稀疏或近似稀疏的信号。我们将看到，这一过程与之前学习的约翰逊-林登斯特劳斯引理密切相关。

从压缩感知定理开始

上一节我们开始了压缩感知这一新章节。我陈述了由Candès、Romberg、Tao以及Donoho提出的定理，该定理指出：

存在一个矩阵和一个在多项式时间内运行的算法，使得对于任意向量 x（可将其视为 k-稀疏向量），如果你仅向算法提供 πx，它就能高效地恢复出一个接近 x 的向量 x̃。

具体来说，如果 x 确实是 k-稀疏的（即尾部范数为0），那么你将精确地恢复出 x。我们将在习题中看到，对于精确稀疏的情况，你可以仅用 2k 次测量，这既是必要的也是充分的。

然而，如果你希望对近似稀疏的向量实现这种近似恢复（即恢复误差取决于向量与真正稀疏状态的距离），那么你确实需要更多的测量次数。实际上存在一个下界，表明要实现这种保证，测量次数确实需要达到这个数量级（相差一个常数因子）。

证明定理与算法设计

今天我想要证明这个定理。你将看到，这个 π 矩阵与约翰逊-林登斯特劳斯引理密切相关。事实上，一旦我们完成设定，JL引理将意味着这样的 π 矩阵存在。

在本节课末，我还将展示一个逆命题：如果一个 π 矩阵满足我们在此所需的条件，那么你可以将其转化为一个JL映射。它们在某种意义上几乎是等价的，尽管在另一个方向上会存在一些损失。

自然的解码算法及其困难

让我们暂时忘记 π。假设 π 具有某种神奇的特性，允许我们进行解压缩。那么，一个自然的解码算法是什么呢？

一个自然的解码算法是：假设 x 确实是 k-稀疏的，那么存在唯一的 k-稀疏解。也就是说，存在唯一的 k-稀疏向量能产生我们得到的测量值 πx。

因此，一个自然的算法将是：
最小化向量 z 的支撑集大小（即 l₀ “范数”），约束条件为 πz = πx。

算法知道 πx，因此我们可以写出这个优化问题。不幸的是，一般来说，这是一个NP难问题。这可以追溯到Gary Johnson的NP难问题目录。这是一个计算上困难的问题，但它是针对任意输入求解时的困难。而我们并不需要解决任意输入的问题。

松弛到可处理的问题：基追踪

因此，我们将把问题松弛为一个非NP难的问题，即一个可处理的问题。我们将使用一个称为“基追踪”的算法。

基追踪算法表述为：
最小化向量 z 的 l₁ 范数，约束条件为 πz = πx。

这个优化问题可以写成一个线性规划问题，而线性规划可以在多项式时间内求解。这里的变量是向量 z 的各个条目。

那么问题就变成了：这种松弛何时能给出一个好的解？ 例如，如果 x 确实是 k-稀疏的，求解这个何时能精确地给出 x？或者如果 x 不是 k-稀疏但接近稀疏，何时能给出那种保证？

基本上，我们需要让这个方法奏效。我们必须说明，如果 π 具有某些良好的性质，那么求解这个松弛问题就会成功。

矩阵 π 所需的关键性质

“良好”意味着什么？如果 x 确实是精确 k-稀疏的，那么我们至少需要能够恢复任何 k-稀疏向量 x。

注意，我们需要的条件是：对于所有不同的 k-稀疏向量 x 和 y，有 πx ≠ πy。这等价于说，对于所有 2k-稀疏向量，π(x - y) ≠ 0。换句话说，没有 2k-稀疏向量位于 π 的零空间中。因为如果发生这种情况，我们就无法恢复了。

更形式化地说，令 I_S 为一个矩阵，其列是集合 S 对应的标准基向量。我们希望对于所有大小为 2k 的集合 S，矩阵 π I_S 是列满秩的。π I_S 就是取出 π 中属于集合 S 的那些列所构成的矩阵。我们希望说，如果你查看 π 的任何 2k 列，它们都是满秩的。如果不是，那么存在某个列集不是满秩的，就意味着零空间中存在我们不希望出现的向量。

因此，我们需要 (π I_S)ᵀ (π I_S) 是可逆的（没有零奇异值）。

更强的保证与限制等距性质

事实证明，如习题所示，存在一种构造此类 π 矩阵的方法，可以让你高效地（多项式时间内）恢复任何 k-稀疏向量。这是来自1700年代的Prony方法。

然而，如果我们想要更强的保证（即你总是拥有这种误差界），我们将要求 π 不仅非奇异，而且是良态的。我们将证明，如果 π 的每个 2k 列子集都是良态的，那么通过求解基追踪，你实际上可以获得那种保证。

定义：矩阵 π 满足 (ε, k)-限制等距性质，通常称为RIP，如果对于所有 k-稀疏向量 x，有：
(1 - ε) ||x||₂² ≤ ||πx||₂² ≤ (1 + ε) ||x||₂²

这看起来很像约翰逊-林登斯特劳斯引理，它基本上是作用于所有 k-稀疏向量集合的JL性质。这也等价于：对于所有大小为 k 的集合 S，π 中对应于 S（即 x 的支撑集）的那些列应该是良态的。这意味着所有奇异值都在 (1 - ε) 和 (1 + ε) 之间。

我们之前写过类似的东西：|| (π I_S)ᵀ (π I_S) - I || ≤ ε。这回到了子空间嵌入领域。基本上，我们希望 π 是 (ε, k)-子空间嵌入，作用于 C(n, k) 个不同的 k 维子空间。

如果你有一个高斯矩阵或随机符号矩阵，要使一个 k 维子空间以失败概率 δ 工作，你需要行数 m 至少为 (k + log(1/δ)) / ε²。如果你希望同时对所有这些不同的子空间都成立，设 δ < 1 / C(n, k)，然后通过并集界限，你将同时对所有子空间都成立。这意味着 m 至少为 (k + k log(n/k)) / ε²。k log(n/k) 是主导项，这就是所需行数来源的原因。

基本上，那些论文（Candès等人和Donoho）表明，如果你满足某个固定常数 ε（例如0.1）的 (ε, O(k))-RIP，那么如果算法是基追踪，你就能得到这种保证。

我们将在课堂上证明这一点。实际上，在那篇论文之后，还有其他论文加强了常数。今天我在课堂上要展示的是Candès后来的一个版本，它指出：为了得到那里的 l₂/l₁ 保证，只要 π 满足 (√2 - 1, 2k)-RIP 就足够了。这个 2k 显然是最优的，原因与我们之前说过的相同：你至少需要对 2k-稀疏向量非奇异。

确定性构造与公开问题

关于RIP的定义，请注意与JL不同，这里没有随机性。它是一个固定对象。如果你给我 ε 和 k，并要求给出一个 (ε, k)-RIP矩阵，你可能希望存在一个确定性构造，就像非相干矩阵和里德-所罗门码一样，没有随机性。

那么，是否存在RIP矩阵的显式构造呢？这是一个随机事物，随机符号矩阵以高概率是RIP的。然后你就使用那个矩阵，希望它是RIP的。你实际上并不知道你得到的是否是RIP。检查它是否RIP的最朴素方法就是检查这个条件，但你必须检查 C(n, k) 个不同的算子范数，这是一个指数级的数量，因此高效地检查RIP是不可行的。事实上，在某些情况下，检查RIP是NP难的。

如今，你必须这样做：如果你想要那么多行，只需随机选择一个矩阵，然后祈祷它实际上是RIP的。

那么显式RIP呢？我们将证明，如果你有一组向量，假设 π 是 (ε/k)-非相干的，那么它意味着 (ε, k)-RIP。如果你还记得我们需要多少行来获得 ε-非相干性，大约是 (1/ε²) log n。这意味着，使用里德-所罗门码（来自某次习题集），我们大致需要 O(k²/ε² * (log n / log log n + log k/ε)) 行。注意是 k²，这并不理想。它比 k 差。

实际上，这是一个重大的公开问题：是否存在一个确定性构造，其行数类似于 k * polylog(n)？目前无人知晓如何做到。最好的确定性构造是由Bourgain等人在2011年提出的。它说，对于 k 大约为 n^0.4999，可以得到行数最多为 k^1.999 * polylog(n) 的确定性构造。这也不容易，它使用了解析数论。

在组合数学和计算机科学中，这种情况经常发生：存在某个对象，通过概率方法可知其存在。如果你从某个均匀分布中随机取一个对象，它以高概率拥有你想要的属性，比如扩展图。但如果你要求我明确列出一个具有此属性的对象族，没人知道怎么做。我认为这就是我能给出的最好答案：没人知道。

非相干性蕴含RIP的证明

在证明那个界限之前，我想先证明另一个定理：Gershgorin圆盘定理。

该定理说，对于一个矩阵 A（元素为 a_ij），A 的所有特征值都位于以每个对角线元素 a_ii 为中心，半径为 ∑_{j≠i} |a_ij| 的复圆盘的并集中。对于我们而言，我们只处理具有实特征值的情况，所以它们就在这些区间内。

简单来说，如果你看 A 作为一个方阵，查看其对角线元素，每个特征值大致是这个值加上或减去该行其他元素的 l₁ 范数。所有特征值位于所有这些区间的并集中。

证明非常简单，基本上遵循特征值的定义。假设 x 是特征向量，特征值为 λ，那么 Ax = λx 意味着 λ x_i = ∑_j a_ij x_j。你可以由此推导出界限。如果需要，可以查阅在线资料。

现在，定理：(ε/k)-非相干性蕴含 (ε, k)-RIP。

我们关心什么特征值？我们关心矩阵 (π I_S)ᵀ (π I_S) 的特征值，其中 S 的大小最多为 k。π 的列是这些非相干向量，它们范数为1，且内积最多为 ε。

那么这个矩阵看起来像什么？对角线元素是1，非对角线元素基本上是 ±α，其中 α 大约是 ε/k。所有特征值都是 1 ± 该行所有其他元素的和。该行其他元素每个最多为 ε/k，且该行其他元素的数量是 k-1。因此，实际上 (ε/(k-1))-非相干性就蕴含 (ε, k)-RIP。这由Gershgorin圆盘定理可得。

证明基追踪的误差界

现在我想证明那个界限。为了得到 l₂/l₁ 保证，只要具有那种程度的RIP就足够了。

我将定义 x̃ = x + h，所以 h 定义为 x̃ - x。注意，基追踪（写在黑板上）是最小化 x̃ 的 l₁ 范数，即最小化 x + h 的 l₁ 范数，约束条件为 π x̃ = π x。这意味着 π (x + h) = π x，所以 h 在 π 的零空间中。

因此，最终我们需要的是一个矩阵 π，使得其零空间中的任何向量都是“平滑的”，即其 l₂ 范数最多是 (1/√k) 乘以其 l₁ 范数。为什么称之为平滑？因为如果一个向量全部集中在一个点上（比如向量 e1），那么它的 l₂ 范数和 l₁ 范数是相同的。但如果你想要在 l₂ 和 l₁ 范数之间有一个大的差距，你必须将质量分散开。

让我们证明这个定理。我将建立一些符号。

定义 h。令 T₀ 是 x 中幅度最大的 k 个坐标的集合（按绝对值）。令 T₁ 是 h 在 T₀ 补集中最大的 k 个坐标的集合。也就是说，查看 h，将 T₀ 中的元素置零，然后保留其余部分，T₁ 是这个向量中最大的 k 个坐标。以此类推，T_j 是 h 在 T₀ ∪ T₁ ∪ ... ∪ T_{j-1} 补集中最大的 k 个坐标。

T₀ 是特殊的，它是 x 中最大的 k 个坐标。然后，你不断剥离 h 中下一个最大的 k 个坐标，但限制在 T₀ 的补集中。

我们想要证明的是：||h||₂ ≤ O(1/√k) * ||x_{T₀^c}||₁。T₀ 是幅度最大的 k 个坐标，所以 T₀^c 是尾部。

注意，||h||₂ ≤ ||h_{T₀ ∪ T₁}||₂ + ||h_{(T₀ ∪ T₁)^c}||₂。我们将分别界定这两项。一旦我们证明了这两项，我们就完成了，因为根据三角不等式。

我们首先需要证明一个引理来处理第二项。

引理：∑_{j≥2} ||h_{T_j}||₂ ≤ (1/√k) * ||h_{T₀^c}||₁。

一旦有了这个引理，注意 h_{(T₀ ∪ T₁)^c} 这个向量正好是 h_{T₂} + h_{T₃} + h_{T₄} + ...。因此，它的范数由这个和界定，而根据引理，这正是我们想要的界定方式。所以，一旦我们证明了这个引理，第一项就由三角不等式完成了。

现在证明这个引理。注意，h_{T_j} 是一个 k-稀疏向量，因为我们投影到一个大小为 k 的集合上。对于一个 k-稀疏向量，其 l₂ 范数最多是 √k 乘以其 l_∞ 范数。最坏情况是所有 k 个条目都和最大的一样大，那么你就得到 √k 乘以 l_∞ 范数。

所以，∑_{j≥2} ||h_{T_j}||₂ ≤ √k * ∑_{j≥2} ||h_{T_j}||_∞。

现在注意，h_{T_j} 中的最大条目小于 h_{T_{j-1}} 中的每个条目，因为 j ≥ 2。T₂ 是什么？T₂ 是到目前为止已移除部分中剩余的最大坐标。但在 T₁ 中，我们已经抽出了 h 中最大的坐标。实际上，更具体地说，T₂ 中没有任何坐标能大于 T₁ 中的任何坐标，否则我们就会把它放在 T₁ 中了。

因此，||h_{T_j}||_∞ ≤ (1/k) * ||h_{T_{j-1}}||₁。因为 T_j 中的最大项小于 T_{j-1} 中的平均项，而平均项是 l₁ 范数除以 k。

所以，这最多是 (1/√k) * ∑_{j≥2} ||h_{T_{j-1}}||₁。而这个和正好是 ||h_{T₀^c}||₁。

现在，我将使用基追踪这一事实。基追踪告诉我，在所有满足 πh = 0 的向量 h 中，这个 h 是使 ||x + h||₁ 最小化的那个。这就是基追踪所做的：它最小化 l₁ 范数。

特别地，根据 l₁ 最小化，我们有 ||x||₁ ≥ ||x + h||₁。因为 x 本身是基追踪的一个可行解，但基追踪选择了 x + h 而不是 x，所以它的范数必须更小。

将 ||x + h||₁ 写成 ||(x + h)_{T₀}||₁ + ||(x + h)_{T₀^c}||₁。根据三角不等式，这至少是 ||x_{T₀}||₁ - ||h_{T₀}||₁ + ||h_{T₀^c}||₁ - ||x_{T₀^c}||₁。

现在重新排列这个不等式。它意味着 ||h_{T₀^c}||₁ ≤ ||x_{T₀^c}||₁ + ||x_{T₀^c}||₁ + ||h_{T₀}||₁ = 2||x_{T₀^c}||₁ + ||h_{T₀}||₁。

现在，h_{T₀} 是一个 k-稀疏向量，所以根据柯西-施瓦茨不等式，其 l₁ 范数最多是 √k 乘以其 `l

20：RIP蕴含JL与迭代硬阈值算法 🧮

在本节课中，我们将要学习两个核心内容。首先，我们将完成RIP蕴含分布JL引理的证明。其次，我们将介绍一种比基追踪更快的信号恢复算法——迭代硬阈值算法，并分析其理论保证。

从RIP到分布JL引理的证明 🔗

上一节我们介绍了Krahmer-Ward定理，它指出：如果一个矩阵 Π 满足 ε-RIP 性质，那么矩阵 ΠDσ（其中 Dσ 是一个对角线为随机符号的对角矩阵）将满足分布JL引理。本节中，我们将完成这个定理的证明。

证明的核心技巧是将目标表达式写成一个关于随机变量 σ 的二次型，然后利用Hanson-Wright不等式和Khintchine不等式等工具证明其集中性。

符号定义与矩阵分块

首先，我们定义一些符号。设 x 是我们希望保留的任意单位向量（|x|₂ = 1）。我们按 x 各分量的绝对值大小对其进行降序排列，这等价于对矩阵 Π 的列进行置换，不影响任何范数。

我们将排序后的 x 分成大小为 K 的块：

xᵢ 表示第 i 个块对应的 K-稀疏向量。
Πᵢ 表示 Π 中与第 i 个块对应的 K 列。

我们的目标是分析表达式 |ΠDσ x|₂²。我们可以将其写为二次型：
|ΠDσ x|₂² = σᵀ X σ，其中 X = Dₓ Πᵀ Π Dₓ。

现在，我们将矩阵 X 按 K×K 的块进行划分：

A：所有对角线块，其余部分为零。
B：第一行（除第一个对角线块外）的所有块，其余部分为零。
C：所有既不在对角线块，也不在 B 或 Bᵀ 中的块。

因此，X = A + B + Bᵀ + C。

我们需要证明，以至少 1 - 2^{-c‘K} 的概率，以下各项同时成立：

σᵀ A σ = 1 ± O(ε)
|σᵀ B σ|, |σᵀ Bᵀ σ|, |σᵀ C σ| ≤ O(ε)

分析对角线部分 A

对于 A，我们有：
σᵀ A σ = Σᵢ (Πᵢ D_{σᵢ} xᵢ)ᵀ (Πᵢ D_{σᵢ} xᵢ) = Σᵢ |Πᵢ D_{σᵢ} xᵢ|₂²

由于 xᵢ 是 K-稀疏的，且 Π 满足 ε-RIP，因此对于每个块 i，有：
|Πᵢ D_{σᵢ} xᵢ|₂² = (1 ± ε) |D_{σᵢ} xᵢ|₂² = (1 ± ε) |xᵢ|₂²

求和后得到：
σᵀ A σ = (1 ± ε) Σᵢ |xᵢ|₂² = (1 ± ε) |x|₂² = 1 ± ε

分析非对角线部分 B 和 Bᵀ

我们以 B 为例。B 的结构使得 σᵀ B σ 可以写为：
σᵀ B σ = σ₁ᵀ (D_{x₁} Π₁ᵀ Π₋₁ D_{x₋₁}) σ₋₁
其中 σ₁ 和 σ₋₁ 是独立的。

关键观察是，向量 v = D_{x₁} Π₁ᵀ Π₋₁ D_{x₋₁} σ₋₁ 的 L₂ 范数很小。我们可以利用 Π 的 RIP 性质和一种称为“壳层分解”的技巧来证明 |v|₂ = O(ε/√K)。

一旦我们固定了 v，表达式 σ₁ᵀ v 就是一个固定向量与随机符号向量的点积。根据 Khintchine不等式，这个值以高概率很小。具体地，其超过 O(ε) 的概率上界为 exp(-Ω(ε² K / ε²)) = exp(-Ω(K))。

Bᵀ 的分析完全对称。

分析剩余部分 C

对于 C，我们直接应用 Hanson-Wright不等式。由于 C 的期望为 0，该不等式表明：
P(|σᵀ C σ| > λ) ≤ exp(-c min(λ²/|C|_F², λ/|C|))。

我们需要估计 C 的弗罗贝尼乌斯范数 |C|_F 和算子范数 |C|。通过再次利用 RIP 性质（证明两个不相交的 K-稀疏向量经 Π 变换后的内积为 O(ε)）和壳层分解技巧，我们可以证明：

|C|_F² = O(ε² / K)
|C| = O(ε / K)

取 λ = Θ(ε)，代入Hanson-Wright不等式，两个指数项都变为 exp(-Ω(K))。

结论

将 A, B, Bᵀ, C 的分析结合起来，我们得到：
|ΠDσ x|₂² = σᵀ X σ = (1 ± O(ε)) ± O(ε)
即以极高的概率（1 - exp(-Ω(K))）有 |ΠDσ x|₂² = 1 ± O(ε)。这正好是分布JL引理的定义。证明完成。

迭代硬阈值算法 ⚡

上一节我们证明了RIP与JL引理之间的深刻联系。本节中，我们来看看如何利用类似的思想设计更高效的信号恢复算法。

基追踪（Basis Pursuit）虽然能在多项式时间内求解，但在高维（n很大）场景下，其计算成本（如 O(n³)）仍然过高。因此，我们需要更快的算法。

迭代硬阈值算法是一种高效的迭代方法，它从一个初始估计（如零向量）开始，通过不断迭代来改进估计值。

算法描述

输入：

观测向量 y
测量矩阵 Π
迭代次数 T
（假设 y = Πx + e，其中 x 近似稀疏，e 是噪声）

算法步骤：

初始化 x₀ = 0
对于 i = 1 到 T，执行：
a_{i} = x_{i-1} + Πᵀ (y - Π x_{i-1})
x_{i} = H_K(a_{i})
其中 H_K(z) 是硬阈值算子，它保留向量 z 中绝对值最大的 K 个分量，其余置零。

理论保证

定理：如果矩阵 Π 满足 ε-3K RIP，且 ε < 1/(4√2)，那么迭代硬阈值算法产生的估计 x_T 满足：
|x_T - x|₂ ≤ C₁ * 2^{-T} * |H_K(x)|₂ + C₂ * |x - H_K(x)|₂ + C₃ / √K * |x - H_K(x)|₁ + C₄ * |e|₂

解读：

第一项 2^{-T} |H_K(x)|₂：这是初始误差的主要部分，但随着迭代次数 T 增加而指数衰减。
第二项 |x - H_K(x)|₂：这是信号尾部的 L₂ 范数。
第三项 (1/√K) |x - H_K(x)|₁：这是信号尾部的 L₁ 范数除以 √K。通过“壳层分解”技巧，可以证明此项与基追踪的经典误差界 （1/√K） L₁-tail* 同阶。
第四项 |e|₂：测量噪声的 L₂ 范数，无法避免。

因此，经过足够多次迭代（T = O(log(1/ε))）后，算法的误差主要受限于信号本身的近似稀疏性和测量噪声，这与最优的基追踪恢复保证相匹配。

证明思路概要

证明的核心是跟踪残差 r_t = x - x_t 的演化。通过将更新步骤代入，并利用 Π 的 3K-RIP 性质（用于分析涉及最多 3K 个列的子矩阵），可以推导出残差范数之间的递归关系：
|r_{t+1}|₂ ≤ (1/2) * |r_t|₂ + C * (|e|₂ + |x - H_K(x)|₂ + (1/√K)|x - H_K(x)|₁)

通过归纳法，这个递归关系直接导出了定理中的误差上界。

总结 📚

本节课中我们一起学习了两个重要内容：

RIP蕴含分布JL引理：我们完成了Krahmer-Ward定理的证明，展示了如何利用矩阵的RIP性质，通过添加随机符号对角线矩阵，构造出满足分布JL引理的变换。这连接了压缩感知和降维中的两个核心概念。
迭代硬阈值算法：我们介绍了一种高效、迭代的信号恢复算法。该算法在测量矩阵满足RIP的条件下，能以线性收敛速度达到与基追踪相近的恢复精度，同时计算复杂度远低于求解线性规划。这为处理大规模稀疏恢复问题提供了实用的工具。

这些内容体现了大数据算法中“理论指导设计，高效实现理论保证”的核心思想。

21：压缩感知与01恢复 🎯

在本节课中，我们将学习压缩感知领域的最后一个核心主题：01恢复。我们将探讨如何利用一种称为稀疏图的结构，在保证快速计算的同时，实现最优的测量次数和L1-L1误差保证。

背景回顾与目标设定 📋

上一讲我们介绍了迭代硬阈值算法。只要我们的测量矩阵 Π 具有良好的RIP常数（例如小于0.1或某个固定常数），我们就能获得L2-L1误差保证。该算法每次迭代的时间取决于矩阵 Π 的矩阵-向量乘法速度。

如果我们能自由选择测量矩阵，通常会选择从离散傅里叶变换中采样行向量。这样，应用 Π 到任意输入向量大致是线性时间。然而，这种方法存在一个缺点：我们无法通过傅里叶采样获得最优的测量次数。已知的最佳边界大约是 k * log n * log² k，而不是像随机高斯矩阵那样的 k * log(n/k)。

在许多实际应用中，k 与 n 处于同一数量级（例如 k = n/100）。我们希望将向量压缩50倍，此时常数因子至关重要。当 k 与 n 相当时，log n 是一个较大的值，而 log(n/k) 则接近常数。因此，我们希望能获得 k * log(n/k) 量级的最优测量次数。

我们将看到，如果我们将误差保证从 L2-L1 放宽到 L1-L1，我们就能实现这一目标，同时保持快速的迭代时间。

L1-L1 误差保证 定义为：我们恢复出的向量 x̃ 满足：

||x - x̃||₁ ≤ C * ||x - x_k||₁

其中 x_k 是 x 的最佳k项近似（尾部）。这是一个比L2-L1更弱的保证，因此更容易实现。

RIP1 矩阵与扩展图 🌉

为了实现上述目标，我们将使用一种称为 RIP1 矩阵的工具。

RIP1 矩阵定义

矩阵 Π 满足 (k, ε)-RIP1，如果对于所有k稀疏向量 x，都有：

(1 - ε) ||x||₁ ≤ ||Πx||₁ ≤ (1 + ε) ||x||₁

为方便起见，我们通常假设 Π 不会过度放大或缩小L1范数。

扩展图

RIP1矩阵可以通过一种称为非平衡扩展图的结构来构造。

定义：一个二分图 G = (U, V, E) 是一个 (N, M, D, k, ε) 扩展器，如果：

|U| = n, |V| = m。
图是左D-正则的（U中每个顶点恰好连接到V中的D个顶点）。
对于任意大小不超过k的子集 S ⊆ U，其邻居集合 N(S) 的大小满足：
```
|N(S)| ≥ (1 - ε) * D * |S|
```

直观上，这意味着对于任意不大的顶点集合S，其在右侧的邻居数量几乎达到了最大值（D * |S|），碰撞损失不超过ε比例。

定理：存在这样的扩展图，其参数为：

M ≈ k * log(n/k) / ε² （这正是我们想要的测量次数）
D ≈ log(n/k) / ε²

该定理可以通过概率方法证明：随机生成一个左D-正则的二分图，并以高概率满足扩展性条件。

从扩展图到测量矩阵

我们定义测量矩阵 Π 为：

Π = (1/D) * A_G

其中 A_G 是扩展图G的邻接矩阵（一个 m×n 的矩阵，每列有D个1，表示该左顶点连接的右顶点）。

这个矩阵具有快速矩阵-向量乘法的特性。因为每列大约有 D ≈ log n 个非零元，所以应用 Π 或 Π^T 到一个向量的时间复杂度约为 O(n log n)。

定理：如果G是一个 (k, ε) 扩展器，那么按上述方式构造的 Π 满足 (k, 2ε)-RIP1。

恢复算法：SSMP 🚀

虽然存在基于扩展图直接分析的算法（如扩展图匹配追踪EMP和稀疏匹配追踪SP），但我们将介绍一个仅依赖于RIP1抽象性质的算法及其分析：顺序稀疏匹配追踪。

SSMP算法是一个迭代算法，其伪代码如下：

初始化: x^0 = 0
参数: 迭代次数 T
for j = 1 to T:
    x^{j,0} = x^{j-1}
    for a = 1 to (C-1)k:
        // 找到最优的1-稀疏向量 z 来添加到当前估计中
        z = argmin_{1-sparse z} ||b - Π(x^{j,a-1} + z)||₁
        x^{j,a} = x^{j,a-1} + z
    // 通过硬阈值化保持解的稀疏性
    x^j = H_k(x^{j, (C-1)k})

其中，b = Πx + e 是观测向量（含噪声e），H_k 是保留最大k个分量的硬阈值算子。

算法的核心思想是：在每轮外层迭代中，我们通过内层循环多次添加一个最优的1-稀疏向量，逐步改进当前估计，最后通过阈值化保持解的稀疏性。

算法分析概要

为了分析算法，我们假设信号 x 本身就是k稀疏的（通过将尾部误差归入噪声项，该假设不失一般性）。我们的目标是证明，每轮外层迭代都能将误差减少一个常数因子。

定义第j步的残差为 y^j = x - x^j。分析将分为四个关键步骤：

内层循环的渐进改进：证明每次内层迭代都能将残差的观测范数减少一个因子。
```
||b - Πx^{j,a}||₁ ≤ (1 - 1/(2k+a))^{1/2} * ||b - Πx^{j,a-1}||₁
```
内层循环的累积效果：经过 (C-1)k 次内层迭代后，观测残差显著减小。
```
||b - Πx^{j, (C-1)k}||₁ ≤ (1/8) * ||b - Πx^{j-1}||₁
```

从观测误差到真实误差：利用RIP1性质，将观测误差的减小转化为真实L1误差的减小。
```
||x - x^{j, (C-1)k}||₁ ≤ (1/4)||x - x^{j-1}||₁ + O(||e||₁)
```
硬阈值化的影响：证明最后的硬阈值化步骤不会显著增大误差。
```
||x - x^j||₁ ≤ (1/2)||x - x^{j-1}||₁ + O(||e||₁)
```

通过归纳法，上述不等式意味着经过T轮迭代后，最终误差为：

||x - x^T||₁ ≤ (1/2)^T * ||x||₁ + O(||e||₁)

步骤1的证明是整个分析中最巧妙的部分，它依赖于几个关于RIP1矩阵和向量对齐性质的引理，而不需要深入扩展图的具体结构。这正是该分析方法的优势所在——它仅依赖于RIP1这一抽象性质。

关于算法的效率：虽然内层循环需要寻找“最优的1-稀疏向量”，这看似需要O(n)时间，但可以通过维护优先队列等数据结构，将每次内层迭代的时间复杂度降低到O(log n)，从而使整个算法保持高效。

总结 📝

本节课我们一起学习了压缩感知中实现01恢复的高级方法。我们主要探讨了：

目标：在保证快速迭代（近线性时间）的前提下，获得最优测量次数 m = O(k log(n/k)) 和 L1-L1 误差保证。
工具：引入了 RIP1 矩阵和非平衡扩展图。扩展图可以构造出满足RIP1且支持快速运算的测量矩阵。
算法：介绍了顺序稀疏匹配追踪算法。该算法通过迭代地添加1-稀疏向量来逼近信号，并利用硬阈值化保持稀疏性。
分析核心：算法的分析（特别是关键的步骤1）仅依赖于测量矩阵的RIP1性质，而不依赖于扩展图的具体结构，体现了抽象分析的威力。

通过本讲内容，我们看到了如何通过巧妙的算法设计和分析，在压缩感知中同时优化测量次数、计算复杂度和恢复精度，这也是大数据算法研究的核心魅力所在。

22：矩阵补全 🧩

在本节课中，我们将要学习矩阵补全问题，有时也被称为Netflix问题。我们将探讨其核心思想、算法以及理论分析。

概述

矩阵补全问题的核心是：给定一个巨大矩阵 M 的一部分随机观测条目，我们如何恢复整个矩阵？我们假设矩阵 M 是低秩的，这意味着它可以分解为两个较小矩阵的乘积。这个假设在许多实际应用中（如用户-电影评分矩阵）是合理的，因为用户偏好通常可以由少数几个“用户类型”的线性组合来解释。

问题定义与算法

我们假设矩阵 M 的维度是 n × n（为简化讨论），并且其真实秩为 r。我们观测到的是 M 中一个随机子集 Ω 上的条目，观测数量为 m。

理想情况下，我们希望找到与观测值一致且秩最小的矩阵。然而，最小化秩是一个NP难问题。因此，类似于压缩感知中将L0范数松弛为L1范数，我们将矩阵的秩松弛为其核范数（即奇异值之和）。

核范数最小化 (NNM) 算法公式如下：

minimize ||X||_*
subject to R_Ω(X) = R_Ω(M)

其中，||X||_* 表示矩阵 X 的核范数，R_Ω 是一个线性算子，它将矩阵限制在索引集 Ω 上，其他位置置零。

这个优化问题可以表述为一个半定规划问题，因此存在多项式时间算法求解。

理论分析的关键条件

上一节我们介绍了NNM算法，本节中我们来看看保证该算法成功恢复原始矩阵 M 所需满足的关键条件。这些条件是确定性的，如果成立，则恢复成功。随后我们将说明，当观测数量 m 足够大时，这些条件以高概率成立。

以下是三个关键条件：

近似矩阵乘法条件：|| (n²/m) * P_T R_Ω P_T - (m/n²) P_T || < 1/2
这个条件确保了在子空间 T（与 M 的行列空间相关的子空间）上，随机采样算子 R_Ω 的行为类似于恒等算子的缩放版本。这本质上是近似矩阵乘法性质的应用。
算子范数边界条件：|| R_Ω || ≤ C * log n
这个条件表明采样矩阵 R_Ω 的算子范数（最大奇异值）不会太大，这可以通过切尔诺夫界等概率工具来证明。
对偶证书存在性条件：存在一个矩阵 Y 位于 R_Ω 的值域中，使得：
- || P_T(Y) - UV^T ||_F ≤ √r / (2n)
- || P_{T^⊥}(Y) || ≤ 1/2
  这个条件是为了构造一个“对偶证书”，用于证明 M 是核范数最小化问题的唯一解。其证明涉及更多层的近似矩阵乘法论证。

其他算法简介

除了核范数最小化，还存在其他更快的迭代算法用于矩阵补全，特别是在处理含噪声数据时。

以下是两种常见的迭代方法：

软阈值迭代算法：这种算法类似于基追踪去噪中的迭代方法。它通过梯度下降步骤更新当前矩阵估计，然后对奇异值进行软阈值操作以促进低秩性。其更新步骤涉及计算SVD，但可以尝试使用本课程中学到的子空间嵌入技术来加速低秩近似步骤。
交替最小化算法：该算法假设 M ≈ XY^T，并交替地优化 X 和 Y。每一步都是一个回归问题，同样可以考虑使用快速近似回归技术来加速。

总结

本节课中我们一起学习了矩阵补全问题。我们从问题定义出发，介绍了基于核范数最小化的主流算法。我们深入探讨了保证算法成功所需的理论条件，特别是近似矩阵乘法在分析中扮演的核心角色。最后，我们简要了解了其他更快速的迭代算法。矩阵补全是一个连接了低秩近似、压缩感知和随机采样理论的丰富领域。

23：外部存储模型 🗂️

在本节课中，我们将要学习一种新的计算模型——外部存储模型。这个模型也被称为磁盘访问模型，由Hder和Bitter提出。我们将探讨当数据量过大，无法全部装入内存（RAM）时，如何设计高效的算法和数据结构来处理存储在磁盘上的数据。

模型介绍与基本原理

上一节我们提到了大数据处理的不同场景，例如流式处理。本节中，我们来看看当数据可以存储但无法全部放入内存时的情况。

外部存储模型的基本思想是，系统包含一个容量无限的磁盘和一个容量有限的内存。内存被划分为多个“页”，每个页的大小为 B，可以容纳 B 个元素。磁盘上的数据也被划分为大小为 B 的块。

内存操作：如果所需数据已经在内存中，读写操作是免费的。
磁盘操作：如果数据不在内存中，则需要从磁盘读取相应的块到内存中。这个操作的成本记为 1次I/O。如果内存已满，还需要决定将哪个内存页的内容写回磁盘。

我们的目标是最小化总的I/O操作次数。

模型合理性

为什么这个模型是合理的？因为读写磁盘（尤其是随机读写）的速度远慢于读写内存。以下是基于传统机械硬盘（HDD）的一些近似数据对比：

顺序读写：速度可达约 140 MB/s。
随机读写：速度约为 140次 I/O/s（每次读写4KB数据块）。
内存（RAM）读写：速度可达 5-10 GB/s。

随机读写和顺序读写之间的速度差异巨大（超过六个数量级）。因此，模型假设每次I/O能读取一个大小为 B 的块是合理的，因为进行一次磁盘寻址（随机访问）的耗时，足以让我们顺带读取 B 个元素，从而分摊了寻址成本。

基础问题示例

在理解了模型之后，我们来看几个基础问题在这个模型下的表现。

顺序读取

读取 N 个连续存储的元素非常简单：进行一次线性扫描即可。所需的I/O次数为 O(N/B + 1)。

矩阵乘法 🧮

接下来，我们讨论一个稍微复杂的问题：矩阵乘法。这是本课程中唯一涉及线性代数的部分。

问题：给定两个 N x N 的矩阵 X 和 Y，计算它们的乘积 Z = X * Y。

假设：作为算法设计者，我们可以决定数据在磁盘上的存储格式。这里我们假设矩阵以块的形式存储。

算法思路：

将矩阵 X 和 Y 分别划分为大小为 √(M/2) x √(M/2) 的块。
结果矩阵 Z 的每个块，可以通过 X 和 Y 中对应行和列的块进行矩阵乘法再求和得到。例如，Z 的左上角块等于 X 的第一行块与 Y 的第一列块的点积之和。
计算时，一次将参与计算的一对块（来自 X 和 Y）读入内存。由于每个块的大小约为 M/2 个元素，一对块可以同时放入内存。
在内存中使用任何矩阵乘法算法进行计算。

I/O复杂度分析：

输出矩阵 Z 有 O((N/√M)^2) 个块。
对于 Z 的每个块，需要计算 O(N/√M) 对块的乘积。
读取每一对块需要 O(M/B) 次I/O。
因此，总I/O次数约为 O((N/√M)^3 * (M/B)) = O(N^3 / (B√M))。

注意：使用更高效的矩阵乘法算法（如Strassen算法）可以进一步改进这个指数。

外部存储模型下的数据结构

现在，让我们看看如何在外部存储模型中实现一些常见的数据结构。

链表 🔗

在内部内存中，链表支持 O(1) 时间的插入/删除和 O(K) 时间的遍历（访问连续 K 个元素）。在外部存储模型中，如果简单地将每个节点随机存储在磁盘块中，遍历 K 个元素可能需要 O(K) 次I/O，这非常低效。

优化思路：将连续的链表节点分组到大小为 B 的“组”中，并维护一个不变式：除了最后一个组，每个组至少半满（即包含至少 B/2 个节点）。

操作：

插入：找到对应组，插入节点。如果组已满（超过 B 个节点），则将其分裂成两个组。
删除：从组中删除节点。如果组变得少于 B/2 个节点，则将其与相邻组合并，必要时再进行分裂以保持组大小在 [B/2, B] 之间。
遍历：由于每个组至少半满，要访问 K 个元素，最多需要读取 O(1 + K/(B/2)) = O(1 + K/B) 个组（即I/O次数）。

这样，插入和删除的摊销I/O复杂度为 O(1)，而遍历的I/O复杂度优化为 O(1 + K/B)。

前驱查询与B树 🌳

前驱查询是一个经典的数据结构问题：维护一个包含 N 个可比较键值对的集合，支持查询给定键 K 的前驱（即集合中小于等于 K 的最大键）。

在内部内存的比较模型中，有序数组+二分查找可以达到 O(log N) 查询时间。在更强大的Word RAM模型中，存在更高效的算法（如vEB树）。

在外部存储模型中，对应的经典数据结构是 B树。

B树本质：B树是一种 (a, b)-树，其中每个节点（页）的大小对应于磁盘块大小 B。通常设置 a = B/2, b = B。这意味着：

每个内部节点至少有 B/2 个子节点（根节点除外）。
每个节点（页）可以一次性读入内存。

性能：由于树的高度为 O(log_B N)，且每次访问一个节点只需一次I/O，因此B树的查询、插入、删除操作的I/O复杂度均为 O(log_B N)。

思考：B树的插入复杂度是 O(log_B N) = O((log N) / (log B))。有没有可能获得接近日志文件（O(1/B) 摊销插入）的插入效率，同时保持近似B树的查询效率？

更高效的数据结构变体

答案是肯定的。以下是两种改进的数据结构。

缓冲存储树 (BRT)

BRT是对 (2,4)-树（即2-3-4树）的简单修改。

结构：

每个内部节点关联一个大小为 B 的缓冲区。
所有数据项要么存储在叶子节点（叶子节点大小在 [B/2, B] 之间），要么存储在某个内部节点的缓冲区中。

操作：

插入：新元素直接插入根节点的缓冲区。如果根缓冲区满，则将其中的元素“冲刷”到合适的子节点缓冲区。这个过程可能递归进行，直到元素到达叶子节点，可能引发叶子分裂。
查询：像普通B树一样搜索，但同时检查路径上每个节点的缓冲区。

性能：

查询：O(log N) 次I/O（树高）。
插入（摊销）：O((log N) / B) 次I/O。分析思路是，每次冲刷涉及 Θ(B) 个元素，但只产生常数次I/O，因此每个元素每次被冲刷只需分摊 O(1/B) 的成本。一个元素最多被冲刷树高次。
插入（最坏情况）：通过每次只向接收元素最多的一个子节点进行冲刷，可以限制为 O(log N) 次I/O。

缓冲B^ε树

这是BRT思想的一个泛化，使用 (B^ε, 2B^ε)-树作为底层结构，其中 ε 是一个常数（如0.5）。

性能：

查询：O((1/ε) * log_B N) 次I/O。
插入（摊销）：O((1/(εB^(1-ε))) * log_B N) 次I/O。
插入（最坏情况）：O((1/ε) * log_B N) 次I/O。

当 ε 为常数时，查询复杂度与B树同量级（仅常数因子更大），但摊销插入复杂度远优于B树（除以 B 的多项式级而非对数级）。

外部排序

排序是算法的基础问题。在外部存储模型中，最优的基于比较的排序算法是 (M/B)-路归并排序。

算法步骤：

如果数据量 N ≤ B，直接读入内存排序。
否则，将数据分成 M/B 个组。
递归地对每个组排序。
合并这 M/B 个有序组。合并时，我们可以在内存中为每个组维护一个当前数据块（大小 B），因此内存足以同时容纳所有组的当前块。

递归式：设 T(N) 为排序 N 个元素所需的I/O次数。

T(N) = O(N/B) + (M/B) * T(N / (M/B))，当 N > B
T(N) = O(1)，当 N ≤ B

解此递归式，得到 T(N) = Θ((N/B) * log_{M/B} (N/B))。

最优性证明（基于比较模型）：

考虑所有满足“每个大小为 B 的块内部已排序”的输入排列。这样的排列总数约为 N! / (B!)^(N/B)。
算法需要识别出具体是哪种排列，信息量约为 N log(N/B) 比特。
每次I/O最多能带来 B log(M/B) 比特的新信息（因为I/O引入的 B 个新元素与内存中已有的 M 个元素交织的可能方式最多为 (M+B choose B) 种）。
因此，至少需要 (N log(N/B)) / (B log(M/B)) = (N/B) * log_{M/B} (N/B) 次I/O。

这证明了该归并排序算法的I/O复杂度在比较模型下是最优的。

总结

本节课中，我们一起学习了外部存储模型（磁盘访问模型）。我们从模型的基本定义和合理性出发，探讨了在该模型下如何分析基础操作（如顺序读取）和经典问题（如矩阵乘法）。接着，我们深入研究了如何为链表和前驱查询设计高效的外部存储数据结构，重点介绍了B树及其两种更高效的变体——缓冲存储树和缓冲B^ε树。最后，我们分析了外部排序问题，并得出了基于比较的最优算法及其复杂度。

核心要点：

模型核心：最小化慢速磁盘的I/O次数，利用块访问（大小 B）来分摊成本。
数据结构设计关键：利用数据的局部性，将相关元素组织在同一磁盘块中。
性能衡量：操作复杂度通常表示为 N（数据量）、M（内存大小）、B（块大小）的函数。

在下节课中，我们将探讨一个更高级的话题：缓存无关算法。这些算法无需预先知道参数 M 和 B，就能自动达到接近知晓这些参数时的最优性能。

24：缓存无关模型 🧠

在本节课中，我们将要学习一种新的计算模型——缓存无关模型。我们将了解它与之前讨论的外部存储模型有何不同，并探索如何在该模型下设计高效的算法，即使算法本身对系统的内存和块大小参数一无所知。

缓存无关模型简介

上一节我们介绍了外部存储模型，其目标是最小化I/O操作次数。在该模型中，我们假设内存大小为M，并被划分为M/B个页面，每个页面可以容纳一个大小为B的块。

本节中，我们来看看缓存无关模型。它与外部存储模型相似，但有两个关键区别：

算法不允许知道系统的参数，如M或B。
算法不控制缓存替换策略。我们假设操作系统处理所有替换，并做出最优选择。

这意味着在我们的分析中，可以假设操作系统进行的驱逐操作至少和我们分析中假设的一样好。

为何需要缓存无关模型？🤔

以下是缓存无关模型具有实际意义的几个原因：

可移植性：算法无需针对不同机器的参数进行调优。
环境适应性：在多进程环境中，分配给算法的有效内存M可能随时间变化。缓存无关算法的性能能自动适应可用的最佳块大小和内存。
多级内存层次结构：现代计算机有磁盘、RAM、L3/L2/L1缓存等多级存储，每级都有自己的M和B参数。缓存无关算法的性能分析可以同时适用于任何相邻的两级存储层次。

关于操作系统“最优选择”的合理性

你可能会质疑：操作系统如何能做出最优的缓存替换选择？这似乎需要预知未来。

这个假设实际上与在线算法领域的一个经典问题——分页问题——密切相关。以下是其核心思想：

问题设定：内存可容纳K = M/B个页面。有一个页面访问请求序列。如果请求的页面在内存中，可以免费访问；否则，需要从磁盘载入（可能需驱逐一个旧页面），这产生一次I/O成本。
最优算法（OPT）：如果知道未来，最优策略是驱逐最晚再次被使用的页面。
实际算法：操作系统使用如LRU或FIFO等无需预知未来的策略。
理论保证：研究表明，LRU和FIFO等算法是K-竞争的。更重要的是，在资源增强分析下，它们对于只有K/2大小内存的“先知”最优算法是2-竞争的。

结论：这意味着，如果一个缓存无关算法在我们的分析（假设最优驱逐）下具有I/O复杂度T(N, M, B)，那么在实际使用LRU等策略的操作系统下，其I/O复杂度为O(T(N, M/2, B))。因此，“操作系统做出最优选择”这个假设在常数因子内是合理的。

缓存无关算法示例

现在，我们开始探讨缓存无关模型下的一些具体算法。

数组遍历与反转

我们之前已经见过一个简单的缓存无关算法：扫描数组。

遍历一个大小为N的连续数组，I/O复杂度为 O(1 + N/B)。
反转数组可以通过从左右两端同时向中间扫描并交换元素来实现，复杂度同样为 O(1 + N/B)。

矩阵乘法 🧮

在外部存储模型中，我们通过将矩阵划分为大小为 √M × √M 的子块来进行高效的矩阵乘法，但该算法需要知道M。

在缓存无关模型中，我们使用递归布局和递归分治算法：

递归布局：将一个 N × N 的矩阵递归地划分为四个 N/2 × N/2 的子矩阵（A11, A12, A21, A22），并将它们在内存中连续排列。

递归算法：计算 C = A × B 时，使用标准的分块矩阵乘法公式递归计算。

# 伪代码示意
def mat_mul(A, B):
    if size(A) <= sqrt(M):  # 基础情况：子矩阵能放入内存
        return naive_mul(A, B)
    else:
        C11 = mat_mul(A11, B11) + mat_mul(A12, B21)
        C12 = mat_mul(A11, B12) + mat_mul(A12, B22)
        C21 = mat_mul(A21, B11) + mat_mul(A22, B21)
        C22 = mat_mul(A21, B12) + mat_mul(A22, B22)
        return combine(C11, C12, C21, C22)

分析：通过分析递归树，当子矩阵大小约等于 √M 时，它们可以完全放入内存。最终得到的I/O复杂度为 O(N²/B + N³/(B√M))，这与外部存储模型中的最优结果一致，但算法本身并不知道M。

链表操作 ⛓️

链表需要支持插入、删除和遍历K个元素的操作。在外部存储模型中，我们通过将节点分组来优化，但这需要知道B。

缓存无关的链表设计思路如下：

存储：所有节点存储在一个大数组中。
插入(P, X)：在数组末尾写入X，并更新指针。I/O复杂度 O(1)。
删除(P)：仅将节点P标记为“已删除”。I/O复杂度 O(1)。
遍历(P, K)：
1. 从P开始，跟随指针读取K个元素。
2. 将这K个元素连续地重写到数组末尾。
3. 将它们原来的位置标记为“已删除”。
定期整理：每进行Θ(N)次操作后（遍历算作K次操作），将整个数据结构重写到一个新的连续数组中，释放旧空间。

摊销分析：通过银行家方法进行摊销分析，可以将非连续读取（“运行”）产生的额外I/O成本“收费”给导致该碎片化的前序操作。最终，摊销后的遍历复杂度为 O(1 + K/B)。

注：存在更复杂的数据结构可以实现最坏情况下的高效操作，但这超出了本简介的范围。

静态前驱查询与van Emde Boas布局 🌳

对于静态前驱查询（例如在有序集中查找比某值小的最大元素），外部存储模型中使用B树，这需要知道B。

缓存无关模型中使用van Emde Boas (vEB) 布局：

递归布局：构建一棵完美的二叉搜索树。递归地将其分割：树高为 h，则分割点设在 h/2 处。将顶部的子树（大小为 √N）和底部的 √N 棵子树分别递归地进行连续内存布局。
查询：像在普通二叉搜索树中一样，从根节点开始向下遍历路径。
分析：考虑递归树中，子树大小首次小于等于B但大于等于 √B 的那一层。在这一层，每棵子树在内存中是连续的，因此最多需要2次I/O即可读入。查询路径需要访问 O(log N) 个这样的子树，因此总I/O复杂度为 O(log_B N)，达到了B树的效果而无需知道B。

排序：懒惰漏斗排序 🌀

排序在外部存储模型中的最优I/O复杂度是 O((N/B) log_{M/B} (N/B))。缓存无关模型中的经典算法是漏斗排序，这里介绍其简化版——懒惰漏斗排序。

该算法需要一个关键构件：K-漏斗。

定义：一个K-漏斗是一个数据结构，它能合并K个已排序的输入列表（列表总大小至多为K³），并输出完整排序结果。
性能声明：K-漏斗可以使用 O(K²) 空间，并以 O(K³/B * log_{M/B}(K³/B) + K) 的I/O复杂度完成合并。

排序算法流程：

将包含N个元素的数组划分为 K = N^{1/3} 组，每组大小为 N^{2/3}。
递归排序每一组。
使用一个 K-漏斗 来合并这K个已排序的组。

递归式与复杂度：通过设立递归式（基础情况是当 N ≤ B² 时，整个数组可读入内存排序），并利用高缓存假设（M ≥ B²，或更一般地 M = Ω(B^{1+γ})），可以解出总I/O复杂度为 O((N/B) log_{M/B} (N/B))，匹配了外部存储模型的最优界。

K-漏斗的实现简述：
K-漏斗本身也采用递归结构：

将一个K-漏斗在中间高度处分割，得到一个顶层的 √K-漏斗和 √K 个底层的 √K-漏斗。
连接顶层和底层漏斗的边上有缓冲区。
合并操作是“懒惰”进行的：当一个节点的输出缓冲区需要填充时，它从其子节点的缓冲区中取最小元素；如果子缓冲区为空，则递归地触发子节点填充其自身的缓冲区。
通过递归布局和仔细分析（特别是聚焦于大小约等于 √M 的递归层），可以证明其声明的性能。

重要：高缓存假设对于在缓存无关模型中达到此排序复杂度是必要的。

总结

本节课中我们一起学习了缓存无关模型。我们了解了它与外部存储模型的核心区别：算法对内存参数M和块大小B不可知，且不控制缓存替换。我们通过在线算法中的分页问题理论，论证了“操作系统做出最优替换”这一分析假设的合理性。接着，我们探讨了在该模型下的一系列算法和数据结构，包括数组扫描、矩阵乘法、链表、静态二叉搜索树（vEB布局）以及排序（懒惰漏斗排序）。这些设计通过巧妙的递归布局和分治策略，使得算法无需知晓底层硬件参数，就能自动适应内存层次结构，并获得接近最优的I/O性能。

25：大规模并行计算 🚀

在本节课中，我们将要学习大规模并行计算模型，特别是MapReduce模型。我们将探讨如何在此模型下设计高效算法，并分析其资源使用情况。课程将涵盖排序、最小生成树和三角形计数等经典问题的MapReduce解决方案。

模型介绍：MapReduce

上一节我们介绍了大规模并行计算的背景，本节中我们来看看一个具体且广泛应用的模型——MapReduce。

MapReduce是一个用于处理海量数据集的编程模型。计算被分解为多轮进行。

以下是MapReduce一轮计算中的三个阶段：

Map阶段：每个输入项由一个映射函数处理。映射函数接收一个键值对 <k, v> 作为输入，并输出一组新的键值对。
```
# 伪代码示例
def map_function(key, value):
    # 处理逻辑
    emit(new_key, new_value)
```

Shuffle阶段：此阶段对程序员透明，由系统内部处理。系统将所有在Map阶段输出的项，根据其键进行分组。具有相同键的项会被发送到同一个Reducer。

Reduce阶段：每个Reducer接收一个键 K 和与该键关联的一系列值 [V1, V2, ..., Vr]。Reducer对这些值进行计算，并输出一组新的项。

在这个模型中，我们的设计目标包括：

使用较少的计算轮数。
每个Reducer使用的内存远小于总数据量 n。
所有机器的总内存使用量远小于 n^2，理想情况下接近线性。
总工作量或最大单机工作量较小。

算法示例一：排序 (TeraSort)

现在，我们来看如何在MapReduce模型中解决排序问题。我们将介绍一个名为TeraSort的两轮算法。

该算法的核心思想是，如果我们有 P 台机器，我们希望将最小的 n/P 个元素发送到机器1，第二小的 n/P 个元素发送到机器2，依此类推。每台机器在本地对自己的数据进行排序，最后将所有结果串联起来。

问题在于我们不知道这些分界点具体是哪些元素。解决方案是通过采样来近似地找到这些分界点。

以下是TeraSort算法的两轮伪代码描述：

第一轮

Map函数：对于每个输入元素 (i, a_i)。
1. 将其发送到由 i mod P 决定的Reducer，以便初步分区。
2. 以概率 T/n 对该元素进行采样（T 是一个参数，例如 T = log(P)/ε^2）。如果被采样，则将该元素的值发送给所有 P 个Reducer。
Reduce函数：每个Reducer j 收到两类数据：初步分到本机的元素 B，以及所有采样元素 S。
1. 对采样集合 S 进行排序。
2. 以排序后的 S 作为参考，为 B 中的每个元素 a_i 确定其最终应归属的Reducer（例如，通过二分查找确定 a_i 在 S 中的排名区间）。
3. 输出新的键值对 (r, (i, a_i))，其中 r 是计算出的目标Reducer编号。

第二轮

Map函数：恒等函数，直接传递 (r, (i, a_i))。
Reduce函数：Reducer j 收到所有标记为 j 的元素，在本地对这些元素进行排序，并输出结果。

算法分析
通过设置参数 ε = 1/P 和 T = log(P)/ε^2，并利用切尔诺夫界限进行概率分析，可以证明：

每台机器在第二轮需要排序的元素数量最多为 O(n/P)。
当机器数量 P 约为 n^(1/3) 时，每台机器所需内存约为 n^(2/3)，实现了亚线性的内存使用和并行时间。

实际的TeraSort实现针对整数键进行了优化（例如使用Trie树加速查找），但上述简化版本揭示了算法的核心思想。

算法示例二：最小生成树 (MST)

接下来，我们探讨一个图论问题：在边数 m = n^(1+c)（c>0，即非稀疏图）的图中寻找最小生成树。目标是使每台机器的内存使用远小于 m。

该算法由Karloff和Vassilvitskii提出，其核心思想是随机划分顶点，然后在所有顶点子集的诱导子图上并行计算最小生成森林。

算法描述

选择一个参数 k（例如 k = n^(c/2)）。
使用哈希函数 h 将顶点随机划分到 k 个组中：V1, V2, ..., Vk。
对于每一对 (i, j)，其中 1 ≤ i, j ≤ k，考虑由 Vi ∪ Vj 诱导的子图 Gij。
为每个 Gij 分配一个Reducer，计算其最小生成森林 Mij。
最后，一个中心Reducer收集所有 Mij 中的边，并在这个边的集合上计算整个图的最小生成树。

MapReduce实现

第一轮 Map：对于每条边 (u, v)，将其发送到键为 (h(u), h(v)) 的Reducer。如果 h(u) == h(v) = i，则需要将该边发送给所有形如 (i, *) 和 (*, i) 的Reducer，以确保子图 Gij 包含 Vi 和 Vj 内部的所有边。
第一轮 Reduce：Reducer (i, j) 收到子图 Gij 的所有边，计算其最小生成森林 Mij，并将 Mij 中的所有边输出到一个公共键（如 $）下。
第二轮：Map阶段传递数据，一个Reducer收集所有 Mij 的边，并计算最终的全图MST。

资源分析

每机内存：每个 Gij 的期望边数为 O(m/k^2)。通过更精细的度分组分析，可以高概率保证每机内存为 O~(n^(1+c/2))，远小于 m。
总内存：所有 Mij 的边数总和约为 O(k^2 * (n/k)) = O(kn)，当 k = n^(c/2) 时，总内存约为 O(n^(1+c/2))，虽然小于 m^2，但并非显著优于简单方法。后续有更优的算法进一步降低了总内存消耗。

算法示例三：三角形计数

最后，我们看一个更复杂的图问题：计算图中三角形的数量。一个简单的顺序算法可以在 O(m^(3/2)) 时间内完成。

顺序算法思路

为顶点定义一个全序（例如，按度排序，度相同则按ID排序）。
遍历每个顶点 v，只考虑其邻接点中序大于 v 的点 u 和 w，然后检查边 (u, w) 是否存在。通过设置度阈值进行平衡，可得 O(m^(3/2)) 的运行时间。

MapReduce实现
该算法可以自然地映射到MapReduce框架：

第一轮：Mapper发送每条边 (u, v)（确保 u 在序中大于 v）到以 v 为键的Reducer。Reducer收到顶点 v 及其“高阶”邻居列表 S，然后检查 S 中是否存在边（可通过将边列表也发送到相应键来实现或预加载），并为每个找到的三角形候选 (v, u, w) 发出信号。
第二轮：Mapper对不同类型的中间结果（边或三角形候选）进行标记和路由。Reducer收集信息，确认三角形并去重计数。
最终可能需要一轮聚合来求和。

通过精心设计键和数据的路由，可以在保证每机内存较小的前提下，在常数轮内完成三角形计数。

总结

本节课中我们一起学习了大规模并行计算中的MapReduce模型。

我们首先了解了MapReduce的三阶段计算模型和设计目标。
接着，我们深入分析了三个经典问题的MapReduce算法：用于排序的TeraSort算法，它通过采样近似分界点来实现负载均衡；用于稠密图的最小生成树算法，它通过随机顶点划分实现并行子图计算；以及三角形计数算法，它将顺序算法巧妙转化为多轮MapReduce作业。
这些案例展示了如何在限制每机内存、控制通信轮数的前提下，设计高效的大规模并行算法。MapReduce及其开源实现（如Hadoop）已成为工业界处理大数据的重要工具，其背后的算法思想对分布式计算至关重要。

posted @ 2026-03-26 12:31 布客飞龙IV 阅读(2) 评论(0) 收藏举报

刷新页面返回顶部

龙哥盟

哈佛-CS229r-大数据算法笔记-全-

哈佛 CS229r 大数据算法笔记（全）

1：课程介绍与近似计数算法 🚀

课程概述与物流信息 📋

什么是大数据算法？ 🤔

课程主题概览 📚

第一个流式算法：近似计数 🧮

莫里斯算法

改进算法：Morris+ 和 Morris++

总结 ✨

2：流式算法与不同元素问题 🧮

回顾与引言

不同元素问题（F0问题）🔍

问题定义

两种简单算法

近似目标

理想化算法：Flajolet-Martin 算法

算法步骤

算法原理与期望

方差分析

从理想到现实：改进算法

FM+ 算法：通过平均降低方差

FM++ 算法：通过中位数提升成功概率

实现真正的有限空间：哈希函数与层次化算法

有限独立哈希函数族

一个常数因子近似算法

提升到 (1+ε) 近似：层次化算法

总结与下节预告

3：下界与转盘模型 🧮

下界：为什么我们需要随机化和近似？ 🚫

确定性精确算法需要线性空间

确定性近似算法也需要线性空间

随机化精确算法也需要线性空间

转盘模型与线性草图 🌀

什么是转盘模型？

核心算法技术：线性草图

矩估计问题 📊

p=2 时的算法：AMS草图

扩展到 0 < p ≤ 2：Indyk算法 🧪

总结 🎯

4：Lp范数估计（下）与Nisan伪随机生成器 🧮

回顾：Indyk算法与p稳定分布

分析Indyk算法的正确性

随机性挑战与Nisan伪随机生成器 🎲

Nisan伪随机生成器简介

将Nisan生成器应用于Indyk算法

对于p > 2的情况

总结

5：P>2范数估计与点查询问题 🧮

P>2范数估计回顾

算法直观理解与游戏计划 🎯

严格分析

定义与引理1：重坐标数量的期望

控制轻坐标的噪声

结论

从范数估计到点查询与频繁项 🔍

确定性点查询算法简介

总结 📝

6：点查询与频繁项检测 🎯

点查询问题 📍

频繁项检测问题 🔍

使用 ε-不相干矩阵求解 🧮

不相干矩阵的构造方法

从点查询到频繁项检测的黑盒转换

随机化方法：Count-Min Sketch 🎲

基于树的快速频繁项检测算法 🌳

更强的保证：L1/L1 稀疏恢复 🏆

总结 📝

7：L0采样与图连通性算法 🧮➡️🌳

L0采样：基本概念与构建 🎯

什么是Lp采样？

什么是L0采样？

如何构建L0采样器？🔨

组件A：恢复严格1-稀疏向量

组件B：检测非1-稀疏向量

组件C：几何采样与最终整合

应用：动态图连通性 🌉

挑战与目标

图草图算法

`p=2` 时的算法：AMS草图

扩展到 `0 < p ≤ 2`：Indyk算法 🧪