# 实验 07：递归模型索引（RMI）

> 用一串很小的数学模型去拟合"键值 → 数组下标"这条曲线，把有序数组上的查找从"折半比较"改成"先预测、再在很小的邻域里修正"。

## 一、实验目标

完成本实验后，你应该能够：

1. 说出学习索引（learned index）与 B 树在"看待索引的方式"上的根本差别：前者把索引当成一个函数拟合问题，后者把索引当成一个"分层导航 + 最后一页精确查找"的问题。
2. 独立写出四种查找算法：二分查找、指数搜索（exponential search / galloping search）、单层 RMI、多层 RMI，以及带误差界（error bound）的 RMI。
3. 用**平均绝对误差**与**最大绝对误差**两个指标量化"模型预测得准不准"，并解释为什么这两个指标关注点不同。
4. 用 ns/次 为单位测量各算法的单查询耗时，并正确解释"为什么我手写的 RMI 跑不过 `bisect`"。
5. 说明 RMI 的适用边界：静态、内存内、有序、一维。

## 二、预备知识

### 2.1 下界查找（lower bound）

有序数组里可能有重复值，所以"查找"通常不是找"有没有"，而是找**第一个不小于目标值 x 的元素下标**。这个位置叫下界（lower bound）。

```
数组 a = [1, 2, 2, 2, 5, 7]

lower_bound(a, 2) = 1      # 第一个 >= 2 的元素在下标 1
lower_bound(a, 3) = 4      # 第一个 >= 3 的元素是 5，在下标 4
lower_bound(a, 8) = 6      # 没有 >= 8 的元素，返回 len(a)
```

记法：本文档用 `N` 表示数组长度，下标从 0 到 `N-1`，`len(a)` 表示"未找到"。

### 2.2 二分查找

每次取中点比较，把搜索区间减半。

```
初始区间  [lo, hi)         宽度 = N
第1次比较 [lo, mid) 或 [mid+1, hi)   宽度 ≈ N/2
第2次比较 ...                        宽度 ≈ N/4
...
第k次比较                            宽度 ≈ N / 2^k
```

比较次数：

```
k = ⌈log₂(N + 1)⌉
```

`N = 200,000` 时，`log₂(200001) ≈ 17.61`，即约 18 次比较。**注意**：这个公式只统计"比较次数"，不统计内存访问。

### 2.3 指数搜索

二分查找的搜索区间是"整段"。如果我们事先知道目标离某个位置 `p` 很近，就先用步长 1、2、4、8… 向外扩，把目标"框"进一个小区间，再在小区间里二分。

```
  预测位置 pred
       |
  ...|  |  |   |     |         |                 |...
     +1 +2  +4   +8     +16        +32               +64
     扩 扩  扩   扩      扩         扩                 扩   ← 步长翻倍
     └──────────── 找到区间后在此二分 ─────────────┘
```

若真实位置与预测位置的距离为 `d`，则扩区间需要约 `⌈log₂ d⌉` 步，区间内二分再需约 `⌈log₂ d⌉` 步，合计约 `2⌈log₂ d⌉ + 1` 次比较。**`d` 越小，指数搜索越划算**，这正是 RMI 要努力减小 `d` 的原因。

### 2.4 线性样条模型（linear spline）

一条直线，穿过点集的首尾两点。给一组输入 `X = [x₀, x₁, ..., x_{m-1}]` 和目标 `Y = [y₀, ..., y_{m-1}]`（本实验中 `Y` 就是下标）：

```
斜率      a = (Y[m-1] - Y[0]) / (X[m-1] - X[0])
截距      b = Y[0] - a · X[0]
预测      ŷ = f(x) = a · x + b
```

符号含义：`X` 为按键升序排列的键值；`Y` 为对应的数组下标；`a` 的单位是"下标/键值"，`b` 的单位是"下标"。当 `X[m-1] == X[0]`（该段所有键相同）时除零，实现里必须特判。

只有两个点参与拟合，所以**训练是 O(1) 的**——这是 RMI 能做得极快的前提。

### 2.5 关键类比：索引就是累积分布函数（CDF）

这是整个实验的思想核心。给定一个升序数组 `D`（长度 `N`），定义它的经验累积分布函数：

```
F(x) = |{ k ∈ D : k ≤ x }| / N        （取值在 [0, 1] 之间）
```

`F(x)` 的含义是"数组里不大于 x 的元素所占的比例"。那么 `x` 在数组中的位置就近似是：

```
pos(x) ≈ F(x) · N
```

举例：`N = 200,000`，若 `F(x) = 0.4`，则 `pos(x) ≈ 80,000`。

于是：

| 传统索引观 | 学习索引观 |
| --- | --- |
| 索引是一棵用于导航的树 | 索引是一个拟合 `F(x)` 的函数 |
| 查找靠逐层比较 | 查找靠一次函数求值 |
| 命中叶节点后在小页内二分 | 预测后用误差邻域内二分 |
| 复杂度与数据分布无关 | 预测误差依赖数据分布 |

这一观点出自 Kraska et al., *The Case for Learned Index Structures*, SIGMOD 2018（arXiv:1712.01208）。

