CS267 · 第 17 课

随机化算法

随机化算法在内部使用随机性来加速计算、提升可扩展性或增强可靠性。本课程分为两部分：第一部分介绍随机化数值线性代数（Randomized NLA），通过线性草图（sketching）将高维问题压缩到低维子空间；第二部分介绍随机化在数据科学与 HPC 中的应用，包括加权采样、数据流中的基数估计，以及 Luby 的并行最大独立集算法。

---

概述

随机化算法通过引入随机性，在保持可控误差的前提下显著降低计算成本。与传统确定性算法不同，随机化算法的正确性通常以概率保证（例如 (ε, δ) 近似）。本课程涵盖两类核心应用：

• 随机化数值线性代数（Randomized NLA）：通过随机投影将矩阵/向量压缩到低维空间，保留关键几何信息，从而加速最小二乘、低秩逼近等问题。
• 随机化数据科学算法：在数据流、采样、并行图算法等场景中利用随机性实现亚线性复杂度或高效并行化。

---

第一部分：随机化数值线性代数

线性草图（Linear Sketching）

线性草图的核心思想是将原始问题通过一个随机线性变换 S 压缩到低维空间。给定一个高维向量 x ∈ ℝⁿ，其草图为 Sx ∈ ℝᵐ，其中 m ≪ n。

典型形式：对于矩阵 A ∈ ℝⁿˣᵈ，构造草图 SA，其中 S ∈ ℝᵐˣⁿ 是一个随机矩阵。草图保留了 A 的列空间的关键几何性质。

子空间嵌入（Subspace Embedding）

一个矩阵 S 被称为 A 的 ε-子空间嵌入（ε-subspace embedding），如果对所有 x ∈ ℝᵈ 满足：

(1 - ε) ||Ax||² ≤ ||SAx||² ≤ (1 + ε) ||Ax||²

这保证了草图 SA 保留了 A 的列空间中所有向量的范数（即内积和距离）至多 ε 倍的相对误差。

Oblivious Subspace Embedding 定理

定理（Oblivious Subspace Embedding, OSE）：设 S ∈ ℝᵐˣⁿ 的每个元素独立地取值为 ±1/√m（等概率）。如果 m = O(d²/ε²)，那么 S 以至少 2/3 的概率成为任意固定 d 维子空间的 ε-子空间嵌入。

如何选择 S：快速结构化矩阵

实际应用中，我们需要在以下性质之间权衡：

矩阵类型	m 的阶	应用代价 SA	并行性
稠密随机 ±1	O(d²/ε²)	O(nmd)	良好
SHRT（Subsampled Randomized Hadamard Transform）	O(d²/ε²)	O(nd log n)	良好
CountSketch	O(d²/ε²)	O(nnz(A))	优秀

Fast Johnson-Lindenstrauss Transform / SHRT

SHRT（Subsampled Randomized Hadamard Transform） 将 S 分解为三个算子的乘积：

S = PHD

• D：对角矩阵，对角线元素为随机 ±1（随机符号翻转，消除方向偏好）。
• H：归一化的 Hadamard 矩阵（利用快速 Hadamard 变换，FHT）。
• P：随机采样矩阵，从 n 行中均匀采样 m 行。

应用代价仅为 O(n log n)（由 FHT 决定），远优于稠密矩阵的 O(nm)。

import numpy as np
from scipy.linalg import hadamard

def shrt_sketch(A, m):
    """应用 SHRT 到矩阵 A，返回 m 行的草图"""
    n, d = A.shape
    # Step 1: 随机符号翻转
    D = np.random.choice([-1, 1], size=n)
    AD = A * D[:, np.newaxis]
    # Step 2: 快速 Hadamard 变换
    H = hadamard(n) / np.sqrt(n)
    AHD = H @ AD
    # Step 3: 随机采样 m 行
    idx = np.random.choice(n, size=m, replace=False)
    return AHD[idx, :] / np.sqrt(m / n)

Block SHRT（BSHRT）：并行化

为了在并行环境中高效应用 SHRT，可以使用 Block SHRT（BSHRT）。核心思想是将矩阵按行分块，每块独立应用 SHRT，然后合并结果。BSHRT 采用 1D 外积算法：

• 每个进程持有矩阵的一个行块 A_i。
• 每个进程独立应用 SHRT 到本地块。
• 通过 all-gather 或 reduce-scatter 合并草图。