### 2.6 RMI 的分层结构

RMI（Recursive Model Index，递归模型索引）由若干"层"组成，每层有若干模型（本实验全部用线性样条）。

```
                        第 0 层（根）
                     ┌──────────────┐
       key ────────► │   model 0    │  预测出的位置 → 决定走哪个下层模型
                     └──────┬───────┘
                            │
          ┌─────────────────┼─────────────────┐
          ▼                 ▼                 ▼
   ┌────────────┐    ┌────────────┐    ┌────────────┐   第 1 层（末层）
   │  model 0   │    │  model 1   │    │  model 2   │   ... 共 M 个
   └─────┬──────┘    └─────┬──────┘    └─────┬──────┘
         │                 │                 │
         ▼                 ▼                 ▼
      预测位置           预测位置          预测位置
         └─────────────────┴─────────────────┘
                          │
                  局部搜索（指数/二分）
                          │
                          ▼
                      真实下标
```

要点：

- 第 0 层恒为 1 个模型（根模型）。
- 用 `layer_sizes` 描述每层模型数，例如 `[1, 1000]` 表示两层、末层 1000 个模型。
- 层数 ≥ 3 时，中间层只负责"路由"，只有末层模型负责给出最终位置预测。
- 每个模型只覆盖一小段键域，段越窄，直线拟合越准 → 误差越小。这就是"加模型能降误差"的直觉来源。

## 三、实验原理

### 3.1 训练：自顶向下分配键

训练是一个逐层进行的过程：

1. 根模型在**全部** `(key, index)` 对上训练。
2. 用刚训练好的模型预测它所见到的每个键的位置 `pred`，再按比例映射到一个下层模型编号：
   ```
   next_model = ⌊ pred / N × layer_sizes[layer+1] ⌋
   ```
3. 把该键分配给 `next_model`。
4. 一层训练完，进入下一层，直到末层。

`pred / N` 是"预测位置占全表的比例"，乘以"下层模型总数"就得到"应该落到第几个下层模型"。理想情况下每个末层模型分到 `N / M` 个键（等深划分）。

这条规则带来一个很好的性质：**查询数据集中真实存在的键时，只会用到确实训练过它的那些模型**。

### 3.2 查询：预测 + 局部修正

```
value
  │
  ├─ 第 0 层：pred₀ = f₀(value)
  │            model = ⌊pred₀ / N × 下层模型数⌋
  ├─ 第 1 层：pred₁ = f₁(value)           （若为末层则停止）
  │            ...
  └─ 末层：    pred = clamp(f_last(value), 0, N-1)
               │
               ├─ 无误差界：在 pred 处做指数搜索
               └─ 有误差界：在 [pred - err, pred + err] 内做二分
```

`clamp(v, 0, N-1)` 表示 `max(0, min(N-1, v))`，防止预测越界。

### 3.3 误差界（error bound）

训练完成后，回扫一遍数据集，为每个末层模型 `j` 记录它训练样本上的**最大绝对误差**：

```
err_j = max | clamp(f_j(key)) − true_index |      （对该模型分到的所有 key 取最大）
```

查询时用 `err_j` 把搜索范围锁死：

```
lo = max(0, pred − err_j)
hi = min(N, pred + err_j + 1)
```

区间宽度 `2·err_j + 1`，二分需要约 `⌈log₂(2·err_j + 1)⌉` 次比较。

对比：末层 1000 个模型时本数据集实测最大绝对误差为 1555（见 6.2 表），区间宽度 3111，`⌈log₂ 3111⌉ = 12` 次比较；而整表二分 `N = 200,000` 需要约 18 次，省 6 次。数据分布越平滑、误差越小，这个差距越大——这就是误差界的价值。

**注意一处易错点**：`err_j` 是在**训练样本**上统计的。如果查询的键不在数据集中（未命中查询），真实误差可能超出 `err_j`，此时区间内二分会返回错误结果。工程实现需要额外处理（例如区间端点处做边界检查，或退化为整表二分）。

### 3.4 为什么本实验里 RMI 跑不过 bisect

必须讲清楚，否则学生会得出错误结论：

- 本实验的 RMI 是**纯 Python**实现。一次预测要经历：方法调用 → 属性查找 → 浮点乘加 → `int()` → `max/min` → 又一次函数调用（指数搜索）→ 循环里的下标访问。Python 解释器的每条字节码都有固定开销。
- 实验里的 `bisect` 是 **C 实现**的 `bisect.bisect_left`，一次调用几乎不进入 Python 字节码循环。
- 所以"RMI 比二分快"这件事，在纯 Python 环境下测不出来。

RMI 的真实价值体现在另外两处，本实验可以观测到其中一处：

1. **预测精度**：误差随末层模型数单调下降（6.2 表），这正是学习索引的核心机制，与实现语言无关。
2. **缓存友好性**：RMI 查找时访问的内存地址是"预测位置附近的一小段"，而不是二分那样跳跃访问 `N/2`、`N/4`、`N/8`… 这些跨度极大的地址。二分的前若干次比较几乎每次都触发 cache miss。这一点在本实验的纯 Python 环境下被解释器开销掩盖了。

论文中报告的"相比 B 树 3× 以上加速"来自优化过的 C++ 实现，本实验不承诺复现该数字。

### 3.5 参考图件

| 图件 | 内容 |
| --- | --- |
| `../figures/fig_rmi_cdf.png` | 键的累积分布曲线 `F(x)` 与 `pos ≈ F(x)·N` 的对应关系 |
| `../figures/fig_rmi_error.png` | 不同末层模型数下预测位置与真实位置的偏差 |
| `../figures/fig_rmi_bench.png` | 各查找算法的 ns/次 对比 |

## 四、实验环境与准备

### 4.1 环境要求

| 项目 | 要求 |
| --- | --- |
| Python | 3.10 及以上 |
| 必装库 | 无（只用标准库 `random` / `bisect` / `statistics` / `time`） |
| 选装库 | `matplotlib`（步骤 8 画图）、`numpy`（步骤 8 画图，可选） |
| 参考实现 | `notebooks/Recursive-Model-Index.ipynb` |

安装命令（若缺）：

```bash
pip install matplotlib numpy
```

### 4.2 目录约定

把本指导书的全部代码放在同一个 `.py` 文件里（例如 `rmi_lab.py`）顺序执行，或在 Jupyter 中按步骤分格执行。步骤 2–7 是"搭积木"，步骤 8–9 是"测量"。

### 4.3 关于数据集规模

本指导书统一使用 `N = 200,000`。第 6 节给出的实测数据均在此规模下取得。若机器较慢，把 `N` 改成 `20,000` 也能跑通全部流程，但**数字会与第 6 节不同**，此时请填写你自己的实测值。

## 五、实验步骤

### 步骤 1：生成合成数据

数据由多项式加噪声生成，密度随键值增大而降低（这是 RMI 能"学到东西"的前提——如果数据完全均匀，一条直线就够了）。

```python
import random

N = 200_000

def poly_fn(x: int) -> int:
    """多项式函数：值随 x 二次增长，因此键域上稀疏、键域下稠密。"""
    return int(1 / 1000 * x ** 2)

random.seed(0)
data = []
for i in range(N):
    noise = random.randint(int(-N * 0.001), int(N * 0.001))
    data.append(poly_fn(i) + noise)

data = sorted(data)          # RMI 只作用于已排序数组
print("N =", len(data))
print("min =", data[0], " max =", data[-1])
```

画直方图看分布（需要 `matplotlib`）：

```python
import matplotlib.pyplot as plt

fig, ax = plt.subplots()
ax.hist(data, bins=20)
ax.set_xlabel("Key")
ax.set_ylabel("Frequency")
plt.show()
```

你会看到左侧频次高、右侧频次低——即"密度随键值增大而降低"。

### 步骤 2：基线算法——二分查找

```python
def lower_bound_binary_search(a, value, start=None, lo=None, hi=None):
    """返回 a 中第一个 >= value 的下标；不存在则返回 len(a)。

    lo/hi 可限定搜索区间 [lo, hi)；start 可指定第一次比较的位置。
    """
    if lo is None:
        lo = 0
    if hi is None:
        hi = len(a)
    assert 0 <= lo <= hi <= len(a)

    if start is not None:
        assert lo <= start < hi
        if a[start] < value:
            lo = start + 1
        else:
            hi = start

    while lo < hi:
        mid = lo + (hi - lo) // 2
        if a[mid] < value:
            lo = mid + 1
        else:
            hi = mid
    return lo


# 自检：与线性查找对拍
def lower_bound_linear_search(a, value):
    for i, v in enumerate(a):
        if v >= value:
            return i
    return len(a)

small = sorted([3, 1, 4, 1, 5, 9, 2, 6, 5, 3, 5])
for x in range(0, 11):
    assert lower_bound_binary_search(small, x) == lower_bound_linear_search(small, x)
print("二分查找自检通过")
```

### 步骤 3：线性样条模型

```python
def is_sorted(lst):
    return all(lst[i] <= lst[i + 1] for i in range(len(lst) - 1))


class LinearSplineModel:
    """穿过 (X[0], Y[0]) 与 (X[-1], Y[-1]) 的直线模型。"""

    def __init__(self, X, Y):
        assert is_sorted(X), "X 必须升序"
        assert len(X) == len(Y)
        if len(X) == 0:
            self.slope_ = 0.0
            self.intercept_ = 0.0
        elif len(X) == 1:
            self.slope_ = 0.0
            self.intercept_ = float(Y[0])
        elif X[-1] == X[0]:          # 防止除零
            self.slope_ = 0.0
            self.intercept_ = float(Y[-1])
        else:
            self.slope_ = (Y[-1] - Y[0]) / (X[-1] - X[0])
            self.intercept_ = Y[0] - self.slope_ * X[0]

    def predict(self, x):
        return self.slope_ * x + self.intercept_


m = LinearSplineModel([0, 10, 20], [0, 1, 2])
print("slope =", m.slope_, "intercept =", m.intercept_)
print("predict(5) =", m.predict(5))          # 期望 0.5
```

### 步骤 4：指数搜索