最小二乘问题：Sketch-and-Solve

对于超定最小二乘问题 min_x ||Ax - b||₂，经典解法需要 O(nd²) 时间。Sketch-and-Solve 方法：

1. 构造草图 SA 和 Sb（S 是子空间嵌入）。
2. 求解压缩后的问题：x* = min_x ||SAx - Sb||₂。
3. 返回 x* 作为原问题的近似解。

(ε, δ) 保证：以至少 1-δ 的概率，||Ax* - b||₂ ≤ (1+ε) · min_x ||Ax - b||₂。

import numpy as np

def sketch_and_solve(A, b, m):
    """Sketch-and-Solve 求解最小二乘问题"""
    n, d = A.shape
    # 构造高斯随机草图矩阵
    S = np.random.randn(m, n) / np.sqrt(m)
    SA = S @ A
    Sb = S @ b
    # 求解压缩后的问题
    x_star, _, _, _ = np.linalg.lstsq(SA, Sb, rcond=None)
    return x_star

---

第二部分：数据科学与 HPC 中的随机化算法

从加权分布中采样

许多应用需要从离散分布 p = (p₁, p₂, ..., pₙ) 中按权重采样。常用方法：

方法	预处理	单次采样	备注
逆变换采样（Inverse Transform）	O(n)	O(log n)（二分查找）	简单通用
拒绝采样（Rejection Sampling）	O(1)	O(1) 期望	需要提议分布
别名方法（Alias Method）	O(n)	O(1)	最优单次采样

别名方法（Alias Method）

别名方法通过将概率分布重构为 n 个"桶"，每个桶最多包含两个元素，实现 O(1) 采样：

import numpy as np

class AliasMethod:
    """别名方法：O(n) 预处理，O(1) 采样"""
    def __init__(self, probs):
        n = len(probs)
        self.prob = np.zeros(n)
        self.alias = np.zeros(n, dtype=int)
        scaled = np.array(probs) * n
        small, large = [], []
        for i, p in enumerate(scaled):
            (small if p < 1.0 else large).append(i)
        while small and large:
            s, l = small.pop(), large.pop()
            self.prob[s] = scaled[s]
            self.alias[s] = l
            scaled[l] = scaled[l] + scaled[s] - 1.0
            (small if scaled[l] < 1.0 else large).append(l)
        for idx in large + small:
            self.prob[idx] = 1.0

    def sample(self):
        i = np.random.randint(len(self.prob))
        return i if np.random.rand() < self.prob[i] else self.alias[i]

数据流中不同元素的计数

给定一个数据流，估计其中不同元素的数量（cardinality estimation），要求空间远小于 n。

PCSA / Flajolet-Martin 算法

Probabilistic Counting with Stochastic Averaging (PCSA) 的核心思想：