```python
def lower_bound_exponential_search(a, value, start=None, lo=None, hi=None):
    """从 start 出发按 1,2,4,8... 扩区间，再在区间内二分。"""
    if lo is None:
        lo = 0
    if hi is None:
        hi = len(a)
    if start is None:
        start = lo
    assert 0 <= lo <= start < hi <= len(a)

    bound = 1
    prev = start

    if a[start] < value:                       # 向右扩
        curr = min(prev + bound, hi)
        while curr < hi and a[curr] < value:
            bound *= 2
            prev = curr
            curr = min(prev + bound, hi)
        return lower_bound_binary_search(a, value, lo=prev, hi=min(curr + 1, len(a)))
    else:                                      # 向左扩
        curr = max(prev - bound, lo)
        while curr > lo and a[curr] >= value:
            bound *= 2
            prev = curr
            curr = max(prev - bound, lo)
        return lower_bound_binary_search(a, value, lo=curr, hi=prev)


for x in range(0, 11):
    for s in range(len(small)):
        assert lower_bound_exponential_search(small, x, start=s) == \
               lower_bound_linear_search(small, x)
print("指数搜索自检通过")
```

### 步骤 5：单层 RMI

```python
class SingleLayerRMI:
    """一个线性样条模型 + 指数搜索修正。"""

    def __init__(self, sorted_list):
        assert len(sorted_list) > 0 and is_sorted(sorted_list)
        self.N = len(sorted_list)
        self.model = LinearSplineModel(sorted_list, list(range(self.N)))

    def _predict(self, value):
        return max(0, min(self.N - 1, int(self.model.predict(value))))

    def lower_bound_search(self, sorted_list, value):
        pred = self._predict(value)
        return lower_bound_exponential_search(sorted_list, value, start=pred)


single = SingleLayerRMI(data)
v = 1337
print("二分:", lower_bound_binary_search(data, v))
print("单层RMI:", single.lower_bound_search(data, v))
assert single.lower_bound_search(data, v) == lower_bound_binary_search(data, v)
```

单次预测误差统计：

```python
import statistics

pred_single = [single._predict(k) for k in data]
err_single = [abs(p - i) for i, p in enumerate(pred_single)]
print("单层 RMI  平均绝对误差 =", round(statistics.mean(err_single), 1))
print("单层 RMI  最大绝对误差 =", max(err_single))
```

### 步骤 6：两层 RMI

```python
class MultiLayerRMI:
    """layer_sizes=[1, M] 表示两层：根模型 + M 个末层模型。"""

    def __init__(self, sorted_list, layer_sizes):
        assert len(sorted_list) > 0 and is_sorted(sorted_list)
        assert len(layer_sizes) > 0 and layer_sizes[0] == 1
        self.N = len(sorted_list)
        indexes = list(range(self.N))

        model_keys = [[[] for _ in range(s)] for s in layer_sizes]
        model_idx = [[[] for _ in range(s)] for s in layer_sizes]
        model_keys[0][0] = sorted_list
        model_idx[0][0] = indexes

        self.models = []
        for layer, layer_size in enumerate(layer_sizes):
            self.models.append([])
            for mi in range(layer_size):
                mdl = LinearSplineModel(model_keys[layer][mi], model_idx[layer][mi])
                self.models[layer].append(mdl)
                if layer < len(layer_sizes) - 1:
                    for key, idx in zip(model_keys[layer][mi], model_idx[layer][mi]):
                        pred = max(0, min(self.N - 1, int(mdl.predict(key))))
                        nxt = int(pred / self.N * layer_sizes[layer + 1])
                        nxt = min(nxt, layer_sizes[layer + 1] - 1)
                        model_keys[layer + 1][nxt].append(key)
                        model_idx[layer + 1][nxt].append(idx)

    def _get_last_layer_model(self, value):
        model, layer = 0, 0
        while layer < len(self.models) - 1:
            pred = max(0, min(self.N - 1, int(self.models[layer][model].predict(value))))
            model = int(pred / self.N * len(self.models[layer + 1]))
            model = min(model, len(self.models[layer + 1]) - 1)
            layer += 1
        return model

    def _predict(self, value):
        j = self._get_last_layer_model(value)
        return max(0, min(self.N - 1, int(self.models[-1][j].predict(value))))

    def lower_bound_search(self, sorted_list, value):
        pred = self._predict(value)
        return lower_bound_exponential_search(sorted_list, value, start=pred)


two = MultiLayerRMI(data, [1, 1000])
assert two.lower_bound_search(data, v) == lower_bound_binary_search(data, v)
print("两层 RMI 自检通过")
```

### 步骤 7：带误差界的两层 RMI

```python
class MultiLayerRMIWithBounds(MultiLayerRMI):
    """为每个末层模型记录训练期最大绝对误差，查询时在 [pred-err, pred+err] 内二分。"""

    def __init__(self, sorted_list, layer_sizes):
        super().__init__(sorted_list, layer_sizes)
        self.errors = [0] * len(self.models[-1])
        for index, key in enumerate(sorted_list):
            j = self._get_last_layer_model(key)
            pred = max(0, min(self.N - 1, int(self.models[-1][j].predict(key))))
            self.errors[j] = max(self.errors[j], abs(pred - index))

    def _predict_with_bounds(self, value):
        j = self._get_last_layer_model(value)
        pred = max(0, min(self.N - 1, int(self.models[-1][j].predict(value))))
        e = self.errors[j]
        return pred, max(0, pred - e), min(self.N, pred + e + 1)

    def lower_bound_search(self, sorted_list, value):
        pred, lo, hi = self._predict_with_bounds(value)
        return lower_bound_binary_search(sorted_list, value, start=pred, lo=lo, hi=hi)


two_b = MultiLayerRMIWithBounds(data, [1, 1000])
assert two_b.lower_bound_search(data, v) == lower_bound_binary_search(data, v)
print("带误差界 RMI 自检通过；末层模型误差上界 =", max(two_b.errors))
```

### 步骤 8：误差随末层模型数的变化

```python
import random as _rnd
import statistics
_rnd.seed(0)
SAMPLE = 20_000
sample_idx = sorted(_rnd.sample(range(N), SAMPLE))

def error_stats(rmi):
    errs = [abs(rmi._predict(data[i]) - i) for i in sample_idx]
    return statistics.mean(errs), max(errs)

rows = []
for M in [1, 10, 100, 500, 1000, 5000, 20000]:
    rmi = MultiLayerRMI(data, [1, M])
    mean_err, max_err = error_stats(rmi)
    rows.append((M, mean_err, max_err))
    print(f"末层模型数 {M:>6}: 平均绝对误差 {mean_err:8.1f}  最大绝对误差 {max_err:6d}")
```

画图（需要 `matplotlib`）：

```python
fig, ax = plt.subplots()
xs = [r[0] for r in rows]
ax.plot(xs, [r[1] for r in rows], "o-", label="平均绝对误差")
ax.plot(xs, [r[2] for r in rows], "s--", label="最大绝对误差")
ax.set_xscale("log")
ax.set_yscale("log")
ax.set_xlabel("末层模型数量")
ax.set_ylabel("预测误差（下标）")
ax.legend()
ax.grid(True, which="both", alpha=0.3)
plt.show()
```

> 说明：为缩短运行时间，误差在 20,000 个随机抽样的键上统计。改成 `sample_idx = range(N)` 即为全量统计，结论一致但耗时约 10 倍。

### 步骤 9：性能基准测试

```python
from bisect import bisect_left
import time

random.seed(0)
NUM_QUERIES = 20_000
queries = [random.choice(data) for _ in range(NUM_QUERIES)]

def benchmark(name, fn):
    checksum = 0
    t0 = time.perf_counter()
    for q in queries:
        checksum += fn(q)
    t1 = time.perf_counter()
    ns = (t1 - t0) * 1_000_000_000 / len(queries)
    print(f"{name:<28} {ns:12.0f} ns/次   (checksum={checksum})")
    return ns

t_bisect = benchmark("bisect (C 实现)", lambda q: bisect_left(data, q))
t_pybin  = benchmark("二分查找 (纯 Python)", lambda q: lower_bound_binary_search(data, q))
t_1l     = benchmark("单层 RMI + 指数搜索", lambda q: single.lower_bound_search(data, q))
t_2l     = benchmark("两层 RMI (1000 模型)", lambda q: two.lower_bound_search(data, q))
t_2lb    = benchmark("两层 RMI + 误差界", lambda q: two_b.lower_bound_search(data, q))
```

线性查找太慢，只在少量查询上测：

```python
lin_q = queries[::100]                       # 200 个查询
t0 = time.perf_counter()
for q in lin_q:
    lower_bound_linear_search(data, q)
t1 = time.perf_counter()
print(f"线性查找: {(t1 - t0) / len(lin_q) * 1e6:.0f} μs/次")
```

## 六、观测记录表

### 6.1 基本信息

| 项目 | 你的值 |
| --- | --- |
| Python 版本 | |
| CPU 型号 | |
| `N`（数组长度） | |
| 查询条数 | |

### 6.2 预测误差 vs 末层模型数（实测，N = 200,000，线性样条）

| 末层模型数 | 平均绝对误差（下标） | 最大绝对误差（下标） | 你的实测值 |
| --- | --- | ---: | ---: |
| 1 | 33313.5 | 50000 | |
| 10 | 3518.1 | 15811 | |
| 100 | 355.5 | 4995 | |
| 500 | 73.0 | 2225 | |
| 1000 | 36.3 | 1555 | |
| 5000 | 7.1 | 650 | |
| 20000 | 1.7 | 248 | |

**两个"最大误差"口径**：本表按 20,000 个抽样键统计（`random.Random(0).sample`）；步骤 7 的 `errors` 是对全量 200,000 个键统计的，M = 1000 时为 1558——与本表的 1555 相差 3，属抽样与全量的口径差，不是错误。

观察结论（请在横线上写一句）：

- 模型数每增加 10 倍，平均绝对误差大约降到原来的 ______。
- 平均绝对误差与最大绝对误差的比值随模型数 ______（增大 / 减小 / 基本不变）。