• 对每个元素应用哈希函数 h，观察 h(x) 的二进制表示中末尾 0 的个数（即 ρ(h(x)）。
• 维护一个位图（bitmap），记录遇到的最大 ρ 值。
• 估计值 ≈ 2^R / φ，其中 R 是最大 ρ 值，φ ≈ 0.77351。

HyperLogLog：调和平均改进

HyperLogLog 使用调和平均（harmonic mean） 替代 PCSA 的几何平均，显著降低异常值的影响：

• 将流分成 m 个子流（使用哈希值的前 log₂m 位作为桶索引）。
• 每个桶维护最大 ρ 值。
• 使用调和平均合并桶的估计值。

误差率：O(1/√m)，空间复杂度：O(m log log n)。

import hashlib
import numpy as np

def rho(x):
    """计算二进制表示中末尾 0 的个数 + 1"""
    if x == 0:
        return 64
    n = 0
    while (x & 1) == 0:
        n += 1
        x >>= 1
    return n + 1

class HyperLogLog:
    """HyperLogLog 基数估计器"""
    def __init__(self, b=10):
        self.b = b          # 桶的指数
        self.m = 1 << b     # 桶的数量
        self.registers = np.zeros(self.m, dtype=np.int32)
        self.alpha = self._correction(self.m)

    def _correction(self, m):
        if m == 16: return 0.673
        if m == 32: return 0.697
        if m == 64: return 0.709
        return 0.7213 / (1.0 + 1.079 / m)

    def add(self, item):
        h = int(hashlib.md5(str(item).encode()).hexdigest()[:16], 16)
        idx = h >> (64 - self.b)
        self.registers[idx] = max(self.registers[idx], rho(h))

    def estimate(self):
        Z = np.sum(2.0 ** (-self.registers))
        E = self.alpha * self.m ** 2 / Z
        return int(E)

PCSA 中的并行化

PCSA 和 HyperLogLog 天然支持并行化：

• 可合并草图（Mergeable Sketches）：每个进程维护本地草图。
• 合并操作：位图的合并等价于按位 OR（bitwise OR）。
• 全局聚合：使用 all-reduce（按位 OR）合并所有进程的草图。

from mpi4py import MPI

def parallel_hll(local_stream, b=10):
    """并行 HyperLogLog 基数估计"""
    comm = MPI.COMM_WORLD
    rank = comm.Get_rank()

    # 每个进程本地估计
    hll = HyperLogLog(b)
    for item in local_stream:
        hll.add(item)
    local_regs = hll.registers.copy()

    # all-reduce: 按位取最大值（等价于位图的 OR）
    global_regs = np.zeros_like(local_regs)
    comm.Allreduce([local_regs, MPI.INT],
                   [global_regs, MPI.INT],
                   op=MPI.MAX)

    # 在 rank 0 上计算最终估计
    if rank == 0:
        hll.registers = global_regs
        return hll.estimate()
    return None

最大独立集（Maximal Independent Set, MIS）

问题：给定无向图 G = (V, E)，找一个极大独立集（即一组互不相邻的顶点，且无法再添加任何顶点）。

顺序算法的困境

贪心 MIS 算法顺序复杂度为 O(n)，但并行跨度（span）为 O(n)——本质上是顺序的。

Luby 的并行随机化算法

Luby 算法（1986） 通过随机化实现 O(log n) 跨度：

1. 每个顶点独立地随机生成一个标签 r(v) ∈ (0, 1)。
2. 如果一个顶点 v 的标签比其所有邻居都小（局部最小值），则 v 加入独立集。
3. 移除 v 及其所有邻居（它们不能再加入）。
4. 在剩余图上重复上述过程。

import random
from collections import defaultdict

def luby_mis(n, adj):
    """Luby 的并行 MIS 算法"""
    in_graph = set(range(n))
    mis = set()

    while in_graph:
        # Step 1: 为每个顶点分配随机标签
        labels = {v: random.random() for v in in_graph}

        # Step 2: 局部最小值加入 MIS
        chosen = set()
        for v in in_graph:
            if all(labels[v] < labels[u] for u in adj[v] if u in in_graph):
                chosen.add(v)

        mis.update(chosen)

        # Step 3: 移除 chosen 及其邻居
        to_remove = set(chosen)
        for v in chosen:
            to_remove.update(u for u in adj[v] if u in in_graph)
        in_graph -= to_remove

    return mis

复杂度分析：

• 每轮期望移除至少一半的边。
• 期望轮数：O(log n)。
• 每轮工作量：O(m)（检查邻居）。
• 总期望工作量：O(m log n)，跨度：O(log n)。

---

Bonus：Balls into Bins

经典的"球入箱"问题展示了随机化的力量：

• 纯随机：将 n 个球随机放入 n 个箱，最重箱的负载期望为 O(ln n / ln ln n)。
• Power of 2 Choices：每个球随机选 2 个箱，放入负载较轻者，最重箱负载降为 O(ln ln n)。

import random

def balls_into_bins(n, two_choices=False):
    """模拟 balls into bins 实验"""
    bins = [0] * n
    for _ in range(n):
        if two_choices:
            i, j = random.randrange(n), random.randrange(n)
            idx = i if bins[i] <= bins[j] else j
        else:
            idx = random.randrange(n)
        bins[idx] += 1
    return max(bins)

# 示例：n=10000 时
n = 10000
print(f"纯随机: {balls_into_bins(n, False)}")
print(f"2-choices: {balls_into_bins(n, True)}")

策略	最重箱负载（n 球入 n 箱）	改进
纯随机	Θ(ln n / ln ln n)	基线
Power of 2 Choices	Θ(ln ln n)	指数改进
Power of d Choices	Θ(ln ln n / ln d)	进一步改进

---

本课要点

---

知识检查