### 6.3 查询耗时（实测，N = 200,000，随机 20,000 次查询）

| 算法 | 实测 ns/次 | 你的实测值 | 相对 bisect 的倍数 |
| --- | ---: | ---: | ---: |
| `bisect` 二分（C 实现） | 523 | | 1.00 |
| 单层 RMI + 指数搜索 | 5798 | | |
| 两层 RMI（1000 模型） | 1817 | | |
| 两层 RMI + 误差界 | 1672 | | |
| 线性查找 | 3,768,000（=3768 μs） | | |

填写要求：把最后一列算出来并保留两位小数。

（耗时类数字：绝对值随机器不同而变化，请以自己的实测为准；趋势——单层最慢、两层 + 误差界最快、`bisect` 与两层接近——在任何机器上都成立。）

### 6.4 现象判断题

在 6.3 表里，哪一个数字最可能让你意外？用一句话写下原因：

```
________________________________________________________________
```

## 七、思考题

1. 给出 `N = 200,000` 时二分查找的比较次数，以及末层 1000 个模型（最大绝对误差 1555）时带误差界 RMI 在区间内二分的比较次数，并说明两者相差多少次。
2. 表 6.2 中，末层模型数从 1 增加到 20000（约 2×10⁴ 倍），平均绝对误差从 33313.5 降到 1.7（也约 2 万倍）——实测误差与模型数近似成反比。为什么是 1/M，而不是分段线性插值理论给出的更快下降？
3. 单层 RMI（5798 ns）比两层 RMI（1817 ns）还慢。既然两层 RMI 的模型更多、每次查询要走两层，为什么它反而更快？
4. 误差界版本（1672 ns）比无误差界版本（1817 ns）快。除了"搜索区间更小"，还有一个更本质的原因是什么？（提示：两种修正方式的算法不同。）
5. 把数据分布从"多项式 + 噪声"改成"完全均匀"，预测误差会怎样变化？此时 RMI 相对二分还有优势吗？
6. 带误差界的 RMI 在查询一个**不在数据集中**的键时，可能返回错误下标。请构造一个能触发该错误的具体例子，并给出修复方案。
7. 为什么 RMI 更擅长静态数据？如果允许插入和删除，会发生什么？
8. 本实验中 RMI 的"缓存友好性"为什么测不出来？设计一个能体现出该优势的实验条件。

## 八、思考题答案与解析

### 8.1 题 1

二分查找比较次数：

```
k_bin = ⌈log₂(N + 1)⌉ = ⌈log₂(200001)⌉
log₂(200001) = ln(200001)/ln(2) ≈ 12.2061/0.6931 ≈ 17.61
k_bin = 18
```

带误差界 RMI，最大绝对误差 `err = 1555`，区间宽度 `2·1555 + 1 = 3111`：

```
k_rmi = ⌈log₂(3111)⌉ = ⌈11.60⌉ = 12
```

相差 `18 − 12 = 6` 次比较。

但要注意：RMI 还要额外付出"两次模型预测 + 属性查找 + clamp"的代价。在纯 Python 里这部分远贵于 12 次整数比较，所以比较次数少并不等于墙钟时间少。

### 8.2 题 2

实测规律：模型数每 ×10，平均绝对误差约 ÷10（33313.5 → 3518.1 → 355.5；1000 → 5000 → 20000 的 36.3 → 7.1 → 1.7 也吻合），即误差 ∝ 1/M。原因有三层：

1. **误差由"段长"主导**。本数据的键序（排序后的键-位置曲线）是一条 √ 形的凸曲线加局部不规则，等深分段的段长 = N/M，每段的拟合偏差与段长近似成正比，所以误差 ∝ 1/M。
2. **数据不是理想光滑的**。±200 的随机噪声与"键域上端稀疏"让误差里始终有一层与段长同阶的不规则贡献。若数据完全光滑（噪声为 0），分段线性插值误差应按段长的平方下降——比实测快一个 M 的幂次。
3. **误差降到底后，加模型只剩开销**。M = 20000 时平均误差 1.7，已接近排序噪声决定的下界；再增模型，预测与参数存储的开销开始占上风。

可以用代码验证第 2 点：把噪声幅度改为 0 再跑一遍步骤 8，误差随 M 的下降会明显快于 1/M。

### 8.3 题 3

关键不在"模型数量"，而在**指数搜索的跨度 `d`**。

指数搜索的代价约为 `2⌈log₂ d⌉ + 1`，其中 `d = |pred − true|`，即预测误差。

- 单层 RMI：平均绝对误差 33313.5 → `d ≈ 33314`，`⌈log₂ 33314⌉ = 16`，约 `2×16+1 = 33` 次比较，且全是 Python 层的循环迭代。
- 两层 RMI（1000 模型）：平均绝对误差 36.3 → `d ≈ 37`，`⌈log₂ 37⌉ = 6`，约 13 次比较。

两层 RMI 多付出的代价只是"一次额外的模型预测 + 一次整数除法取模编号"（O(1) 的浮点运算），省下的却是指数搜索里约 20 次 Python 循环迭代。后者贵得多，所以净收益为正。

对照 6.3 表：`5798 → 1817`，节省约 69%（耗时类数字随机器变化，看量级与方向）。

### 8.4 题 4

更本质的原因是**修正算法换了**：

| 版本 | 修正算法 | 代价构成 |
| --- | --- | --- |
| 无误差界 | 指数搜索 | 外层"扩区间"的 Python while 循环 + 内层二分 |
| 有误差界 | 一次定区间的二分 | 只有二分 |

指数搜索的"扩区间"阶段是一个 Python 层的 `while` 循环，每次迭代要做列表下标访问、比较、`bound *= 2`。这个循环本身就比一次直接算出 `lo/hi` 再二分要贵。误差界版本直接给出 `[pred−err, pred+err]`，省掉了整个扩区间阶段。

所以收益 = （区间变小）+（去掉一层 Python 循环）。后者通常是主要贡献。

验证方法：把误差界人为放大 100 倍（`e = self.errors[j] * 100`），此时区间不再小，但仍然是纯二分。若耗时仍明显低于无误差界版本（约 1817 ns），就说明"去掉扩区间循环"是主因。

### 8.5 题 5

数据均匀时：

- 预测误差会大幅下降。均匀分布的 `F(x)` 本身就是一条直线（严格线性），单层 RMI 的一条直线就能几乎完美拟合，误差接近 0（仅剩离散化误差，量级为 1）。
- 此时"加模型"几乎不再降低误差，6.2 表会变成一条接近水平的线。

RMI 是否还有优势：

- 在**比较次数**上会有优势：`d ≈ 1`，指数搜索约 3 次比较，远少于二分的 18 次。
- 在**墙钟时间**上（本实验的纯 Python 环境）不会有优势：一次 Python 函数调用的开销大于十几次整数比较。
- 在**工程价值**上：均匀数据恰恰是 B 树最容易处理的情况，且真实系统很少是均匀分布的。RMI 的意义在于处理"有规律但非均匀"的分布，例如按时间戳、按 ID 前缀、按地理编码排序的数据。

可以用代码验证：把步骤 1 的 `poly_fn` 改成 `lambda x: 3 * x`，重跑步骤 8。

### 8.6 题 6

构造：令末层某模型 `j` 在训练样本上的最大误差为 `err_j = 5`。查询一个落在模型 `j` 负责区间内、但不在数据集中的键 `x_missing`。此时 `x_missing` 的真实下界位置与 `pred` 的差可能大于 5。

具体例子：

```python
# 取数据集中相邻两个键之间插入一个不存在的键
i = 100000
x_missing = data[i] + 1
while x_missing == data[i + 1]:
    x_missing += 1
print("不存在的键:", x_missing)
print("带误差界 RMI 结果:", two_b.lower_bound_search(data, x_missing))
print("二分查找结果    :", lower_bound_binary_search(data, x_missing))
```

两个结果不一致即触发了该错误。

修复方案（任选）：

1. **边界检查**：若 `lo == 0` 或 `hi == N`，说明区间贴到了数组端点，退化到整表二分。
2. **端点哨兵检查**：检查 `a[lo-1] < value` 与 `a[hi] >= value` 是否成立，不成立则扩大区间重试。
3. **保守放大误差界**：用 `err_j` 时乘以一个安全系数（例如 2），并为未命中情形预留回退路径。真实实现通常结合 1、2 两条。

### 8.7 题 7

RMI 的模型是在**全部数据的秩（rank）**上训练的。一个键的秩等于"排在它前面的元素个数"。

- 插入一个新键后，所有比它大的键的秩都 +1。这意味着**一次插入可能让所有模型的预测整体偏移**，误差界随之失效。
- 删除同理，比被删键大的元素秩都 −1。
- 分布漂移（新数据的键域与训练时不同）也会放大误差：`err_j` 是在旧样本上测的，新样本的真实误差可能远超它。

工程上的应对：

- 静态数据：一次性训练，误差界可靠。
- 需要更新：采用 PG M-index（arXiv:1910.06169）这类带对数级误差保证、支持增删的结构，或 ALEX（arXiv:1905.08898）这类把模型做成可分裂节点的自适应结构。
- 或者定期重训练，并在两次重训练之间监控误差，超阈值即触发重建。

### 8.8 题 8

测不出来的原因：本实验的数据是 **Python 的 `list`**，每个元素是一个 Python 对象（整数对象），访问 `data[mid]` 的成本是"指针解引用 + 对象头开销"，而 Python 解释器本身的字节码分派开销比这个内存访问开销还大。也就是说，cache miss 的时间被解释器开销"稀释"了。

能体现缓存优势的实验条件：

1. 用 `array` 模块或 `numpy` 的紧凑数组（连续内存、每个元素 4/8 字节）替代 Python list，去掉对象开销。
2. 数据量增大到远超 L2/L3 缓存（例如千万级、每元素 8 字节 → 80 MB）。
3. 用 C/C++/Rust 实现，或至少用 `numpy` 的向量化操作完成内层循环。
4. 对比指标改为 cache miss 次数（Linux 下用 `perf stat -e cache-misses`，Windows 下需自行实测确认可用工具）。

在这样的条件下，二分查找前 5~6 次比较每次都会 miss（跨度从 N/2 到 N/64，远超缓存行局部性），而 RMI 只在预测位置附近的几十个字节内访问，几乎全部命中。

## 九、扩展任务（选做）

### 9.1 三层及以上的 RMI

```python
def compute_layer_sizes(num_layers: int, num_models: int) -> list[int]:
    """把 num_models 个模型分到 num_layers 层，除首层外逐层大致翻倍。"""
    sizes = [2 * i + 1 for i in range(num_layers)]      # 最小可行解 [1,3,5,...]
    need = sum(sizes)
    if need > num_models:
        return sizes
    rest = num_models - need
    j = num_layers - 1
    while rest > j:
        sizes[j] += j
        rest -= j
        j = (j - 1) % num_layers
    sizes[j] += rest
    return sizes


print(compute_layer_sizes(3, 20))    # 例：[1, 6, 13] 一类的结果
print(compute_layer_sizes(4, 20))
```

按同样思路构造 `MultiLayerRMI(data, compute_layer_sizes(3, 1000))`，比较"相同总模型数、不同层数"下的误差与耗时。通常层数增加会降低误差，但预测本身的路径变长，二者存在折中。

### 9.2 换成其他模型

把 `LinearSplineModel` 换成"段内最小二乘拟合"的线性回归（用 `numpy.polyfit(x, y, 1)`），比较误差。注意训练成本会从 O(1) 变成 O(段长)。

### 9.3 对照 B 树

自己实现一个最简单的 B 树（或直接用实验材料中 B 树相关 notebook 的实现），在同一数据集上测量 ns/次，与 RMI 对比。

### 9.4 误差界的分布

画出 `two_b.errors` 的直方图，观察 1000 个末层模型的误差是否均匀。若不均匀，思考：哪些键域区间的模型误差最大？为什么？

### 9.5 阅读

精读 Kraska et al., *The Case for Learned Index Structures*, SIGMOD 2018（arXiv:1712.01208），重点看第 3 节 RMI 的 staged model 设计与第 4 节关于字符串索引的讨论。

## 十、常见错误与排查

| 现象 | 原因 | 处理 |
| --- | --- | --- |
| `AssertionError: Inputs must be sorted` | 传给模型的数据未排序，或在构造 RMI 后又修改了 `data` | RMI 构造前先 `sorted()`；构造后不要再改动数组 |
| `ZeroDivisionError` 于 `LinearSplineModel.__init__` | 某段内所有键相同，`X[-1] - X[0] == 0` | 在 `__init__` 中特判 `X[-1] == X[0]`，令斜率为 0 |
| 预测位置越界（`IndexError`） | `predict` 返回值超出 `[0, N-1]` | 用 `clamp`：`max(0, min(N-1, int(pred)))` |
| `IndexError: list index out of range` 于下层模型分配 | `int(pred / N * M)` 在 `pred == N` 时等于 `M`，越界 | 分配与查询时都做 `min(nxt, M - 1)` |
| 带误差界版本结果与二分不一致 | 查询了数据集中不存在的键（见题 6） | 只对存在的键查询，或按题 6 方案加边界检查与回退 |
| 构造 20000 个末层模型时内存暴涨 | 训练期 `model_keys` 保存了每个模型的键列表副本，总量约 `2N` 个整数 | 减小 `N`，或改用"记录每个模型的键区间 [start, end)"代替拷贝键列表 |
| 步骤 8 跑得极慢 | 对每个配置都在全量 `N` 个键上做预测 | 用 20,000 的抽样（`sample_idx`）先跑通，再改全量 |
| 线性查找把程序卡住 | 查询条数太多，每条 O(N) | 用 `queries[::100]` 抽样，且不要在 `N` 很大时全量跑 |
| `bisect` 与手写二分结果不同 | 两者定义一致（都是下界），若不同说明手写实现有误 | 用第 2 步的对拍代码逐值校验 |
| 每次运行的 ns/次 差异很大 | 后台进程干扰、`time.time()` 精度不足 | 用 `time.perf_counter()`；关闭其他程序；多跑几次取中位数 |
| 误以为"RMI 一定比二分快" | 见 3.4 节 | 明确区分"比较次数"与"墙钟时间"，明确实现语言的影响 |

## 附：本实验引用的文献

| 编号 | 出处 | 本实验用到什么 |
| --- | --- | --- |
| [1] | Kraska et al., *The Case for Learned Index Structures*, SIGMOD 2018（arXiv:1712.01208） | CDF 类比、`pos ≈ F(x)·N`、RMI 分层结构 |
| [2] | Ferragina & Vinciguerra, PG M-index（arXiv:1910.06169） | 带误差保证、支持更新的学习索引（扩展阅读） |
| [3] | Ding et al., ALEX（arXiv:1905.08898） | 自适应、可更新的学习索引（扩展阅读） |

局限性归纳（出自 [1] 及其后续工作）：早期 RMI 更适合**静态、内存内、有序、一维**的数据；插入删除会改变秩；分布漂移会放大误差。