Sklearn-源码解析-书-v1-0-九-

Sklearn 源码解析（书）v1.0（九）

逐行注释：

1. 定义类型：y、w 为连续浮点数组，target 为 intp 数组初始化为 [0, 1, ..., n-1]

2. target 编码块结构：块 [i..j] 满足 target[i]=j 且 target[j]=i

3. 活跃块起始索引 i 存储块聚合值：w[i] 为权重和，y[i] 为加权平均

4. 主循环：i 从 0 遍历到 n-1，k = target[i]+1 为下一块起始

5. 若 y[i] < y[k] 单调递增，前进 i = k 继续

6. 否则进入违规合并：累积 sum_wy、sum_w，跳转 k = target[k]+1 继续合并

7. 直到 k==n 或 prev_y < y[k]（违规序列结束）

8. 更新块起始 i 的聚合值：y[i]=sum_wy/sum_w，w[i]=sum_w，更新 target 指针

9. 回溯关键：i = target[i-1] 回退到前一块末尾，保证单次遍历 O(n)

10. 重构解：遍历 target 链表，将块内所有值设为块均值 y[i]

代码总结：这段代码实现了 PAVA 算法的 Cython 优化版本。target 数组巧妙编码了块的双向链表结构，回溯机制 i = target[i-1] 确保每个元素最常被访问常数次，实现了真正的 O(n) 单次遍历。with nogil 释放 GIL 实现并行友好，原地修改 y、w 节省内存。这是保序回归高性能的基石。

19.8.1.2 Cython 去重实现

# 第 19 章 —— sklearn/_isotonic.pyx - _make_unique (第62-110行)
def _make_unique(const floating[::1] X,
                 const floating[::1] y,
                 const floating[::1] sample_weights):
    """Average targets for duplicate X, drop duplicates."""
    unique_values = len(np.unique(X))

    if floating is float:
        dtype = np.float32
    else:
        dtype = np.float64

    cdef floating[::1] y_out = np.empty(unique_values, dtype=dtype)
    cdef floating[::1] x_out = np.empty_like(y_out)
    cdef floating[::1] weights_out = np.empty_like(y_out)

    cdef floating current_x = X[0]
    cdef floating current_y = 0
    cdef floating current_weight = 0
    cdef int i = 0
    cdef int j
    cdef floating x
    cdef int n_samples = len(X)
    cdef floating eps = np.finfo(dtype).resolution

    for j in range(n_samples):
        x = X[j]
        if x - current_x >= eps:
            # next unique value
            x_out[i] = current_x
            weights_out[i] = current_weight
            y_out[i] = current_y / current_weight
            i += 1
            current_x = x
            current_weight = sample_weights[j]
            current_y = y[j] * sample_weights[j]
        else:
            current_weight += sample_weights[j]
            current_y += y[j] * sample_weights[j]

    x_out[i] = current_x
    weights_out[i] = current_weight
    y_out[i] = current_y / current_weight
    return(
        np.asarray(x_out[:i+1]),
        np.asarray(y_out[:i+1]),
        np.asarray(weights_out[:i+1]),
    )

逐行注释：

1. 计算唯一值数量预分配数组，根据 floating 类型模板选择 float32/float64

2. 获取浮点数分辨率 eps = np.finfo(dtype).resolution 作为去重容差

3. 遍历已排序的 X，若 x - current_x >= eps 视为新值，输出当前聚合块

4. 否则累积权重 current_weight 与加权目标 current_y

5. 循环结束输出最后一块，返回切片后的数组

代码总结：这段代码实现了 Cython 级别的重复 X 值合并。使用 np.finfo(dtype).resolution 作为容差，体现了浮点数比较的工程智慧：float32 分辨率约 1e-7，float64 约 1e-16，自动适配精度。加权平均 y_out = sum(w*y)/sum(w) 保证了样本权重正确传递。模板化设计避免了为每种 dtype 重复代码。

---

19.9 IsotonicRegression 完整流程 —— 从拟合到插值预测

IsotonicRegression 将 PAVA 算法封装为完整的 scikit-learn 估计器。fit 阶段验证输入、自动判断单调方向、排序去重、PAVA 拟合、trim_duplicates 精简插值点。_build_f 构建 interp1d 线性插值函数，处理越界模式。序列化通过存储阈值数组重建插值函数。

flowchart TD subgraph IsotonicFit["IsotonicRegression.fit() 流程"] direction TB validate["validate_data: 检查输入形状(1D/单特征2D)\ncheck_array dtype=float32/64"] check_inc["increasing='auto'?\n是: check_increasing(Spearman+Fisher)\n否: 使用用户指定"] filter["过滤零权重样本\nmask = sample_weight > 0"] sort["lexsort((y, X)): 先按X后按y排序\n保证相同X下y有序"] unique["_make_unique: 合并重复X\n加权平均y, resolution容差"] pava["isotonic_regression: PAVA求解\nSciPy>=1.12用优化器,否则Cython"] bounds["记录 X_min_, X_max_ 训练域边界"] trim["trim_duplicates=True:\n移除y等于前后邻居的冗余点\n保留首尾,加速interp1d"] build_f["_build_f: 构建interp1d线性插值\nbounds_error='raise'模式\n单点退化为常数函数"] store["存储 X_thresholds_, y_thresholds_\n用于pickle序列化重建"] end validate --> check_inc check_inc --> filter filter --> sort sort --> unique unique --> pava pava --> bounds bounds --> trim trim --> build_f build_f --> store

19.9.1 逐行解析关键函数

19.9.1.1 核心拟合流程

# 第 19 章 —— sklearn/isotonic.py - IsotonicRegression._build_y (第250-300行)
def _build_y(self, X, y, sample_weight, trim_duplicates=True):
    """Build the y_ IsotonicRegression."""
    self._check_input_data_shape(X)
    X = X.reshape(-1)  # use 1d view

    # Determine increasing if auto-determination requested
    if self.increasing == "auto":
        self.increasing_ = check_increasing(X, y)
    else:
        self.increasing_ = self.increasing

    # If sample_weights is passed, removed zero-weight values and clean
    # order
    sample_weight = _check_sample_weight(sample_weight, X, dtype=X.dtype)
    mask = sample_weight > 0
    X, y, sample_weight = X[mask], y[mask], sample_weight[mask]

    order = np.lexsort((y, X))
    X, y, sample_weight = [array[order] for array in [X, y, sample_weight]]
    unique_X, unique_y, unique_sample_weight = _make_unique(X, y, sample_weight)

    X = unique_X
    y = isotonic_regression(
        unique_y,
        sample_weight=unique_sample_weight,
        y_min=self.y_min,
        y_max=self.y_max,
        increasing=self.increasing_,
    )

    # Handle the left and right bounds on X
    self.X_min_, self.X_max_ = np.min(X), np.max(X)

    if trim_duplicates:
        # Remove unnecessary points for faster prediction
        keep_data = np.ones((len(y),), dtype=bool)
        # Aside from the 1st and last point, remove points whose y values
        # are equal to both the point before and the point after it.
        keep_data[1:-1] = np.logical_or(
            np.not_equal(y[1:-1], y[:-2]), np.not_equal(y[1:-1], y[2:])
        )
        return X[keep_data], y[keep_data]
    else:
        return X, y

逐行注释：

1. 验证输入形状（1D 或单特征 2D），展平为 1D

2. 自动单调方向：increasing='auto' 调用 check_increasing 基于 Spearman 相关判断

3. 过滤零权重样本，按 (X, y) 字典序排序（lexsort 先按 X 后按 y）

4. 调用 _make_unique 合并重复 X 值（加权平均 y）

5. 调用 isotonic_regression 统一入口（SciPy>=1.12 用优化器，旧版回退 Cython PAVA）

6. 记录训练域边界 X_min_、X_max_

7. trim_duplicates 优化：移除 y 值等于前后邻居的冗余点（保留首尾），加速插值预测

8. 返回精简后的阈值数组

代码总结：这段代码实现了保序回归的核心拟合流水线。关键设计包括：lexsort((y, X)) 确保相同 X 下 y 有序，配合 _make_unique 正确处理并列值；trim_duplicates 移除线性插值中冗余的共线点（y 等于前后邻居），显著减少 interp1d 的断点数量，加速预测且不改变插值结果。自动单调方向检测降低了用户使用门槛。

19.9.1.2 统一入口函数

# 第 19 章 —— sklearn/isotonic.py - isotonic_regression (第112-155行)
@validate_params(
    {
        "y": ["array-like"],
        "sample_weight": ["array-like", None],
        "y_min": [Interval(Real, None, None, closed="both"), None],
        "y_max": [Interval(Real, None, None, closed="both"), None],
        "increasing": ["boolean"],
    },
    prefer_skip_nested_validation=True,
)
def isotonic_regression(
    y, *, sample_weight=None, y_min=None, y_max=None, increasing=True
):
    """Solve the isotonic regression model."""
    y = check_array(y, ensure_2d=False, input_name="y", dtype=[np.float64, np.float32])
    if sp_base_version >= parse_version("1.12.0"):
        res = optimize.isotonic_regression(
            y=y, weights=sample_weight, increasing=increasing
        )
        y = np.asarray(res.x, dtype=y.dtype)
    else:
        # TODO: remove this branch when Scipy 1.12 is the minimum supported version
        # Also remove _inplace_contiguous_isotonic_regression.
        order = np.s_[:] if increasing else np.s_[::-1]
        y = np.array(y[order], dtype=y.dtype)
        sample_weight = _check_sample_weight(sample_weight, y, dtype=y.dtype, copy=True)
        sample_weight = np.ascontiguousarray(sample_weight[order])
        _inplace_contiguous_isotonic_regression(y, sample_weight)
        y = y[order]

    if y_min is not None or y_max is not None:
        # Older versions of np.clip don't accept None as a bound, so use np.inf
        if y_min is None:
            y_min = -np.inf
        if y_max is None:
            y_max = np.inf
        np.clip(y, y_min, y_max, y)
    return y

逐行注释：

1. 参数验证：y 为 1D 数组，dtype 限制 float32/64

2. SciPy >= 1.12：调用 scipy.optimize.isotonic_regression 优化实现

3. 旧版回退：按 increasing 方向切片，调用 Cython PAVA 核心 _inplace_contiguous_isotonic_regression

4. 结果按原顺序恢复

5. 应用 y_min/y_max 边界裁剪

6. 返回单调拟合后的 y

代码总结：这段代码实现了 isotonic_regression 统一入口函数，体现了 scikit-learn "新库优先、旧版兼容" 的版本适配策略。SciPy 1.12+ 提供了原生优化的 PAVA 实现，旧版回退到自研 Cython 版本。该函数被 IsotonicRegression._build_y 调用，也是用户可直接调用的函数式接口。

19.9.1.3 插值函数构建

# 第 19 章 —— sklearn/isotonic.py - IsotonicRegression._build_f (第230-250行)
def _build_f(self, X, y):
    """Build the f_ interp1d function."""

    bounds_error = self.out_of_bounds == "raise"
    if len(y) == 1:
        # single y, constant prediction
        self.f_ = lambda x: y.repeat(x.shape)
    else:
        self.f_ = interpolate.interp1d(
            X, y, kind="linear", bounds_error=bounds_error
        )

逐行注释：

1. bounds_error=True 对应 out_of_bounds='raise' 抛出异常

2. 单点特殊处理：返回常数函数 lambda x: y.repeat(x.shape)

3. 多点情况：使用 scipy.interpolate.interp1d 线性插值

代码总结：这段代码构建了预测用的插值函数。interp1d 的 bounds_error 参数直接映射 'raise' 模式，'clip' 和 'nan' 模式在 _transform 中通过 np.clip 和 NaN 返回处理。单点退化为常数函数避免 interp1d 报错。

19.9.1.4 统一预测/变换逻辑

# 第 19 章 —— sklearn/isotonic.py - IsotonicRegression._transform (第350-380行)
def _transform(self, T):
    """`_transform` is called by both `transform` and `predict` methods.

    Since `transform` is wrapped to output arrays of specific types (e.g.
    NumPy arrays, pandas DataFrame), we cannot make `predict` call `transform`
    directly.

    The above behaviour could be changed in the future, if we decide to output
    other type of arrays when calling `predict`.
    """
    if hasattr(self, "X_thresholds_"):
        dtype = self.X_thresholds_.dtype
    else:
        dtype = np.float64

    T = check_array(T, dtype=dtype, ensure_2d=False)

    self._check_input_data_shape(T)
    T = T.reshape(-1)  # use 1d view

    if self.out_of_bounds == "clip":
        T = np.clip(T, self.X_min_, self.X_max_)

    res = self.f_(T)

    # on scipy 0.17, interp1d up-casts to float64, so we cast back
    res = res.astype(T.dtype)

    return res

逐行注释：

1. 获取训练数据 dtype，默认 float64

2. 验证输入形状，展平为 1D

3. out_of_bounds='clip'：np.clip 限制在训练域 [X_min_, X_max_]

4. 调用插值函数 f_ 计算预测值

5. 结果转回输入 dtype（应对旧版 interp1d 向上转型 float64）

代码总结：这段代码实现了 predict 和 transform 共享的核心预测逻辑。out_of_bounds 三种模式：'raise' 由 interp1d 内部 bounds_error 抛出异常，'clip' 由 np.clip 预处理输入，'nan' 由 interp1d 默认行为返回 NaN。dtype 保持一致性体现了对用户数据类型的尊重。

19.9.1.5 公共预测接口

# 第 19 章 —— sklearn/isotonic.py - IsotonicRegression.predict / transform (第382-410行)
def transform(self, T):
    """Transform new data by linear interpolation."""
    return self._transform(T)

def predict(self, T):
    """Predict new data by linear interpolation."""
    return self._transform(T)

逐行注释：

1. transform 调用 _transform，返回插值结果

2. predict 调用 _transform，返回相同结果（回归任务预测即变换）

代码总结：这段代码体现了 RegressorMixin 和 TransformerMixin 的双重身份：predict 用于回归预测，transform 用于流水线特征变换，底层共享 _transform 逻辑。设计说明中提到未来可能区分两者输出类型（如 transform 输出 DataFrame），当前保持一致。

19.9.1.6 输入形状验证

# 第 19 章 —— sklearn/isotonic.py - IsotonicRegression._check_input_data_shape (第215-225行)
def _check_input_data_shape(self, X):
    if not (X.ndim == 1 or (X.ndim == 2 and X.shape[1] == 1)):
        msg = (
            "Isotonic regression input X should be a 1d array or "
            "2d array with 1 feature"
        )
        raise ValueError(msg)

逐行注释：

1. 检查 X.ndim 为 1 或 (2, 1)，即 1D 数组或单特征 2D 数组

2. 不满足条件抛出 ValueError

代码总结：这段代码严格限制输入形状，保证 IsotonicRegression 仅处理一维特征。fit、predict、transform 均调用此验证，确保一致性。

19.9.1.7 序列化支持

# 第 19 章 —— sklearn/isotonic.py - IsotonicRegression.__getstate__ (第410-420行)
def __getstate__(self):
    """Pickle-protocol - return state of the estimator."""
    state = super().__getstate__()
    # remove interpolation method
    state.pop("f_", None)
    return state

# 第 19 章 —— sklearn/isotonic.py - IsotonicRegression.__setstate__ (第420-430行)
def __setstate__(self, state):
    """Pickle-protocol - set state of the estimator.

    We need to rebuild the interpolation function.
    """
    super().__setstate__(state)
    if hasattr(self, "X_thresholds_") and hasattr(self, "y_thresholds_"):
        self._build_f(self.X_thresholds_, self.y_thresholds_)

逐行注释：

1. getstate：调用父类方法获取状态字典，移除不可 pickle 的 f_ (interp1d 对象)

2. setstate：恢复父类状态，若存在阈值数组则调用 _build_f 重建插值函数

代码总结：这段代码实现了保序回归模型的序列化支持。scipy.interpolate.interp1d 对象包含复杂内部状态无法直接 pickle，设计选择存储关键阈值数组 X_thresholds_、y_thresholds_，反序列化时重建插值函数。这种"存参数、重建函数"模式是 scikit-learn 处理不可序列化科学计算对象的通用范式。

19.9.1.8 Array API 兼容标记

# 第 19 章 —— sklearn/isotonic.py - IsotonicRegression.__sklearn_tags__ (第432-440行)
def __sklearn_tags__(self):
    tags = super().__sklearn_tags__()
    tags.input_tags.one_d_array = True
    tags.input_tags.two_d_array = False
    return tags

逐行注释：

1. 调用父类方法获取标签字典

2. 设置 one_d_array=True、two_d_array=False，标记仅接受 1D 数组输入

3. 返回更新后的标签

代码总结：这段代码声明了 IsotonicRegression 的输入标签，明确仅支持 1D 数组（或单特征 2D，内部展平）。这是 scikit-learn 元数据路由与输入验证体系的关键配置。

---

19.10 单调性自动检测 —— Spearman 相关性与 Fisher 变换

check_increasing 基于 Spearman 秩相关系数自动判断单调方向，通过 Fisher z 变换计算 95% 置信区间，若 CI 跨越零点则发出警告，提示判断可能不可靠。

flowchart TD subgraph CheckIncreasing["check_increasing() 单调方向自动检测"] direction TB spearman["Spearman秩相关系数 rho\n捕捉单调关系(不限线性)"] sign["increasing_bool = rho >= 0\n符号决定单调方向"] fisher["Fisher变换: F = 0.5*log((1+rho)/(1-rho))\n相关系数->近似正态分布"] se["标准误: F_se = 1/sqrt(n-3)"] ci["95% CI: F ± 1.96*F_se\n反变换: tanh() 回相关系数域"] warn["CI跨零? sign(rho_0) != sign(rho_1)\n是: UserWarning 提示判断可能不可靠"] edge_cases["边界情况:\nrho=±1 无标准误跳过CI\nn<=3 自由度不足跳过CI"] end spearman --> sign spearman --> fisher fisher --> se se --> ci ci --> warn warn -.-> edge_cases

19.10.1 逐行解析关键函数

# 第 19 章 —— sklearn/isotonic.py - check_increasing (第35-85行)
def check_increasing(x, y):
    """Determine whether y is monotonically correlated with x."""
    # Calculate Spearman rho estimate and set return accordingly.
    rho, _ = spearmanr(x, y)
    increasing_bool = rho >= 0

    # Run Fisher transform to get the rho CI, but handle rho=+/-1
    if rho not in [-1.0, 1.0] and len(x) > 3:
        F = 0.5 * math.log((1.0 + rho) / (1.0 - rho))
        F_se = 1 / math.sqrt(len(x) - 3)

        # Use a 95% CI, i.e., +/-1.96 S.E.
        rho_0 = math.tanh(F - 1.96 * F_se)
        rho_1 = math.tanh(F + 1.96 * F_se)

        # Warn if the CI spans zero.
        if np.sign(rho_0) != np.sign(rho_1):
            warnings.warn(
                "Confidence interval of the Spearman "
                "correlation coefficient spans zero. "
                "Determination of ``increasing`` may be "
                "suspect."
            )

    return increasing_bool

逐行注释：

1. 计算 Spearman 秩相关系数 rho，符号决定单调方向 increasing_bool = rho >= 0

2. Fisher 变换：F = 0.5 * log((1+rho)/(1-rho)) 将相关系数变换到近似正态分布

3. 标准误 F_se = 1/sqrt(n-3)，95% CI 为 F ± 1.96*F_se

4. 反变换回相关系数域：rho_0 = tanh(F - 1.96F_se)，rho_1 = tanh(F + 1.96F_se)

5. 若 CI 跨越零点（sign(rho_0) != sign(rho_1)），发出 UserWarning

6. 边界情况：rho=±1 无标准误跳过 CI；n<=3 自由度不足跳过 CI

代码总结：这段代码实现了统计学严谨的单调方向自动检测。Spearman 秩相关系数捕捉单调关系（不限线性），Fisher 变换提供置信区间量化不确定性。CI 跨零警告保护用户免受弱相关数据误导。这种"统计检验+置信区间+警告"的设计模式，体现了 scikit-learn 在自动化决策中平衡易用性与统计严谨性的工程哲学。

---

19.11 LDA/QDA 测试矩阵 —— 从玩具数据到统计验证

测试套件覆盖求解器×收缩参数全组合、概率校准数学验证、先验处理、系数一致性、降维正交性、QDA 正则化必要性等核心场景。

19.11.1 核心测试解析

19.11.1.1 求解器全组合测试

# 第 19 章 —— sklearn/tests/test_discriminant_analysis.py - test_lda_predict (第50-95行)
solver_shrinkage = [
    ("svd", None),
    ("lsqr", None),
    ("eigen", None),
    ("lsqr", "auto"),
    ("lsqr", 0),
    ("lsqr", 0.43),
    ("eigen", "auto"),
    ("eigen", 0),
    ("eigen", 0.43),
]

def test_lda_predict():
    # Test LDA classification.
    # This checks that LDA implements fit and predict and returns correct
    # values for simple toy data.
    for test_case in solver_shrinkage:
        solver, shrinkage = test_case
        clf = LinearDiscriminantAnalysis(solver=solver, shrinkage=shrinkage)
        y_pred = clf.fit(X, y).predict(X)
        assert_array_equal(y_pred, y, "solver %s" % solver)

        # Assert that it works with 1D data
        y_pred1 = clf.fit(X1, y).predict(X1)
        assert_array_equal(y_pred1, y, "solver %s" % solver)

        # Test probability estimates
        y_proba_pred1 = clf.predict_proba(X1)
        assert_array_equal((y_proba_pred1[:, 1] > 0.5) + 1, y, "solver %s" % solver)
        y_log_proba_pred1 = clf.predict_log_proba(X1)
        assert_allclose(
            np.exp(y_log_proba_pred1),
            y_proba_pred1,
            rtol=1e-6,
            atol=1e-6,
            err_msg="solver %s" % solver,
        )

        # Primarily test for commit 2f34950 -- "reuse" of priors
        y_pred3 = clf.fit(X, y3).predict(X)
        # LDA shouldn't be able to separate those
        assert np.any(y_pred3 != y3), "solver %s" % solver

    clf = LinearDiscriminantAnalysis(solver="svd", shrinkage="auto")
    with pytest.raises(NotImplementedError):
        clf.fit(X, y)

    clf = LinearDiscriminantAnalysis(
        solver="lsqr", shrinkage=0.1, covariance_estimator=ShrunkCovariance()
    )
    with pytest.raises(
        ValueError,
        match=(
            "covariance_estimator and shrinkage "
            "parameters are not None. "
            "Only one of the two can be set."
        ),
    ):
        clf.fit(X, y)

    # test bad solver with covariance_estimator
    clf = LinearDiscriminantAnalysis(solver="svd", covariance_estimator=LedoitWolf())
    with pytest.raises(
        ValueError, match="covariance estimator is not supported with svd"
    ):
        clf.fit(X, y)

    # test bad covariance estimator
    clf = LinearDiscriminantAnalysis(
        solver="lsqr", covariance_estimator=KMeans(n_clusters=2, n_init="auto")
    )
    with pytest.raises(ValueError):
        clf.fit(X, y)

代码总结：这段代码定义了 9 种 solver×shrinkage 组合的参数化测试矩阵，覆盖了所有合法参数组合。测试验证 fit/predict/predict_proba/predict_log_proba 在玩具数据上的正确性，包括 1D 数据、二分类/多分类、概率归一化一致性。svd 不支持 shrinkage 的组合预期抛出 NotImplementedError。互斥参数组合（shrinkage 与 covariance_estimator 同时设置）预期抛出 ValueError。错误的协方差估计器（无 covariance_ 属性）也被捕获。

19.11.1.2 概率校准数学验证

# 第 19 章 —— sklearn/tests/test_discriminant_analysis.py - test_lda_predict_proba (第97-160行)
@pytest.mark.parametrize("n_classes", [2, 3])
@pytest.mark.parametrize("solver", ["svd", "lsqr", "eigen"])
def test_lda_predict_proba(solver, n_classes):
    def generate_dataset(n_samples, centers, covariances, random_state=None):
        """Generate a multivariate normal data given some centers and
        covariances"""
        rng = check_random_state(random_state)
        X = np.vstack(
            [
                rng.multivariate_normal(mean, cov, size=n_samples // len(centers))
                for mean, cov in zip(centers, covariances)
            ]
        )
        y = np.hstack(
            [[clazz] * (n_samples // len(centers)) for clazz in range(len(centers))]
        )
        return X, y

    blob_centers = np.array([[0, 0], [-10, 40], [-30, 30]])[:n_classes]
    blob_stds = np.array([[[10, 10], [10, 100]]] * len(blob_centers))
    X, y = generate_dataset(
        n_samples=90000, centers=blob_centers, covariances=blob_stds, random_state=42
    )
    lda = LinearDiscriminantAnalysis(
        solver=solver, store_covariance=True, shrinkage=None
    ).fit(X, y)
    # check that the empirical means and covariances are close enough to the
    # one used to generate the data
    assert_allclose(lda.means_, blob_centers, atol=1e-1)
    assert_allclose(lda.covariance_, blob_stds[0], atol=1)

    # implement the method to compute the probability given in The Elements
    # of Statistical Learning (cf. p.127, Sect. 4.4.5 "Logistic Regression
    # or LDA?")
    precision = linalg.inv(blob_stds[0])
    alpha_k = []
    alpha_k_0 = []
    for clazz in range(len(blob_centers) - 1):
        alpha_k.append(
            np.dot(precision, (blob_centers[clazz] - blob_centers[-1])[:, np.newaxis])
        )
        alpha_k_0.append(
            np.dot(
                -0.5 * (blob_centers[clazz] + blob_centers[-1])[np.newaxis, :],
                alpha_k[-1],
            )
        )

    sample = np.array([[-22, 22]])

    def discriminant_func(sample, coef, intercept, clazz):
        return np.exp(intercept[clazz] + np.dot(sample, coef[clazz])).item()

    prob = np.array(
        [
            float(
                discriminant_func(sample, alpha_k, alpha_k_0, clazz)
                / (
                    1
                    + sum(
                        [
                            discriminant_func(sample, alpha_k, alpha_k_0, clazz)
                            for clazz in range(n_classes - 1)
                        ]
                    )
                )
            )
            for clazz in range(n_classes - 1)
        ]
    )

    prob_ref = 1 - np.sum(prob)

    # check the consistency of the computed probability
    # all probabilities should sum to one
    prob_ref_2 = float(
        1
        / (
            1
            + sum(
                [
                    discriminant_func(sample, alpha_k, alpha_k_0, clazz)
                    for clazz in range(n_classes - 1)
                ]
            )
        )
    )

    assert prob_ref == pytest.approx(prob_ref_2)
    # check that the probability of LDA are close to the theoretical
    # probabilities
    assert_allclose(
        lda.predict_proba(sample), np.hstack([prob, prob_ref])[np.newaxis], atol=1e-2
    )

代码总结：这段代码实现了基于《统计学习基础》(ESL) 书 P.127 公式的理论后验概率计算，生成已知协方差的高斯混合数据，验证 LDA 概率输出与理论值的一致性（容差 1e-2）。这种"生成真实分布数据 → 理论计算 ground truth → 对比模型输出"的测试模式，是验证概率校准正确性的黄金标准。参数化测试覆盖二分类/三分类与三大求解器的全组合。

19.11.1.3 降维正交性验证

# 第 19 章 —— sklearn/tests/test_discriminant_analysis.py - test_lda_orthogonality (第230-270行)
def test_lda_orthogonality():
    # arrange four classes with their means in a kite-shaped pattern
    means = np.array([[0, 0, -1], [0, 2, 0], [0, -2, 0], [0, 0, 5]])

    # We construct perfectly symmetric distributions, so the LDA can estimate
    # precise means.
    scatter = np.array(
        [
            [0.1, 0, 0],
            [-0.1, 0, 0],
            [0, 0.1, 0],
            [0, -0.1, 0],
            [0, 0, 0.1],
            [0, 0, -0.1],
        ]
    )

    X = (means[:, np.newaxis, :] + scatter[np.newaxis, :, :]).reshape((-1, 3))
    y = np.repeat(np.arange(means.shape[0]), scatter.shape[0])

    # Fit LDA and transform the means
    clf = LinearDiscriminantAnalysis(solver="svd").fit(X, y)
    means_transformed = clf.transform(means)

    d1 = means_transformed[3] - means_transformed[0]
    d2 = means_transformed[2] - means_transformed[1]
    d1 /= np.sqrt(np.sum(d1**2))
    d2 /= np.sqrt(np.sum(d2**2))

    # the transformed within-class covariance should be the identity matrix
    assert_almost_equal(np.cov(clf.transform(scatter).T), np.eye(2))

    # the means of classes 0 and 3 should lie on the first component
    assert_almost_equal(np.abs(np.dot(d1[:2], [1, 0])), 1.0)

    # the means of classes 1 and 2 should lie on the second component
    assert_almost_equal(np.abs(np.dot(d2[:2], [0, 1])), 1.0)

代码总结：这段代码设计了风筝形类均值分布，验证 LDA 降维的正交性与方差最大化特性。构造完美球形类内分布（协方差为各向同性），变换后类内协方差应为单位矩阵（白化性质）。长轴类均值差应映射到第一判别成分（最大方差方向），短轴映射到第二成分。这种几何构造测试精准验证了 Fisher 判别准则的数学实现。

19.11.1.4 QDA 正则化必要性测试

# 第 19 章 —— sklearn/tests/test_discriminant_analysis.py - test_qda_regularization (第370-430行)
@pytest.mark.parametrize("solver", ["svd", "eigen"])
def test_qda_regularization(global_random_seed, solver):
    # The default is reg_param=0. and will cause issues when there is a
    # constant variable.
    rng = np.random.default_rng(global_random_seed)

    # Fitting on data with constant variable without regularization
    # triggers a LinAlgError.
    msg = r"The covariance matrix of class .+ is not full rank."
    clf = QuadraticDiscriminantAnalysis(solver=solver)
    with pytest.raises(linalg.LinAlgError, match=msg):
        clf.fit(X2, y6)

    with pytest.raises(AttributeError):
        y_pred = clf.predict(X2)

    # Adding a little regularization fixes the fit time error.
    if solver == "svd":
        clf = QuadraticDiscriminantAnalysis(solver=solver, reg_param=0.01)
    elif solver == "eigen":
        clf = QuadraticDiscriminantAnalysis(solver=solver, shrinkage=0.01)
    with warnings.catch_warnings():
        warnings.simplefilter("error")
    clf.fit(X2, y6)
    y_pred = clf.predict(X2)
    assert_array_equal(y_pred, y6)

    # LinAlgError should also be there for the n_samples_in_a_class <
    # n_features case.
    X = rng.normal(size=(9, 4))
    y = np.array([1, 1, 1, 1, 1, 1, 2, 2, 2])

    clf = QuadraticDiscriminantAnalysis(solver=solver)
    if solver == "svd":
        msg2 = msg + " When using `solver='svd'`"
    elif solver == "eigen":
        msg2 = msg

    with pytest.raises(linalg.LinAlgError, match=msg2):
        clf.fit(X, y)

    # The error will persist even with regularization for SVD
    # because the number of singular values is limited by n_samples_in_a_class.
    if solver == "svd":
        clf = QuadraticDiscriminantAnalysis(solver=solver, reg_param=0.3)
        with pytest.raises(linalg.LinAlgError, match=msg2):
            clf.fit(X, y)
    # The warning will be gone for Eigen with regularization, because
    # the covariance matrix will be full-rank.
    elif solver == "eigen":
        clf = QuadraticDiscriminantAnalysis(solver=solver, shrinkage=0.3)
        clf.fit(X, y)

代码总结：这段代码系统验证了 QDA 正则化机制的边界。常数特征导致协方差奇异，reg_param/shrinkage 可修复。但 n_samples < n_features 时，SVD 求解器奇异值数量受限于样本数，根本性秩缺失无法通过正则化修复（报错建议改用 eigen）；Eigen 求解器配合 shrinkage 通过收缩向单位矩阵靠拢，可修复秩缺失。测试精准界定了两种求解器的能力边界。

---

19.12 IsotonicRegression 测试全景 —— 从边界条件到性能基准

测试覆盖排列不变性、PAVA 算法对标 R isotone 包、边界处理、序列化、快速预测优化、数据类型一致性等全维度验证。

19.12.1 核心测试解析

19.12.1.1 排列不变性与 PAVA 正确性

# 第 19 章 —— sklearn/tests/test_isotonic.py - test_permutation_invariance (第15-30行)
def test_permutation_invariance():
    # check that fit is permutation invariant.
    # regression test of missing sorting of sample-weights
    ir = IsotonicRegression()
    x = [1, 2, 3, 4, 5, 6, 7]
    y = [1, 41, 51, 1, 2, 5, 24]
    sample_weight = [1, 2, 3, 4, 5, 6, 7]
    x_s, y_s, sample_weight_s = shuffle(x, y, sample_weight, random_state=0)
    y_transformed = ir.fit_transform(x, y, sample_weight=sample_weight)
    y_transformed_s = ir.fit(x_s, y_s, sample_weight=sample_weight_s).transform(x)

    assert_array_equal(y_transformed, y_transformed_s)

代码总结：这段代码验证了拟合过程对样本顺序的不变性。打乱 (X, y, sample_weight) 后拟合，在原始 X 上预测应得到相同结果。这是回归测试，修复了早期版本样本权重未随数据排序导致的 Bug。

19.12.1.2 对标 R isotone 包

# 第 19 章 —— sklearn/tests/test_isotonic.py - test_isotonic_regression_ties_secondary_ (第80-110行)
def test_isotonic_regression_ties_secondary_():
    """
    Test isotonic regression fit, transform  and fit_transform
    against the "secondary" ties method and "pituitary" data from R
     "isotone" package, as detailed in: J. d. Leeuw, K. Hornik, P. Mair,
     Isotone Optimization in R: Pool-Adjacent-Violators Algorithm
    (PAVA) and Active Set Methods

    Set values based on pituitary example and
     the following R command detailed in the paper above:
    > library("isotone")
    > data("pituitary")
    > res1 <- gpava(pituitary$age, pituitary$size, ties="secondary")
    > res1$x

    `isotone` version: 1.0-2, 2014-09-07
    R version: R version 3.1.1 (2014-07-10)
    """
    x = [8, 8, 8, 10, 10, 10, 12, 12, 12, 14, 14]
    y = [21, 23.5, 23, 24, 21, 25, 21.5, 22, 19, 23.5, 25]
    y_true = [
        22.22222,
        22.22222,
        22.22222,
        22.22222,
        22.22222,
        22.22222,
        22.22222,
        22.22222,
        22.22222,
        24.25,
        24.25,
    ]

    # Check fit, transform and fit_transform
    ir = IsotonicRegression()
    ir.fit(x, y)
    assert_array_almost_equal(ir.transform(x), y_true, 4)
    assert_array_almost_equal(ir.fit_transform(x, y), y_true, 4)

代码总结：这段代码使用 R isotone 包 gpava(..., ties="secondary") 的权威结果作为 ground truth，验证 scikit-learn 对并列值的 secondary 处理方法一致性。期望值精确到小数点后 4 位，体现了算法数值实现的严格正确性要求。

19.12.1.3 快速预测优化验证

# 第 19 章 —— sklearn/tests/test_isotonic.py - test_fast_predict (第320-360行)
def test_fast_predict():
    # test that the faster prediction change doesn't
    # affect out-of-sample predictions:
    # https://github.com/scikit-learn/scikit-learn/pull/6206
    rng = np.random.RandomState(123)
    n_samples = 10**3
    # X values over the -10,10 range
    X_train = 20.0 * rng.rand(n_samples) - 10
    y_train = (
        np.less(rng.rand(n_samples), expit(X_train)).astype("int64").astype("float64")
    )

    weights = rng.rand(n_samples)
    # we also want to test that everything still works when some weights are 0
    weights[rng.rand(n_samples) < 0.1] = 0

    slow_model = IsotonicRegression(y_min=0, y_max=1, out_of_bounds="clip")
    fast_model = IsotonicRegression(y_min=0, y_max=1, out_of_bounds="clip")

    # Build interpolation function with ALL input data, not just the
    # non-redundant subset. The following 2 lines are taken from the
    # .fit() method, without removing unnecessary points
    X_train_fit, y_train_fit = slow_model._build_y(
        X_train, y_train, sample_weight=weights, trim_duplicates=False
    )
    slow_model._build_f(X_train_fit, y_train_fit)

    # fit with just the necessary data
    fast_model.fit(X_train, y_train, sample_weight=weights)

    X_test = 20.0 * rng.rand(n_samples) - 10
    y_pred_slow = slow_model.predict(X_test)
    y_pred_fast = fast_model.predict(X_test)

    assert_array_equal(y_pred_slow, y_pred_fast)

代码总结：这段代码验证了 trim_duplicates=True 优化不改变预测结果。慢模型保留所有 PAVA 输出点构建插值，快模型移除冗余共线点。1000 样本随机数据测试证明两者预测完全相等，确认了优化的正确性与性能收益。

---

19.13 设计中的取舍

为什么 LDA 不在 svd 求解器中支持 shrinkage？

SVD 求解器的核心优势在于不显式计算协方差矩阵，通过两阶段 SVD 直接从中心化数据得到判别方向，计算复杂度主要取决于样本数而非特征数，极其适合高维数据。引入 shrinkage 需要显式估计协方差矩阵（或其收缩版本），这会破坏 SVD 求解器"避免协方差矩阵"的设计初衷，且在 n_features > n_samples 时协方差矩阵本身就是奇异的，收缩估计需先降维或正则化，引入额外复杂度。用户需收缩估计，应使用 lsqr 或 eigen 求解器，它们天然基于协方差矩阵运算，可无缝集成 _cov 的收缩机制。

为什么 QDA 的 svd 求解器无法修复 n_samples <= n_features 的秩缺失？

SVD 分解 Xc = U S V^T 中，奇异值数量 min(n_samples, n_features)。当 n_samples <= n_features 时，最多只有 n_samples 个非零奇异值，协方差矩阵秩至多为 n_samples-1（中心化后），根本性秩缺失。reg_param 只是将现有奇异值平方加常数，无法凭空创造缺失的奇异值方向。而 Eigen 求解器配合 shrinkage='auto' 使用 Ledoit-Wolf 收缩估计，数学上等价于 Σ_shrunk = (1-λ)Σ_emp + λ*I，显式添加单位矩阵成分，可使协方差矩阵满秩。这是两种求解器数学基础（SVD vs 特征分解+收缩）决定的本质差异。

为什么 IsotonicRegression 存储阈值数组而非直接 pickle interp1d 对象？

scipy.interpolate.interp1d 对象包含复杂的内部状态（如分段多项式系数、缓存的搜索结构），跨 Python 版本、SciPy 版本、平台的 pickle 兼容性极差，且序列化体积大。存储 X_thresholds_、y_thresholds_ 两个一维数组（仅含非冗余断点）极其紧凑、稳定、跨版本兼容。反序列化时调用 _build_f 重建 interp1d 仅需微秒级开销。这是 scikit-learn 处理不可序列化科学计算对象的标准模式：存核心参数，重建计算图。

为什么 PAVA 算法使用 target 数组编码块结构而非显式链表？

target 数组利用数组索引隐式表示双向链表：块 [i..j] 满足 target[i]=j 且 target[j]=i。这种设计在 Cython 中具有极大优势：(1) 连续内存访问，缓存友好；(2) 无指针开销，intp 数组极其紧凑；(3) 回溯 i = target[i-1] 仅需一次数组访问；(4) 原地修改 y、w 存储聚合值，无需额外内存分配。显式链表在 Python/Cython 中需对象分配、指针追踪，开销大数量级。target 数组是算法逻辑与数据结构深度融合的典范。

---

19.14 动手练习

1. 对比LDA三大求解器的数值行为与适用边界

   • 阅读 sklearn/discriminant_analysis.py 中 _solve_svd、_solve_lstsq、_solve_eigen 实现（第142-300行）。

   • 回答问题：

     • 为什么 solver='lsqr' 不支持 transform() 降维？

     • solver='svd' 如何通过两阶段 SVD 避免显式计算协方差矩阵？

     • solver='eigen' 求解广义特征值问题 Sb v = λ Sw v 时，Sb 和 Sw 分别如何计算？

     • 在 n_features > n_samples 高维场景下，哪个求解器最稳健？为什么？

2. 分析QDA正则化机制对秩缺失协方差的修复能力

   • 阅读 sklearn/discriminant_analysis.py 中 QuadraticDiscriminantAnalysis.fit 与 _solve_svd、_solve_eigen（第500-610行）。

   • 回答问题：

     • solver='svd' 时 reg_param 如何作用于奇异值平方？为什么无法修复 n_samples <= n_features 导致的根本秩缺失？

     • solver='eigen' 时 shrinkage 参数（'auto' 或 float）如何通过 _cov -> ledoit_wolf 或 shrunk_covariance 修复协方差秩缺失？

     • fit() 中秩缺失检测逻辑 np.sum(scaling_class > self.tol) 的阈值 tol 默认值多少？如何区分 SVD 与 Eigen 模式下的报错信息？

3. 深入PAVA算法与IsotonicRegression完整流程

   • 阅读 sklearn/_isotonic.pyx 中 _inplace_contiguous_isotonic_regression（第15-60行）与 sklearn/isotonic.py 中 IsotonicRegression._build_y、_build_f（第250-300行）。

   • 回答问题：

     • PAVA 算法如何利用 target 数组表示块结构？回溯机制 i = target[i-1] 如何保证 O(n) 单次遍历？

     • _make_unique 合并重复 X 值时，如何用 np.finfo(dtype).resolution 作为浮点去重容差？float32 与 float64 行为有何差异？

     • trim_duplicates=True 移除哪些冗余点？为何能加速预测且不改变插值结果？

     • out_of_bounds='clip' 与 'nan' 在 _build_f 中如何通过 interp1d 的 bounds_error 和 np.clip 实现？

---

19.15 本章小结

本章小结表格汇总了判别分析与保序回归的核心概念、关键类与方法、算法原理及工程实现要点。

| 概念 | 解释 |

|------|------|

| LinearDiscriminantAnalysis | 线性判别分析，假设类别共享协方差矩阵，线性决策边界，支持降维 |

| QuadraticDiscriminantAnalysis | 二次判别分析，类别独立协方差矩阵，二次决策边界，需正则化防秩缺失 |

| _solve_svd (LDA) | SVD求解器，不显式计算协方差，适合高维数据，支持Array API，两阶段SVD |

| _solve_lstsq (LDA) | 最小二乘求解器，直接求解线性系统，支持任意协方差估计器，不支持降维 |

| _solve_eigen (LDA) | 特征值分解求解器，求解广义特征值问题，同时支持分类和降维，计算explained_variance_ratio |

| _solve_svd (QDA) | 逐类SVD，奇异值平方为方差，reg_param正则化防止奇异，要求n_samples > n_features |

| _solve_eigen (QDA) | 逐类特征分解，支持shrinkage/自定义协方差估计器，正则化可修复秩缺失 |

| DiscriminantAnalysisPredictionMixin | 统一预测接口，LDA线性决策函数，QDA二次决策函数（马氏距离+log先验） |

| IsotonicRegression | 保序回归，单调约束拟合，PAVA算法O(n)求解，线性插值预测，支持越界处理 |

| _inplace_contiguous_isotonic_regression | PAVA核心Cython实现，回溯合并块保证线性时间，原地修改数组 |

| _make_unique | 合并重复X值，加权平均y，按浮点精度容差去重，Cython/NumPy双实现 |

| check_increasing | Spearman秩相关系数+Fisher变换置信区间，自动判断单调方向，CI跨零报警 |

| trim_duplicates | 移除冗余插值点（y等于前后邻居），加速预测不改变结果 |

| getstate/setstate | 序列化支持，interp1d对象不可pickle，存储阈值数组重建插值函数 |

下一章中，我们将学习随机投影与 Dummy 估计器 —— 品味“化繁为简的实用主义”，探索 Johnson-Lindenstrauss 引理如何确保低维嵌入保持成对距离，以及基线估计器如何为复杂模型提供性能基准。

19.16 生活类比

想象判别分析与保序回归是统计建模中的“两大导航系统”： LDA = 共享地图的线性导航：所有车辆（类别）共用同一张路况图（共享协方差），规划直线路径（线性决策边界），三种路线规划算法（SVD/最小二乘/特征分解）适应不同地形（高维/需降维/任意协方差）。 QDA = 专属地图的曲线导航：每车辆（类别）持有私有地图（类别协方差），规划曲线路径（二次决策边界），但地图可能不全（秩缺失），需正则化（reg_param/shrinkage）填补或换算法（eigen+shrinkage）。 IsotonicRegression = 单调攀登的“台阶修路工”：PAVA算法像推土机，将乱序高程（y）推平成不降/不升台阶（单调拟合），合并相邻违规段（Pool Adjacent Violators），再铺设线性插值桥梁（interp1d）供新车通行，自动识别上坡/下坡（Spearman+Fisher）。 就像导航系统需权衡地图精度与计算速度，判别分析在共享/独立协方差、线性/二次边界、三大求解器间权衡；保序回归在单调约束与数据拟合间平衡，Cython PAVA保证O(n)极速铺路。

第 20 章 —— 随机投影与 Dummy 估计器 —— 品味“化繁为简的实用主义”

20.1 学习目标

• 难度：★★★☆☆（3/5）

• 预备知识：Python 基础、面向对象编程与 Markdown/代码阅读基础

• 理解 Johnson-Lindenstrauss 引理在随机投影中的理论指导作用

• 掌握高斯随机矩阵与稀疏随机矩阵的生成原理及统计特性

• 理解 BaseRandomProjection 基类如何统一管理拟合、变换、逆变换与标签系统

• 掌握 GaussianRandomProjection 与 SparseRandomProjection 的差异化实现与适用场景

• 理解 DummyClassifier 五大策略的预测逻辑、多输出/稀疏目标处理与样本权重感知机制

• 掌握 DummyRegressor 四大策略对应的统计量计算、加权分位数支持与多输出形状统一

• 了解随机投影与 Dummy 估计器在 scikit-learn 中的测试验证体系

20.2 生活类比

想象随机投影是一台高维数据的“压缩相机”。JL 引理好比光学定律，它告诉你至少需要多少像素（n_components）才能让照片里的物体距离不失真（eps）；高斯随机矩阵相当于标准镜头，每个像素独立采光，成像稠密、质量稳定，但存储开销大；稀疏随机矩阵则像稀疏镜头，只保留关键像素（密度 density），大部分置零，存储极小、成像极快，质量几乎不损。BaseRandomProjection 则是相机机身，统一了“拍照”（fit）、“投影”（transform）、“还原”（inverse_transform）接口，并能自动对焦（n_components_）。GaussianRandomProjection 是旗舰机，适合通用场景的稠密成像；SparseRandomProjection 是运动相机，适合大规模、内存敏感场景的稀疏成像，并支持密集输出模式（dense_output）。想象 Dummy 估计器是“懒惰的基线裁判”。DummyClassifier 好比只看比分不看球的裁判：most_frequent 策略永远判赢家赢，prior 策略判赢家赢但给出历史胜率作为置信度，stratified 策略按历史胜率随机判赢负，uniform 策略瞎蒙、各队概率均等，constant 策略则是死忠粉，永远判指定队赢。DummyRegressor 则像只背统计数字不看比赛的解说员：mean 策略总报平均分，median 策略总报中位数分，quantile 策略总报指定分位数分，constant 策略总报固定分。它们完全忽略比赛过程（X）、只看历史结果（y），同时支持样本权重（加权历史）、多输出（多项赛事）、稀疏目标（稀疏比分表）。

20.3 源码地图

sklearn/random_projection.py

├── johnson_lindenstrauss_min_dim() # JL引理计算最小投影维度

├── _check_density() # 稀疏矩阵密度校验与自动推断

├── _check_input_size() # 随机矩阵维度合法性校验

├── _gaussian_random_matrix() # 稠密高斯随机矩阵生成 N(0, 1/n_components)

├── _sparse_random_matrix() # 稀疏 Achlioptas/Li 随机矩阵生成 (CSR格式)

├── BaseRandomProjection (ABC)

│ ├── init() # 参数初始化: n_components, eps, compute_inverse_components, random_state

│ ├── _make_random_matrix() # 抽象方法: 子类实现具体矩阵生成

│ ├── _compute_inverse_components() # 计算伪逆用于逆变换

│ ├── fit() # 核心拟合: 自动推断n_components_, 生成components_, 可选预计算伪逆

│ ├── inverse_transform() # 逆变换: 预计算伪逆或即时计算, 输出始终密集

│ ├── sklearn_tags() # 标签系统: preserves_dtype, input_tags.sparse=True

├── GaussianRandomProjection

│ ├── init() # 继承基类参数

│ ├── _make_random_matrix() # 调用 _gaussian_random_matrix 生成稠密矩阵

│ ├── transform() # 稠密矩阵乘法: X @ components_.T

├── SparseRandomProjection

│ ├── init() # 扩展参数: density, dense_output

│ ├── _make_random_matrix() # 调用 sparse_random_matrix 生成CSR稀疏矩阵, 记录density

│ ├── transform() # 稀疏感知矩阵乘法: safe_sparse_dot, 支持dense_output控制

sklearn/dummy.py

├── DummyClassifier

│ ├── init() # 策略选择: most_frequent/prior/stratified/uniform/constant

│ ├── fit() # 核心拟合: class_distribution计算先验, 稀疏目标处理, 常数策略校验, 样本权重支持

│ ├── predict() # 五策略预测实现: 稀疏输出用_random_choice_csc, 密集输出用numpy索引/采样

│ ├── predict_proba() # 五策略概率输出: one-hot/先验常数/多项采样/均匀/常数one-hot

│ ├── predict_log_proba() # 概率取对数

│ ├── sklearn_tags() # 标签: poor_score=True, no_validation=True, input_tags.sparse=True

│ ├── score() # 允许X=None, 构造虚拟X复用父类逻辑

├── DummyRegressor

│ ├── init() # 策略选择: mean/median/quantile/constant

│ ├── fit() # 统计量计算: np.average/np.median/np.percentile/_weighted_percentile, 常数校验, 形状统一为(1, n_outputs)

│ ├── predict() # np.full广播生成预测, 支持return_std返回全零标准差

│ ├── sklearn_tags() # 标签: poor_score=True, no_validation=True, input_tags.sparse=True

│ ├── score() # 允许X=None, 构造虚拟X复用父类R^2逻辑

sklearn/tests/test_random_projection.py

├── test_invalid_jl_domain() # JL引理参数域校验

├── test_input_size_jl_min_dim() # JL引理广播语义测试

├── test_basic_property_of_random_matrix() # 随矩阵零均值单位范数验证

├── test_gaussian_random_matrix() # 高斯矩阵分布统计特性验证

├── test_sparse_random_matrix() # 稀疏矩阵三值集合与频率/方差验证

├── test_random_projection_embedding_quality()# JL引理契约验证: 距离比率落在[1-eps, 1+eps]

├── test_SparseRandomProj_output_representation() # 稀疏/密集输入输出格式控制验证

├── test_correct_RandomProjection_dimensions_embedding() # 自动/手动维度、密度、确定性复现验证

├── test_inverse_transform() # 逆变换往返一致性验证(预计算/即时计算、稀疏/密集)

├── test_random_projection_feature_names_out()# 特征名生成验证

├── test_random_projection_dtype_match() # dtype保持验证(float32/float64/int->float64)

├── test_random_projection_numerical_consistency() # float32/float64数值一致性验证

├── test_random_projection_transformer_invalid_input() # 无效输入拟合校验

├── test_try_to_transform_before_fit() # 未拟合调用transform异常验证

├── test_too_many_samples_to_find_a_safe_embedding() # 样本过多导致目标维度超限验证

├── test_warning_n_components_greater_than_n_features() # 组件数超过特征数警告验证

├── test_works_with_sparse_data() # 稀疏/密集输入数值等价性验证

├── test_johnson_lindenstrauss_min_dim() # JL引理小eps回归测试

sklearn/tests/test_dummy.py

├── test_most_frequent_and_prior_strategy() # 最频繁/先验策略预测与概率验证

├── test_stratified_strategy() # 分层采样策略预测分布验证

├── test_uniform_strategy() # 均匀采样策略预测分布验证

├── test_constant_strategy() # 常数策略预测与合法性校验

├── test_classifier_score_with_None() # X=None评分验证

├── test_classifier_prediction_independent_of_X() # 预测与X无关验证

├── test_mean_strategy_regressor() # 均值策略回归验证

├── test_median_strategy_regressor() # 中位数策略回归验证

├── test_quantile_strategy_regressor() # 分位数策略回归验证(含边界0/1)

├── test_constant_strategy_regressor() # 常数策略回归验证

├── test_dummy_regressor_sample_weight() # 样本权重加权统计量验证

├── test_dummy_regressor_return_std() # return_std全零验证

├── test_regressor_prediction_independent_of_X() # 回归预测与X无关验证

├── test_constant_strategy_sparse_target() # 分类器稀疏目标常数策略验证

├── test_stratified_strategy_sparse_target() # 分类器稀疏目标分层策略验证

├── test_uniform_strategy_sparse_target_warning() # 分类器稀疏目标均匀策略警告验证

├── test_most_frequent_and_prior_strategy_sparse_target() # 分类器稀疏目标最频繁/先验验证

├── test_quantile_strategy_multioutput_regressor() # 回归多分位数多输出验证

├── test_mean_strategy_multioutput_regressor() # 回归均值多输出验证

├── test_constant_strategy_multioutput_regressor() # 回归常数多输出验证

├── test_y_mean_attribute_regressor() # 回归mean属性验证

├── test_string_labels() # 字符串标签分类验证

├── test_median_strategy_multioutput_regressor() # 回归中位数多输出验证

├── test_constants_not_specified_regressor() # 回归常数未指定异常验证

├── test_quantile_invalid() # 回归分位数未指定异常验证

├── test_quantile_strategy_empty_train() # 回归分位数空训练异常验证

├── test_constant_strategy_multioutput() # 分类器常数多输出验证

├── test_constant_strategy_exceptions() # 分类器常数策略异常验证

├── test_classification_sample_weight() # 分类样本权重加权先验验证

├── test_dummy_regressor_on_3D_array() # 回归3D输入验证

├── test_dummy_classifier_on_3D_array() # 分类3D输入验证

├── test_regressor_score_with_None() # 回归X=None评分验证

├── test_feature_names_in_and_n_features_in_() # 特征名与特征数属性验证

├── test_dtype_of_classifier_probas() # 分类概率dtype验证

├── test_most_frequent_and_prior_strategy_with_2d_column_y() # 2D列向量y兼容性验证

├── test_most_frequent_and_prior_strategy_multioutput() # 分类最频繁/先验多输出验证

├── test_stratified_strategy_multioutput() # 分类分层策略多输出验证

├── test_uniform_strategy_multioutput() # 分类均匀策略多输出验证

20.4 Johnson-Lindenstrauss 引理 —— 随机投影的理论基石

JL 引理的核心结论在于：高维空间中的点集可以嵌入到低维空间，同时近似保持成对距离。目标维度 n_components 仅取决于样本数 n_samples 和失真率 eps，与原始特征数 n_features 无关。其公式为 n_components >= 4 * log(n_samples) / (eps^2 / 2 - eps^3 / 3)。

函数 johnson_lindenstrauss_min_dim 是该引理的工程实现。它首先对输入进行校验，要求 n_samples > 0 且 eps ∈ (0, 1)，同时支持标量与数组广播计算。分母计算采用 denominator = eps^2 / 2 - eps^3 / 3 的形式，这是为了在 eps 很小的时候避免直接计算 eps^2/2 可能导致的数值不稳定。最终返回向上取整的 int64 类型最小维度。在边界情况下，若 eps <= 0 或 >= 1，或 n_samples <= 0，均会抛出 ValueError；同时函数自动处理数组输入的广播语义，便于批量超参搜索。

源码路径：sklearn/random_projection.py - johnson_lindenstrauss_min_dim()（第 76-140 行）

@validate_params(
    {
        "n_samples": ["array-like", Interval(Real, 1, None, closed="left")],
        "eps": ["array-like", Interval(Real, 0, 1, closed="neither")],
    },
    prefer_skip_nested_validation=True,
)
def johnson_lindenstrauss_min_dim(n_samples, *, eps=0.1):
    # 将输入转换为 numpy 数组，支持标量和数组输入
    eps = np.asarray(eps)
    n_samples = np.asarray(n_samples)

    # 检查 eps 是否在 (0, 1) 范围内
    if np.any(eps <= 0.0) or np.any(eps >= 1):
        raise ValueError("The JL bound is defined for eps in ]0, 1[, got %r" % eps)

    # 检查 n_samples 是否大于 0
    if np.any(n_samples <= 0):
        raise ValueError(
            "The JL bound is defined for n_samples greater than zero, got %r"
            % n_samples
        )

    # 计算分母：eps^2 / 2 - eps^3 / 3
    # 这种形式避免了当 eps 很小的时候直接计算 eps^2/2 可能导致的数值不稳定
    denominator = (eps**2 / 2) - (eps**3 / 3)
    # 根据 JL 引理公式计算最小维度，并转换为 int64 类型
    return (4 * np.log(n_samples) / denominator).astype(np.int64)

点击查看 JL 引理维度计算流程图

flowchart TD A[输入 n_samples, eps] --> B{校验 eps ∈ (0,1) 和 n_samples > 0} B -- 不合法 --> C[抛出 ValueError] B -- 合法 --> D[转为 numpy 数组 支持广播] D --> E[计算分母 denominator = eps²/2 - eps³/3] E --> F[计算 n_components = ceil(4 * log(n_samples) / denominator)] F --> G[返回 int64 类型结果]

上述代码实现了 JL 引理的核心公式计算。首先通过 np.asarray 将输入标准化为数组，从而天然支持标量与数组的广播运算。校验逻辑使用 np.any 对数组进行逐元素判断，确保所有 eps 严格位于 (0, 1) 且所有 n_samples 为正。分母采用 eps**2 / 2 - eps**3 / 3 的展开式而非单纯的 eps**2 / 2，这是因为当 eps 趋近于 0 时，高阶项 eps**3 / 3 虽小但能提供更精确的泰勒近似，避免分母过小导致的数值溢出或精度损失。最后通过 astype(np.int64) 向上取整（因除法结果为浮点，转整数即隐含向上取整语义），返回满足 JL 契约的最小投影维度。

20.5 高斯随机矩阵 —— 稠密投影的“标准正态建造师”

高斯随机矩阵的数学分布与归一化策略是：矩阵元素 i.i.d. 采样自 N(0, 1 / n_components)，使得列向量期望模长为 1（零均值、单位范数），从而满足 JL 引理要求。实现上使用 rng.normal(loc=0.0, scale=1.0 / np.sqrt(n_components), ...) 直接生成。

函数 _gaussian_random_matrix 实现细节包括：参数校验由 _check_input_size 保证 n_components, n_features > 0；随机状态管理通过 check_random_state 统一处理 seed/RandomState/None；返回密集 ndarray，形状为 (n_components, n_features)。从统计特性验证来看（测试视角），经验均值约等于 0，经验方差约等于 1 / n_components，列向量范数集中在 1 附近（单位球面投影）。

源码路径：sklearn/random_projection.py - _gaussian_random_matrix()（第 187-215 行）

def _gaussian_random_matrix(n_components, n_features, random_state=None):
    # 检查输入维度的合法性
    _check_input_size(n_components, n_features)
    # 获取随机数生成器
    rng = check_random_state(random_state)
    # 生成服从 N(0, 1/n_components) 分布的矩阵
    # scale 参数设置为 1/sqrt(n_components) 确保方差为 1/n_components
    components = rng.normal(
        loc=0.0, scale=1.0 / np.sqrt(n_components), size=(n_components, n_features)
    )
    return components

点击查看高斯随机矩阵生成流程图

flowchart TD A[输入 n_components, n_features, random_state] --> B[校验维度 > 0] B --> C[获取随机数生成器 rng] C --> D[从 N(0, 1/√n_components) 采样生成矩阵] D --> E[返回稠密 ndarray (n_components, n_features)]

该函数仅三行核心逻辑：校验维度、获取随机数生成器、按指定分布采样。scale=1.0 / np.sqrt(n_components) 是关键：因为正态分布的方差为 scale²，设为 1/n_components 使得每列由 n_components 个独立同分布元素组成，其期望平方和（即模长平方）为 n_components * (1/n_components) = 1，从而满足 JL 引理对正交投影矩阵列向量单位范数的要求。返回的稠密 ndarray 直接用于后续的矩阵乘法投影。

20.6 稀疏随机矩阵 —— Achlioptas 与 Li 等人的“稀疏变奏曲”

稀疏随机矩阵采用三点分布的参数化设计：设 s = 1 / density，非零元素为 ±sqrt(s) / sqrt(n_components)，零元素概率为 1 - 1/s。当 density='auto' 时取 1 / sqrt(n_features)（Li et al. 推荐最小密度）；density=1/3 时退化为 Achlioptas 原始矩阵（±1/sqrt(n_components) 各 1/6，0 为 2/3）。

CSR 稀疏矩阵的高效构建流程为：分行采样，每行非零个数服从二项分布 Binomial(n_features, density)；列索引无放回采样，使用 sample_without_replacement 保证行内列不重复；符号随机化，非零位置 50%/50% 分配 ±1；组装 CSR，indptr 记录行偏移，indices 存列号，data 存缩放后的值。当 density == 1.0 时，走密集退化路径，直接生成密集矩阵再转 CSR，避免索引开销。统计特性验证方面（测试视角），需校验三值集合 {0, +sqrt(s)/sqrt(k), -sqrt(s)/sqrt(k)}，频率校验 P(0)=1-1/s, P(±)=1/(2s)，方差校验符合二项分布。

源码路径：sklearn/random_projection.py - _sparse_random_matrix()（第 217-285 行）

def _sparse_random_matrix(n_components, n_features, density="auto", random_state=None):
    # 检查输入维度的合法性
    _check_input_size(n_components, n_features)
    # 处理 density 参数，支持'auto'和显式浮点数
    density = _check_density(density, n_features)
    rng = check_random_state(random_state)

    if density == 1:
        # 当密度为1时，直接生成密集矩阵（所有元素非零）
        # 使用二项分布生成符号，然后缩放
        components = rng.binomial(1, 0.5, (n_components, n_features)) * 2 - 1
        return 1 / np.sqrt(n_components) * components

    else:
        # 生成非零元素的位置
        indices = []
        offset = 0
        indptr = [offset]
        for _ in range(n_components):
            # 当前行的非零元素数量服从二项分布
            n_nonzero_i = rng.binomial(n_features, density)
            # 无放回采样列索引
            indices_i = sample_without_replacement(
                n_features, n_nonzero_i, random_state=rng
            )
            indices.append(indices_i)
            offset += n_nonzero_i
            indptr.append(offset)

        # 展平所有行的列索引
        indices = np.concatenate(indices)

        # 非零元素的符号：50% 正，50% 负
        data = rng.binomial(1, 0.5, size=np.size(indices)) * 2 - 1

        # 构建 CSR 稀疏矩阵
        components = sp.csr_matrix(
            (data, indices, indptr), shape=(n_components, n_features)
        )

        # 按照理论缩放，使得列向量期望模长为1
        return np.sqrt(1 / density) / np.sqrt(n_components) * components

点击查看稀疏随机矩阵生成流程图

flowchart TD A[输入 n_components, n_features, density, random_state] --> B[校验维度 > 0] B --> C[解析 density: auto -> 1/√n_features] C --> D{density == 1.0?} D -- 是 --> E[生成密集符号矩阵 ±1 并缩放 1/√n_components] D -- 否 --> F[逐行采样非零个数 ~ Binomial(n_features, density)] F --> G[无放回采样列索引 sample_without_replacement] G --> H[符号随机化 ±1 各 50%] H --> I[组装 CSR: data, indices, indptr] I --> J[缩放 √(1/density) / √n_components] E --> K[返回 CSR 稀疏矩阵] J --> K

代码分为密度为 1 的快速路径与稀疏路径。稀疏路径中，rng.binomial(n_features, density) 为每行生成非零元素个数，符合二项分布期望；sample_without_replacement 保证同一行内列索引不重复，这是 CSR 格式合法性的前提。data 数组存储 ±1 符号，后乘以缩放因子 np.sqrt(1 / density) / np.sqrt(n_components)。该因子推导自：非零元素值为 ±√s，方差为 s；每列期望非零数为 density * n_components = n_components / s，列方差和为 (n_components / s) * s = n_components，除以 n_components 归一化后列模长期望为 1。CSR 三元组 data, indices, indptr 分别存非零值、列索引、行起始偏移，是 scipy.sparse 标准格式。

20.7 基类 BaseRandomProjection —— 统一拟合、逆变换与标签系统

BaseRandomProjection 基类实现了参数约束与自动维度推断：当 n_components='auto' 时调用 johnson_lindenstrauss_min_dim 计算 n_components_；校验 n_components_ 不超过 n_features，否则报错；手动指定过大仅警告。compute_inverse_components 控制是否在 fit 阶段预计算伪逆。

核心拟合流程 fit 使用 validate_data 接受 csr/csc 稀疏矩阵，仅用于获取 shape 与 dtype；生成随机矩阵时调用子类 _make_random_matrix，并转 dtype 避免拷贝；可选预计算伪逆 linalg.pinv(components.toarray())，密集存储；设置 _n_features_out 供 get_feature_names_out 使用。

逆变换 inverse_transform 的实现策略为：预计算伪逆时直接矩阵乘法 X @ inverse_components_.T；延迟计算时每次调用临时 pinv，节省内存但重复计算；输入校验使用 check_array 接受稀疏/稠密，输出始终为密集数组。

标签系统与输出特性声明：preserves_dtype=['float64','float32'] 保持精度；input_tags.sparse=True 声明支持稀疏输入；继承 ClassNamePrefixFeaturesOutMixin 自动生成特征名如 gaussianrandomprojection0。

源码路径：sklearn/random_projection.py - BaseRandomProjection.__init__()（第 287-305 行）

@abstractmethod
def __init__(
    self,
    n_components="auto",
    *,
    eps=0.1,
    compute_inverse_components=False,
    random_state=None,
):
    # 初始化核心参数
    self.n_components = n_components
    self.eps = eps
    self.compute_inverse_components = compute_inverse_components
    self.random_state = random_state

源码路径：sklearn/random_projection.py - BaseRandomProjection.fit()（第 307-345 行）

@_fit_context(prefer_skip_nested_validation=True)
def fit(self, X, y=None):
    # 验证输入数据，仅接受CSR和CSC格式的稀疏矩阵，限制dtype为float32/float64
    X = validate_data(
        self, X, accept_sparse=["csr", "csc"], dtype=[np.float64, np.float32]
    )

    n_samples, n_features = X.shape

    if self.n_components == "auto":
        # 自动计算目标维度
        self.n_components_ = johnson_lindenstrauss_min_dim(
            n_samples=n_samples, eps=self.eps
        )

        if self.n_components_ <= 0:
            raise ValueError(
                "eps=%f and n_samples=%d lead to a target dimension of "
                "%d which is invalid" % (self.eps, n_samples, self.n_components_)
            )

        elif self.n_components_ > n_features:
            raise ValueError(
                "eps=%f and n_samples=%d lead to a target dimension of "
                "%d which is larger than the original space with "
                "n_features=%d"
                % (self.eps, n_samples, self.n_components_, n_features)
            )
    else:
        if self.n_components > n_features:
            warnings.warn(
                "The number of components is higher than the number of"
                " features: n_features < n_components (%s < %s)."
                "The dimensionality of the problem will not be reduced."
                % (n_features, self.n_components),
                DataDimensionalityWarning,
            )

        self.n_components_ = self.n_components

    # 生成投影矩阵并转换为输入数据的dtype，避免拷贝
    self.components_ = self._make_random_matrix(
        self.n_components_, n_features
    ).astype(X.dtype, copy=False)

    if self.compute_inverse_components:
        self.inverse_components_ = self._compute_inverse_components()

    # 设置特征输出数量，供 get_feature_names_out 使用
    self._n_features_out = self.n_components

    return self

源码路径：sklearn/random_projection.py - BaseRandomProjection.inverse_transform()（第 347-370 行）

def inverse_transform(self, X):
    # 检查是否已拟合
    check_is_fitted(self)

    # 验证输入数据，接受稀疏/稠密，限制dtype
    X = check_array(X, dtype=[np.float64, np.float32], accept_sparse=("csr", "csc"))

    if self.compute_inverse_components:
        # 使用预计算的伪逆
        return X @ self.inverse_components_.T

    # 即时计算伪逆
    inverse_components = self._compute_inverse_components()
    return X @ inverse_components.T

源码路径：sklearn/random_projection.py - BaseRandomProjection.__sklearn_tags__()（第 372-377 行）

def __sklearn_tags__(self):
    tags = super().__sklearn_tags__()
    tags.transformer_tags.preserves_dtype = ["float64", "float32"]
    tags.input_tags.sparse = True
    return tags

点击查看 BaseRandomProjection 核心流程图

flowchart TD A[fit(X, y)] --> B[validate_data: 仅取shape/dtype, 支持csr/csc] B --> C{n_components == 'auto'?} C -- 是 --> D[johnson_lindenstrauss_min_dim 计算 n_components_] C -- 否 --> E[检查 n_components > n_features 发警告] D --> F[校验 n_components_ <= n_features] E --> G[调用子类 _make_random_matrix 生成 components_] F --> G G --> H[astype(X.dtype, copy=False) 保持精度] H --> I{compute_inverse_components?} I -- 是 --> J[预计算 inverse_components_ = pinv(components_)] I -- 否 --> K[延迟到 inverse_transform 时计算] J --> L[设置 _n_features_out] K --> L L --> M[返回 self] N[inverse_transform(X)] --> O[check_is_fitted] O --> P[check_array 接受稀疏/稠密] P --> Q{预计算了 inverse_components_?} Q -- 是 --> R[X @ inverse_components_.T] Q -- 否 --> S[_compute_inverse_components 即时计算 pinv] S --> R R --> T[返回密集数组]

fit 方法中 validate_data 设置 accept_sparse=["csr", "csc"] 与 dtype=[np.float64, np.float32]，是因为随机投影仅需输入的形状与数值类型，不修改数据本身；限制 dtype 可避免后续矩阵乘法产生不必要的类型提升。components_.astype(X.dtype, copy=False) 实现零拷贝 dtype 对齐。inverse_transform 输出始终为密集数组，因为伪逆 pinv 结果本身是密集的，且重构原始空间通常需要完整表示；预计算伪逆以空间换时间（存 inverse_components_），即时计算以时间换空间（每次 pinv），用户可通过 compute_inverse_components 参数权衡。

20.8 GaussianRandomProjection —— 稠密投影的“开箱即用”封装

GaussianRandomProjection 通过继承与最小实现完成封装：仅需实现 _make_random_matrix 调用 _gaussian_random_matrix；transform 直接稠密矩阵乘法 X @ components_.T；继承基类的 fit、inverse_transform、参数校验、标签系统。

数据类型一致性保障体现在：fit 时 components_.astype(X.dtype, copy=False) 保持 float32/float64；transform 经 validate_data 统一 dtype，输出 dtype 与输入一致。嵌入质量保证方面（测试视角），JL 引理契约要求投影后成对距离比率落在 [1-eps, 1+eps] 区间；确定性复现要求相同 random_state 产生相同 components_ 与投影结果。

源码路径：sklearn/random_projection.py - GaussianRandomProjection.__init__()（第 315-325 行）

def __init__(
    self,
    n_components="auto",
    *,
    eps=0.1,
    compute_inverse_components=False,
    random_state=None,
):
    # 直接透传参数给基类
    super().__init__(
        n_components=n_components,
        eps=eps,
        compute_inverse_components=compute_inverse_components,
        random_state=random_state,
    )

源码路径：sklearn/random_projection.py - GaussianRandomProjection._make_random_matrix()（第 330-345 行）

def _make_random_matrix(self, n_components, n_features):
    # 生成高斯随机矩阵
    random_state = check_random_state(self.random_state)
    return _gaussian_random_matrix(
        n_components, n_features, random_state=random_state
    )

源码路径：sklearn/random_projection.py - GaussianRandomProjection.transform()（第 347-360 行）

def transform(self, X):
    # 检查是否已拟合
    check_is_fitted(self)
    # 验证输入数据，接受稀疏/稠密，重置校验状态，限制dtype
    X = validate_data(
        self,
        X,
        accept_sparse=["csr", "csc"],
        reset=False,
        dtype=[np.float64, np.float32],
    )

    # 稠密矩阵乘法实现投影
    return X @ self.components_.T

点击查看 GaussianRandomProjection 变换流程图

flowchart TD A[transform(X)] --> B[check_is_fitted] B --> C[validate_data: 接受csr/csc, reset=False, dtype限制] C --> D[稠密矩阵乘法: X @ components_.T] D --> E[返回稠密投影结果]

transform 方法极其简洁：check_is_fitted 确保已拟合；validate_data 以 reset=False 复用 fit 时的校验器，仅检查样本数与特征数一致性；稠密矩阵乘法 @ 自动处理稀疏输入（scipy.sparse 矩阵支持 @ 运算符），结果始终为稠密 ndarray，符合 JL 引理嵌入的稠密输出预期。

20.9 SparseRandomProjection —— 稀疏投影的“内存与速度双赢”封装

SparseRandomProjection 扩展了参数 density 与 dense_output：density 控制稀疏度，支持 'auto' 与显式浮点数；dense_output 控制稀疏输入时输出是否强制密集（小 n_components 时更快）。

稀疏感知的矩阵乘法使用 safe_sparse_dot(X, components_.T, dense_output=...)，自动处理稠密@稀疏、稀疏@稀疏、稀疏@稠密，按 dense_output 决定输出格式。密度属性持久化体现在 fit 中 self.density_ = _check_density(...) 记录实际使用密度，便于审计与复现（如测试断言 density_ ≈ 0.03）。稀疏目标与逆变换兼容性上，继承基类 inverse_transform，稀疏 components_ 会在 pinv 时自动转密集；测试覆盖稀疏/稠密输入、dense_output=True/False、逆变换往返一致性。

源码路径：sklearn/random_projection.py - SparseRandomProjection.__init__()（第 435-452 行）

def __init__(
    self,
    n_components="auto",
    *,
    density="auto",
    eps=0.1,
    dense_output=False,
    compute_inverse_components=False,
    random_state=None,
):
    super().__init__(
        n_components=n_components,
        eps=eps,
        compute_inverse_components=compute_inverse_components,
        random_state=random_state,
    )

    self.dense_output = dense_output
    self.density = density

源码路径：sklearn/random_projection.py - SparseRandomProjection._make_random_matrix()（第 463-478 行）

def _make_random_matrix(self, n_components, n_features):
    # 生成稀疏随机矩阵并记录实际使用的密度
    random_state = check_random_state(self.random_state)
    self.density_ = _check_density(self.density, n_features)
    return _sparse_random_matrix(
        n_components, n_features, density=self.density_, random_state=random_state
    )

源码路径：sklearn/random_projection.py - SparseRandomProjection.transform()（第 480-495 行）

def transform(self, X):
    # 检查是否已拟合
    check_is_fitted(self)
    # 验证输入数据
    X = validate_data(
        self,
        X,
        accept_sparse=["csr", "csc"],
        reset=False,
        dtype=[np.float64, np.float32],
    )

    # 使用安全的稀疏矩阵乘法，根据dense_output参数控制输出格式
    return safe_sparse_dot(X, self.components_.T, dense_output=self.dense_output)

点击查看 SparseRandomProjection 变换流程图

flowchart TD A[transform(X)] --> B[check_is_fitted] B --> C[validate_data: 接受csr/csc, reset=False, dtype限制] C --> D[safe_sparse_dot(X, components_.T, dense_output=self.dense_output)] D --> E{dense_output=True 或 输入稠密?} E -- 是 --> F[返回稠密 ndarray] E -- 否 --> G[返回稀疏 CSR 矩阵]

safe_sparse_dot 是关键：它会根据输入矩阵的稀疏性与 dense_output 参数智能选择乘法内核。当输入稀疏且 dense_output=False 时，输出保持 CSR 稀疏，极大节省内存；当 dense_output=True 或输入稠密时，输出稠密 ndarray，避免稀疏矩阵乘法的索引开销。self.density_ 在 _make_random_matrix 中记录，便于后续检查或调试时确认实际稀疏度。

20.10 DummyClassifier —— 分类基线的“策略工厂”

DummyClassifier 实现了五大策略的预测逻辑对比：most_frequent 预测众类，predict_proba 为 one-hot；prior 预测众类，predict_proba 为经验先验分布（常数行）；stratified 按先验分布多项采样，predict 取采样类，predict_proba 为 one-hot 采样；uniform 类别均匀采样，predict_proba 为均匀分布；constant 用户指定常数标签，predict_proba 为对应 one-hot。

多输出与稀疏目标的统一处理体现在：fit 中 class_distribution 逐列计算 classes_, n_classes_, class_prior_；稀疏目标 y 时 sparse_output_=True，预测输出 CSC 稀疏矩阵；uniform 策略不支持稀疏目标（会强制转密集并警告）。常数策略的合法性校验在 fit 时检查 constant 是否出现在每列的 classes_ 中，形状校验要求 constant 为 (n_outputs,) 兼容形状。样本权重感知的先验估计通过 class_distribution 内部使用 sample_weight 加权计算频率实现，stratified 采样、prior 概率均反映加权分布。与 X 无关的预测特性表现为 predict/predict_proba 完全忽略 X 的值，仅使用 n_samples；score 允许 X=None，内部构造虚拟 X 复用父类逻辑。标签系统声明 poor_score=True 标记为基线模型，no_validation=True 跳过输入校验。

源码路径：sklearn/dummy.py - DummyClassifier.__init__()（第 114-120 行）

def __init__(self, *, strategy="prior", random_state=None, constant=None):
    self.strategy = strategy
    self.random_state = random_state
    self.constant = constant

源码路径：sklearn/dummy.py - DummyClassifier.fit()（第 138-210 行）

@_fit_context(prefer_skip_nested_validation=True)
def fit(self, X, y, sample_weight=None):
    # 验证输入数据，跳过check_array，仅检查一致性
    validate_data(self, X, skip_check_array=True)

    self._strategy = self.strategy

    # 处理稀疏目标的特殊情况
    if self._strategy == "uniform" and sp.issparse(y):
        y = y.toarray()
        warnings.warn(
            (
                "A local copy of the target data has been converted "
                "to a numpy array. Predicting on sparse target data "
                "with the uniform strategy would not save memory "
                "and would be slower."
            ),
            UserWarning,
        )

    self.sparse_output_ = sp.issparse(y)

    if not self.sparse_output_:
        y = np.asarray(y)
        y = np.atleast_1d(y)

    if y.ndim == 1:
        y = np.reshape(y, (-1, 1))

    self.n_outputs_ = y.shape[1]

    check_consistent_length(X, y)

    if sample_weight is not None:
        sample_weight = _check_sample_weight(sample_weight, X)

    if self._strategy == "constant":
        if self.constant is None:
            raise ValueError(
                "Constant target value has to be specified "
                "when the constant strategy is used."
            )
        else:
            constant = np.reshape(np.atleast_1d(self.constant), (-1, 1))
            if constant.shape[0] != self.n_outputs_:
                raise ValueError(
                    "Constant target value should have shape (%d, 1)."
                    % self.n_outputs_
                )

    # 计算类别分布（支持样本权重）
    (self.classes_, self.n_classes_, self.class_prior_) = class_distribution(
        y, sample_weight
    )

    if self._strategy == "constant":
        for k in range(self.n_outputs_):
            if not any(constant[k][0] == c for c in self.classes_[k]):
                # 检查常数值是否存在于训练数据中
                err_msg = (
                    "The constant target value must be present in "
                    "the training data. You provided constant={}. "
                    "Possible values are: {}.".format(
                        self.constant, self.classes_[k].tolist()
                    )
                )
                raise ValueError(err_msg)

    if self.n_outputs_ == 1:
        self.n_classes_ = self.n_classes_[0]
        self.classes_ = self.classes_[0]
        self.class_prior_ = self.class_prior_[0]

    return self

源码路径：sklearn/dummy.py - DummyClassifier.predict()（第 212-270 行）

def predict(self, X):
    """Perform classification on test vectors X."""
    check_is_fitted(self)

    # 获取样本数量
    n_samples = _num_samples(X)
    rs = check_random_state(self.random_state)

    # 处理多输出情况下的统一表示
    n_classes_ = self.n_classes_
    classes_ = self.classes_
    class_prior_ = self.class_prior_
    constant = self.constant
    if self.n_outputs_ == 1:
        # 确保即使是单输出也使用列表形式便于统一处理
        n_classes_ = [n_classes_]
        classes_ = [classes_]
        class_prior_ = [class_prior_]
        constant = [constant]
    # 预先计算概率（仅在需要时）
    if self._strategy == "stratified":
        proba = self.predict_proba(X)
        if self.n_outputs_ == 1:
            proba = [proba]

    if self.sparse_output_:
        class_prob = None
        if self._strategy in ("most_frequent", "prior"):
            classes_ = [np.array([cp.argmax()]) for cp in class_prior_]

        elif self._strategy == "stratified":
            class_prob = class_prior_

        elif self._strategy == "uniform":
            raise ValueError(
                "Sparse target prediction is not "
                "supported with the uniform strategy"
            )

        elif self._strategy == "constant":
            classes_ = [np.array([c]) for c in constant]

        y = _random_choice_csc(n_samples, classes_, class_prob, self.random_state)
    else:
        if self._strategy in ("most_frequent", "prior"):
            y = np.tile(
                [
                    classes_[k][class_prior_[k].argmax()]
                    for k in range(self.n_outputs_)
                ],
                [n_samples, 1],
            )

        elif self._strategy == "stratified":
            y = np.vstack(
                [
                    classes_[k][proba[k].argmax(axis=1)]
                    for k in range(self.n_outputs_)
                ]
            ).T

        elif self._strategy == "uniform":
            ret = [
                classes_[k][rs.randint(n_classes_[k], size=n_samples)]
                for k in range(self.n_outputs_)
            ]
            y = np.vstack(ret).T

        elif self._strategy == "constant":
            y = np.tile(self.constant, (n_samples, 1))

            if self.n_outputs_ == 1:
                y = np.ravel(y)

        return y

源码路径：sklearn/dummy.py - DummyClassifier.predict_proba()（第 272-320 行）

def predict_proba(self, X):
    """
    Return probability estimates for the test vectors X.
    """
    check_is_fitted(self)

    # 获取样本数量
    n_samples = _num_samples(X)
    rs = check_random_state(self.random_state)

    # 处理多输出情况下的统一表示
    n_classes_ = self.n_classes_
    classes_ = self.classes_
    class_prior_ = self.class_prior_
    constant = self.constant
    if self.n_outputs_ == 1:
        # 确保即使是单输出也使用列表形式便于统一处理
        n_classes_ = [n_classes_]
        classes_ = [classes_]
        class_prior_ = [class_prior_]
        constant = [constant]

    P = []
    for k in range(self.n_outputs_):
        if self._strategy == "most_frequent":
            ind = class_prior_[k].argmax()
            out = np.zeros((n_samples, n_classes_[k]), dtype=np.float64)
            out[:, ind] = 1.0
        elif self._strategy == "prior":
            out = np.ones((n_samples, 1)) * class_prior_[k]

        elif self._strategy == "stratified":
            out = rs.multinomial(1, class_prior_[k], size=n_samples)
            out = out.astype(np.float64)

        elif self._strategy == "uniform":
            out = np.ones((n_samples, n_classes_[k]), dtype=np.float64)
            out /= n_classes_[k]

        elif self._strategy == "constant":
            ind = np.where(classes_[k] == constant[k])
            out = np.zeros((n_samples, n_classes_[k]), dtype=np.float64)
            out[:, ind] = 1.0

        P.append(out)

    if self.n_outputs_ == 1:
        P = P[0]

    return P

源码路径：sklearn/dummy.py - DummyClassifier.__sklearn_tags__()（第 338-343 行）

def __sklearn_tags__(self):
    tags = super().__sklearn_tags__()
    tags.input_tags.sparse = True
    tags.classifier_tags.poor_score = True
    tags.no_validation = True
    return tags

源码路径：sklearn/dummy.py - DummyClassifier.score()（第 345-360 行）

def score(self, X, y, sample_weight=None):
    """Return the mean accuracy on the given test data and labels."""
    if X is None:
        X = np.zeros(shape=(len(y), 1))
    return super().score(X, y, sample_weight)

点击查看 DummyClassifier 核心流程图

flowchart TD A[fit(X, y, sample_weight)] --> B[validate_data skip_check_array] B --> C{uniform + 稀疏y?} C -- 是 --> D[转密集并警告] C -- 否 --> E[记录 sparse_output_] E --> F[规范化 y 为 (n_samples, n_outputs)] F --> G[检查长度一致性] G --> H{constant 策略?} H -- 是 --> I[校验 constant 形状与合法性] H -- 否 --> J[class_distribution 计算 classes_, n_classes_, class_prior_] I --> J J --> K[单输出时展平属性] K --> L[返回 self] M[predict(X)] --> N[check_is_fitted] N --> O[获取 n_samples, random_state] O --> P[统一多输出为列表形式] P --> Q{stratified?} Q -- 是 --> R[调用 predict_proba 获取 proba] Q -- 否 --> S{稀疏输出?} R --> S S -- 是 --> T[_random_choice_csc 生成 CSC 稀疏预测] S -- 否 --> U{策略分发} U -- most_frequent/prior --> V[np.tile 众类/先验最大类] U -- stratified --> W[proba.argmax 转类别] U -- uniform --> X[rs.randint 均匀采样] U -- constant --> Y[np.tile 常数标签] V --> Z[单输出时 ravel] W --> Z X --> Z Y --> Z Z --> AA[返回预测] AB[predict_proba(X)] --> AC[check_is_fitted] AC --> AD[获取 n_samples, random_state] AD --> AE[统一多输出为列表] AE --> AF[逐输出构造概率矩阵] AF --> AG{策略分发} AG -- most_frequent --> AH[one-hot 众类] AG -- prior --> AI[广播 class_prior_] AG -- stratified --> AJ[rs.multinomial 多项采样] AG -- uniform --> AK[均匀分布 1/n_classes] AG -- constant --> AL[one-hot 常数类] AH --> AM[合并列表/单输出] AI --> AM AJ --> AM AK --> AM AL --> AM AM --> AN[返回概率]

fit 方法中 validate_data(..., skip_check_array=True) 跳过对 X 的详细校验（因 no_validation=True），仅检查长度一致性。class_distribution 返回三元组：classes_ 为每列唯一类别数组列表，n_classes_ 为每列类别数列表，class_prior_ 为每列归一化频率（含样本权重）列表。稀疏目标时 sparse_output_=True，预测走 _random_choice_csc 生成 CSC 矩阵；uniform 策略强制转密集并警告，因为均匀采样无法利用稀疏结构加速。predict 中 stratified 策略先调用 predict_proba 获取采样的 one-hot 矩阵，再 argmax 还原类别，保证 predict 与 predict_proba 采样一致。predict_proba 逐输出构造概率矩阵，最终合并为列表或单数组。

20.11 DummyRegressor —— 回归基线的“统计量速查表”

DummyRegressor 实现了四大策略的统计量对应关系：mean 对应 np.average(y, axis=0, weights=sample_weight)；median 无权重用 np.median，有权重用 _weighted_percentile(y, w, 50)；quantile 无权重用 np.percentile(q=quantile*100)，有权重用 _weighted_percentile；constant 用户提供常数，支持标量、数组、稀疏矩阵。

多输出与形状统一机制为：y 统一重塑为 (n_samples, n_outputs)，constant_ 存为 (1, n_outputs)；predict 通过 np.full 广播生成 (n_samples, n_outputs) 预测；return_std=True 返回全零标准差（符合回归器接口约定）。样本权重的加权分位数支持体现在 median/quantile 策略在有权重时调用 _weighted_percentile，该函数实现加权分位数计算（见 sklearn/utils/stats.py）。

异常与边界处理包括：quantile 策略必须指定 quantile ∈ [0,1]，否则 fit 报错；constant 策略必须提供 constant，否则 TypeError；空 y 在 check_array 阶段拦截，quantile 空训练抛 IndexError。与 X 无关的预测特性表现为 predict 忽略 X 内容，仅用 _num_samples(X) 获取行数；score 允许 X=None，构造虚拟 X 复用父类 RegressorMixin.score。标签系统声明 poor_score=True、no_validation=True、input_tags.sparse=True。

源码路径：sklearn/dummy.py - DummyRegressor.__init__()（第 395-400 行）

def __init__(self, *, strategy="mean", constant=None, quantile=None):
    self.strategy = strategy
    self.constant = constant
    self.quantile = quantile

源码路径：sklearn/dummy.py - DummyRegressor.fit()（第 410-480 行）

@_fit_context(prefer_skip_nested_validation=True)
def fit(self, X, y, sample_weight=None):
    """Fit the baseline regressor."""
    validate_data(self, X, skip_check_array=True)

    y = check_array(y, ensure_2d=False, input_name="y")
    if len(y) == 0:
        raise ValueError("y must not be empty.")

    if y.ndim == 1:
        y = np.reshape(y, (-1, 1))
    self.n_outputs_ = y.shape[1]

    check_consistent_length(X, y, sample_weight)

    if sample_weight is not None:
        sample_weight = _check_sample_weight(sample_weight, X)

    if self.strategy == "mean":
        self.constant_ = np.average(y, axis=0, weights=sample_weight)

    elif self.strategy == "median":
        if sample_weight is None:
            self.constant_ = np.median(y, axis=0)
        else:
            self.constant_ = _weighted_percentile(
                y, sample_weight, percentile_rank=50.0
            )

    elif self.strategy == "quantile":
        if self.quantile is None:
            raise ValueError(
                "When using `strategy='quantile', you have to specify the desired "
                "quantile in the range [0, 1]."
            )
        percentile_rank = self.quantile * 100.0
        if sample_weight is None:
            self.constant_ = np.percentile(y, axis=0, q=percentile_rank)
        else:
            self.constant_ = _weighted_percentile(
                y, sample_weight, percentile_rank=percentile_rank
            )

    elif self.strategy == "constant":
        if self.constant is None:
            raise TypeError(
                "Constant target value has to be specified "
                "when the constant strategy is used."
            )

        self.constant_ = check_array(
            self.constant,
            accept_sparse=["csr", "csc", "coo"],
            ensure_2d=False,
            ensure_min_samples=0,
        )

        if self.n_outputs_ != 1 and self.constant_.shape[0] != y.shape[1]:
            raise ValueError(
                "Constant target value should have shape (%d, 1)." % y.shape[1]
            )

    self.constant_ = np.reshape(self.constant_, (1, -1))
    return self

源码路径：sklearn/dummy.py - DummyRegressor.predict()（第 482-505 行）

def predict(self, X, return_std=False):
    """Perform classification on test vectors X."""
    check_is_fitted(self)
    n_samples = _num_samples(X)

    y = np.full(
        (n_samples, self.n_outputs_),
        self.constant_,
        dtype=np.array(self.constant_).dtype,
    )
    y_std = np.zeros((n_samples, self.n_outputs_))

    if self.n_outputs_ == 1:
        y = np.ravel(y)
        y_std = np.ravel(y_std)

    return (y, y_std) if return_std else y

源码路径：sklearn/dummy.py - DummyRegressor.__sklearn_tags__()（第 507-512 行）

def __sklearn_tags__(self):
    tags = super().__sklearn_tags__()
    tags.input_tags.sparse = True
    tags.regressor_tags.poor_score = True
    tags.no_validation = True
    return tags

源码路径：sklearn/dummy.py - DummyRegressor.score()（第 514-540 行）

def score(self, X, y, sample_weight=None):
    """Return the coefficient of determination R^2 of the prediction."""
    if X is None:
        X = np.zeros(shape=(len(y), 1))
    return super().score(X, y, sample_weight)

点击查看 DummyRegressor 核心流程图

flowchart TD A[fit(X, y, sample_weight)] --> B[validate_data skip_check_array] B --> C[check_array 规范化 y 为 (n_samples, n_outputs)] C --> D[检查 y 非空] D --> E[检查长度一致性] E --> F[处理 sample_weight] F --> G{策略分发} G -- mean --> H[np.average 加权均值] G -- median --> I{sample_weight?} I -- 无 --> J[np.median] I -- 有 --> K[_weighted_percentile 50%] G -- quantile --> L[校验 quantile ∈ [0,1]] L --> M{sample_weight?} M -- 无 --> N[np.percentile q*100] M -- 有 --> O[_weighted_percentile q*100] G -- constant --> P[校验 constant 非空] P --> Q[check_array 接受稀疏 constant] Q --> R[校验形状匹配 n_outputs] H --> S[重塑 constant_ 为 (1, n_outputs)] K --> S N --> S O --> S R --> S S --> T[返回 self] U[predict(X, return_std)] --> V[check_is_fitted] V --> W[_num_samples 获取 n_samples] W --> X[np.full 广播 constant_ 为 (n_samples, n_outputs)] X --> Y[构造全零 y_std] Y --> Z{单输出?} Z -- 是 --> AA[ravel y 和 y_std] Z -- 否 --> BB[保持二维] AA --> CC{return_std?} BB --> CC CC -- 是 --> DD[返回 (y, y_std)] CC -- 否 --> EE[返回 y]

fit 中 check_array(y, ensure_2d=False) 允许一维或二维 y，随后 np.reshape(y, (-1, 1)) 统一为列向量矩阵。constant_ 最终 reshape(1, -1) 形成 (1, n_outputs) 行向量，predict 中 np.full((n_samples, n_outputs), self.constant_) 利用 numpy 广播机制将该行向量复制 n_samples 行，极其高效。return_std 返回全零数组，因为 DummyRegressor 无不确定性估计，但接口兼容性要求返回标准差（全零表示确定性预测）。score 允许 X=None 是为了方便仅用 y 评估基线性能，内部构造形状匹配的零矩阵代替。

20.12 设计中的取舍

为什么不用确定性矩阵代替随机矩阵？虽然可以使用哈达玛矩阵或其他确定性结构，但随机矩阵在理论上提供了更强的保证：Johnson-Lindenstrauss 引理适用于任意点集，而确定性结构通常需要数据满足特定性质（如正交性）才能保证距离保持。随机矩阵的另一个优势是不需要提前知道数据的特征结构。这种设计的 trade-off 是：随机投影牺牲了少量的距离精度（由 eps 参数控制），以换取显著的计算和存储效率提升。相比于 PCA 等需要特征值分解的确定性方法，随机投影的时间复杂度从 O(n²d) 降到 O(ndk)，其中 k 是目标维度，尤其在高维稀疏数据场景下优势明显。

对于 Dummy 估计器，设计取舍在于：完全放弃对输入特征 X 的建模，仅基于目标变量 y 的简单统计量进行预测。这种“懒惰”设计的优势是极简、极快、无需训练、天然支持多输出/稀疏目标/样本权重，且作为基线模型的语义清晰（poor_score=True 标记不应追求高分）。劣势显而易见：预测能力极弱，无法捕捉任何特征与目标的关系。no_validation=True 跳过输入校验进一步减少开销，但也意味着用户需自行保证数据质量。策略选择（most_frequent/prior/stratified/uniform/constant 等）提供了不同“懒惰程度”的基线，覆盖从确定性多数类到随机采样再到用户指定常数的谱系，满足不同基线对比需求。

20.13 动手练习

• 阅读 sklearn/random_projection.py 中以下函数/方法的实现：

  1. johnson_lindenstrauss_min_dim() (第 76-140 行) - JL引理维度计算

  2. _gaussian_random_matrix() (第 187-215 行) - 高斯矩阵生成

  3. _sparse_random_matrix() (第 217-285 行) - 稀疏矩阵生成

  4. BaseRandomProjection.fit() (第 307-345 行) - 统一拟合流程

  5. BaseRandomProjection.inverse_transform() (第 347-370 行) - 逆变换实现

  回答问题：

  • johnson_lindenstrauss_min_dim 如何处理数组输入的广播计算？分母为何设计为 eps^2/2 - eps^3/3？

  • _sparse_random_matrix 中 density='auto' 时的取值依据是什么？CSR矩阵构建的三个数组 data, indices, indptr 分别存什么？

  • BaseRandomProjection.fit 中 validate_data 为何设置 accept_sparse=['csr', 'csc'] 且 dtype=[np.float64, np.float32]？

  • inverse_transform 为何输出始终为密集数组？预计算伪逆与即时计算的内存/计算权衡是什么？

• 对比高斯与稀疏随机投影的 transform 实现

  • 阅读 sklearn/random_projection.py 中两个子类的 transform 方法：

    1. GaussianRandomProjection.transform() (第 347-360 行)

    2. SparseRandomProjection.transform() (第 480-495 行)

  回答问题：

  • 两者矩阵乘法的核心区别是什么？safe_sparse_dot 相比 @ 多了什么能力？

  • SparseRandomProjection 的 dense_output 参数如何影响不同输入(稀疏/稠密)下的输出格式？

  • 为何 GaussianRandomProjection 的 components_ 是稠密数组而 SparseRandomProjection 是 CSR 稀疏矩阵？这如何影响 inverse_transform 的实现？

• 深入 DummyClassifier 五大策略实现

  • 阅读 sklearn/dummy.py 中 DummyClassifier 的核心方法 (第 138-320 行)：

    1. fit() - 统计量计算与校验

    2. predict() - 五策略预测分发

    3. predict_proba() - 五策略概率输出

  回答问题：

  • fit 中 class_distribution 返回的三个值分别是什么？多输出时如何组织？

  • predict 中稀疏目标(sparse_output_=True) 时为何使用 _random_choice_csc 而密集目标用 numpy 操作？

  • stratified 策略下 predict 与 predict_proba 的关系是什么？为何 predict 要先调用 predict_proba？

  • uniform 策略为何不支持稀疏目标？代码中如何处理的？

  • constant 策略在 fit 和 predict/predict_proba 中如何处理多输出形状？

• 分析 DummyRegressor 加权分位数与多输出机制

  • 阅读 sklearn/dummy.py 中 DummyRegressor 的核心方法 (第 410-505 行)：

    1. fit() - 四策略统计量计算

    2. predict() - 预测生成与 return_std

  回答问题：

  • median 和 quantile 策略在有/无 sample_weight 时分别调用什么函数？_weighted_percentile 的作用是什么？

  • constant_ 属性为何统一重塑为 (1, n_outputs) 形状？predict 中如何利用该形状配合 np.full 实现广播？

  • return_std=True 时返回什么？为何设计为全零数组？

  • score 方法为何允许 X=None？内部如何处理？

• 设计随机投影与 Dummy 估计器的扩展测试

  • 基于 sklearn/tests/test_random_projection.py 和 sklearn/tests/test_dummy.py 的测试模式，设计以下新测试用例（仅描述测试思路，不写代码）：

    1. 验证 SparseRandomProjection 在 density=1.0 时与 GaussianRandomProjection 的嵌入质量等价性

    2. 验证 DummyClassifier 的 stratified 策略在多输出、样本权重、稀疏目标三重条件下的预测分布正确性

    3. 验证 DummyRegressor 的 quantile 策略在 quantile=0 和 quantile=1 边界下等价于 min/max 策略

    4. 验证随机投影的 inverse_transform 在 compute_inverse_components=False 时的数值稳定性(重复调用一致性)

    5. 验证 GaussianRandomProjection 和 SparseRandomProjection 在 n_components > n_features 时发出 DataDimensionalityWarning 但仍能工作

20.14 本章小结

这一章中我们学习了随机投影与 Dummy 估计器的实现原理。首先我们探讨了 Johnson-Lindenstrauss 引理作为随机投影的理论基础，理解了如何根据样本数和失真率计算最小投影维度；其次我们深入分析了高斯随机矩阵和稀疏随机矩阵的生成机制，掌握了它们的统计特性和数值实现；接着我们研究了 BaseRandomProjection 基类如何统一管理拟合、变换、逆变换和标签系统；然后我们分别考察了 GaussianRandomProjection 和 SparseRandomProjection 的具体实现，了解了它们在处理稠密和稀疏数据时的差异；最后我们学习了 DummyClassifier 和 DummyRegressor 如何通过简单的统计策略提供基线性能，支持多输出、稀疏目标和样本权重等高级特性。

同时把 summary 以表格方式总结：

| 概念 | 解释 |

|------|------|

| johnson_lindenstrauss_min_dim | JL引理工程实现，根据样本数和失真率计算最小投影维度，与原始特征数无关 |

| _gaussian_random_matrix | 稠密高斯随机矩阵生成，元素服从 N(0, 1/n_components)，列向量期望模长为1 |

| _sparse_random_matrix | 稀疏 Achlioptas/Li 随机矩阵生成，三点分布 {0, ±sqrt(s)/sqrt(k)}，CSR格式，密度可自动推断 1/sqrt(n_features) |

| BaseRandomProjection.fit | 统一拟合流程：自动推断n_components_、生成components_、可选预计算伪逆、设置_n_features_out |

| BaseRandomProjection.inverse_transform | 逆变换实现：预计算伪逆或即时pinv，稀疏components_自动转密集，输出始终密集 |

| GaussianRandomProjection.transform | 稠密矩阵乘法投影，保持输入dtype，接受稀疏/稠密输入 |

| SparseRandomProjection.transform | 稀疏感知矩阵乘法 safe_sparse_dot，dense_output控制输出格式，记录实际密度density_ |

| DummyClassifier.fit | 五策略统一拟合：class_distribution计算先验，支持样本权重、多输出、稀疏目标、常数策略合法性校验 |

| DummyClassifier.predict/predict_proba | 五策略预测分发：稀疏输出用_random_choice_csc，密集输出用numpy索引/多项采样/均匀采样 |

| DummyRegressor.fit | 四策略统计量计算：mean/median/quantile(加权分位数)/constant，形状统一为(1, n_outputs) |

| DummyRegressor.predict | np.full广播生成预测，支持return_std返回全零标准差，忽略X内容仅用行数 |

| sklearn_tags (Dummy) | poor_score=True标记基线模型，no_validation=True跳过输入校验，input_tags.sparse=True支持稀疏 |

| 测试验证体系 | 覆盖JL契约、嵌入质量、逆变换一致性、dtype保持、数值一致性、稀疏/密集互操作、策略分布正确性、样本权重、边界异常 |

下一章中，我们将学习 scikit-learn 概览 —— 数据科学家的“机器学习军火库”。

第 21 章 —— 随机投影与 Dummy 估计器 —— 品味“化繁为简的实用主义”

21.1 学习目标

• 难度：★★★☆☆（3/5）

• 预备知识：Python 基础、面向对象编程与 Markdown/代码阅读基础

• 理解 Johnson-Lindenstrauss 引理及其在随机投影中的理论保障

• 掌握 GaussianRandomProjection 和 SparseRandomProjection 的实现机制

• 了解 DummyClassifier 和 DummyRegressor 作为基线估计器的使用方式

• 能够运行并验证 random_projection 和 dummy 模块的基本功能

21.2 生活类比

想象一下，你正站在一个堆满高维数据的巨大仓库前，每一箱零件都有成百上千个特征，信息极其复杂。随机投影就像一台神奇的维度压缩机：它不需要拆解零件、不需要分析结构，只需用一套随机生成的压缩模具（投影矩阵）瞬间将高维零件压出低维影子。最令人惊叹的是，Johnson-Lindenstrauss 引理给出了数学承诺：只要压缩后的维度足够高，零件间的相对距离形状几乎不变，也就是说，压缩后的“影子”依然忠实保留了原始数据的几何结构。高斯投影模具材质细腻、压缩精准，适合稠密数据；稀疏投影模具则镂空轻量、压缩极快，适合高维稀疏或超大规模数据。而 Dummy 估计器则像是工厂里的基线质量标尺：不管原料多乱、工艺多复杂，它总能不看特征、只凭直觉（如众数、均值、分位数或常量）给出一个基线预测。这个标尺不追求精准，只负责划定“比乱猜强”的及格线——真模型若连这条线都跨不过，便毫无实战价值。

21.3 源码地图

sklearn/
├── random_projection.py
│   ├── class GaussianRandomProjection
│   │   ├── __init__(self, n_components='auto', eps=0.1, random_state=None)
│   │   ├── fit(self, X, y=None)
│   │   ├── transform(self, X)
│   │   └── components_
│   ├── class SparseRandomProjection
│   │   ├── __init__(self, n_components='auto', density='auto', eps=0.1, dense_output=False, random_state=None)
│   │   ├── fit(self, X, y=None)
│   │   ├── transform(self, X)
│   │   └── components_
│   ├── function johnson_lindenstrauss_min_dim
│   │   └── (n_samples, eps=0.1)
│   ├── function _gaussian_random_matrix
│   │   └── (n_components, n_features, random_state=None)
│   └── function _sparse_random_matrix
│       └── (n_components, n_features, density='auto', random_state=None)
├── dummy.py
│   ├── class DummyClassifier
│   │   ├── __init__(self, *, strategy="prior", random_state=None, constant=None)
│   │   ├── fit(self, X, y, sample_weight=None)
│   │   ├── predict(self, X)
│   │   ├── predict_proba(self, X)
│   │   ├── predict_log_proba(self, X)
│   │   ├── score(self, X, y, sample_weight=None)
│   │   └── classes_
│   ├── class DummyRegressor
│   │   ├── __init__(self, *, strategy="mean", quantile=None, random_state=None, constant=None)
│   │   ├── fit(self, X, y, sample_weight=None)
│   │   ├── predict(self, X)
│   │   ├── score(self, X, y, sample_weight=None)
│   │   └── constant_
│   └── __all__
└── tests/
├── test_random_projection.py
│   ├── test_gaussian_random_projection()
│   └── test_sparse_random_projection()
└── test_dummy.py
├── test_dummy_classifier()
├── test_dummy_regressor()
└── test_dummy_classifier_prior()

21.4 随机投影：化繁为简的“维度魔法师”

21.4.1 核心概念：Johnson-Lindenstrauss 引理

如何用随机矩阵安全降维而不破坏距离结构？Johnson-Lindenstrauss 引理提供理论保障：足够高的随机投影维度可近似保持点间欧氏距离。GaussianRandomProjection 使用服从正态分布的随机矩阵，适用于稠密数据场景；SparseRandomProjection 采用稀疏 Achlioptas 分布（值在 {-1, 0, +1} 中），实现更快的乘法和更低的存储开销。两种投影均在 fit 阶段生成 components_（投影矩阵），transform 仅做矩阵乘法，无需重训练。通过 johnson_lindenstrauss_min_dim 函数可计算达到误差容忍度所需的理论最小维度。

21.4.2 类型定义详解：投影矩阵生成函数

在深入类之前，我们先看两个核心工厂函数，它们是投影矩阵的“设计图纸”。

源码路径：sklearn/random_projection.py - _gaussian_random_matrix（第 130-159 行）

def _gaussian_random_matrix(n_components, n_features, random_state=None):
    """Generate a dense Gaussian random matrix.

    The components of the random matrix are drawn from

        N(0, 1.0 / n_components).

    Parameters
    ----------
    n_components : int,
        Dimensionality of the target projection space.

    n_features : int,
        Dimensionality of the original source space.

    random_state : int, RandomState instance or None, default=None
        Controls the pseudo random number generator used to generate the matrix
        at fit time.

    Returns
    -------
    components : ndarray of shape (n_components, n_features)
        The generated Gaussian random matrix.
    """
    _check_input_size(n_components, n_features)  # ① 检查维度合法性：n_components>0 且 n_features>0
    rng = check_random_state(random_state)       # ② 统一随机状态：接受 int/RandomState/None，返回 RandomState 实例
    components = rng.normal(                     # ③ 从 N(0, 1/n_components) 采样生成稠密矩阵
        loc=0.0, scale=1.0 / np.sqrt(n_components), size=(n_components, n_features)
    )                                             #    scale 缩放因子保证投影后方差守恒，是 JL 引理成立的数学前提
    return components                             # ④ 返回稠密 ndarray 矩阵，形状 (n_components, n_features)

概述：该函数生成服从高斯分布的稠密随机投影矩阵。核心在于 scale=1.0 / np.sqrt(n_components)，该缩放因子确保投影后各维度方差之和等于原始空间方差之和，满足 JL 引理对各向同性投影的数学要求。

源码路径：sklearn/random_projection.py - _sparse_random_matrix（第 161-239 行）

def _sparse_random_matrix(n_components, n_features, density="auto", random_state=None):
    """Generalized Achlioptas random sparse matrix for random projection.

    If we note :math:`s = 1 / density`, the components of the random matrix are
    drawn from:
      - -sqrt(s) / sqrt(n_components)   with probability 1 / 2s
      -  0                              with probability 1 - 1 / s
      - +sqrt(s) / sqrt(n_components)   with probability 1 / 2s
    """
    _check_input_size(n_components, n_features)  # ① 检查维度合法性
    density = _check_density(density, n_features) # ② 解析密度参数：'auto' 时取 1/sqrt(n_features)，这是理论最优稀疏度
    rng = check_random_state(random_state)        # ③ 统一随机状态

    if density == 1:                              # ④ 退化为稠密 ±1 矩阵：二项分布采样后映射到 {-1, +1}
        components = rng.binomial(1, 0.5, (n_components, n_features)) * 2 - 1
        return 1 / np.sqrt(n_components) * components

    else:
        # Generate location of non zero elements
        indices = []                              # ⑤ 存储所有行非零元素的列索引
        offset = 0
        indptr = [offset]                         # ⑥ CSR 格式的行指针数组，首元素为 0
        for _ in range(n_components):             # ⑦ 逐行生成非零模式
            n_nonzero_i = rng.binomial(n_features, density)  # 该行非零个数服从二项分布
            indices_i = sample_without_replacement(           # 无放回采样列位置，保证不重复
                n_features, n_nonzero_i, random_state=rng
            )
            indices.append(indices_i)
            offset += n_nonzero_i
            indptr.append(offset)

        indices = np.concatenate(indices)         # ⑧ 合并所有行索引为一维数组

        # Among non zero components the probability of the sign is 50%/50%
        data = rng.binomial(1, 0.5, size=np.size(indices)) * 2 - 1  # ⑨ 生成 ±1 符号

        # build the CSR structure by concatenating the rows
        components = sp.csr_matrix(               # ⑩ 构建 CSR 稀疏矩阵 (data, indices, indptr)
            (data, indices, indptr), shape=(n_components, n_features)
        )

        return np.sqrt(1 / density) / np.sqrt(n_components) * components  # ⑪ 缩放：保证非零元素方差期望为 1/n_components

概述：该函数构建广义 Achlioptas 稀疏随机矩阵。density='auto' 时默认 1/sqrt(n_features)，这是 Ping Li 等人证明的理论最优稀疏度。缩放因子 sqrt(1/density) / sqrt(n_components) 保证非零元素方差期望为 1/n_components，与高斯矩阵统一。通过 CSR 格式仅存储非零元素，大幅降低内存占用。

源码路径：sklearn/random_projection.py - johnson_lindenstrauss_min_dim（第 50-120 行）

@validate_params(
    {
        "n_samples": ["array-like", Interval(Real, 1, None, closed="left")],
        "eps": ["array-like", Interval(Real, 0, 1, closed="neither")],
    },
    prefer_skip_nested_validation=True,
)
def johnson_lindenstrauss_min_dim(n_samples, *, eps=0.1):
    """Find a 'safe' number of components to randomly project to."""
    eps = np.asarray(eps)                         # ① 转为数组支持批量计算
    n_samples = np.asarray(n_samples)             # ① 转为数组支持批量计算

    if np.any(eps <= 0.0) or np.any(eps >= 1):    # ② eps 必须在 (0, 1) 开区间
        raise ValueError("The JL bound is defined for eps in ]0, 1[, got %r" % eps)

    if np.any(n_samples <= 0):                    # ③ 样本数必须 > 0
        raise ValueError(
            "The JL bound is defined for n_samples greater than zero, got %r"
            % n_samples
        )

    denominator = (eps**2 / 2) - (eps**3 / 3)     # ④ 理论分母：eps²/2 - eps³/3，比常见 eps²/2 更紧致
    return (4 * np.log(n_samples) / denominator).astype(np.int64)  # ⑤ 计算公式并取整返回

概述：该函数直接实现 JL 引理核心公式 n_components >= 4 * log(n_samples) / (eps²/2 - eps³/3)。分母中的 eps³/3 项是更紧致的界，优于常见的 eps²/2 近似。函数支持 n_samples 和 eps 为数组的批量计算，便于超参数搜索时并行评估多个配置。

架构图：投影矩阵生成流程

graph TD A[输入参数 n_components, n_features, density, random_state] --> B{_gaussian_random_matrix} A --> C{_sparse_random_matrix} B --> D[check_random_state 统一随机状态] D --> E[rng.normal 采样 N(0, 1/n_components)] E --> F[返回稠密 ndarray] C --> G[check_random_state 统一随机状态] G --> H{density == 1?} H -- 是 --> I[二项分布生成稠密 ±1 矩阵] I --> J[缩放 1/sqrt(n_components)] H -- 否 --> K[_check_density 解析密度] K --> L[逐行二项分布决定非零数量] L --> M[sample_without_replacement 无放回采样列索引] M --> N[生成 CSR 结构: data, indices, indptr] N --> O[构建 sp.csr_matrix] O --> P[缩放 sqrt(1/density)/sqrt(n_components)] P --> Q[返回稀疏 CSR 矩阵]

21.4.3 基类设计：BaseRandomProjection

源码路径：sklearn/random_projection.py - BaseRandomProjection（第 241-334 行）

class BaseRandomProjection(
    ClassNamePrefixFeaturesOutMixin, TransformerMixin, BaseEstimator, metaclass=ABCMeta
):
    """Base class for random projections."""

    _parameter_constraints: dict = {
        "n_components": [
            Interval(Integral, 1, None, closed="left"),
            StrOptions({"auto"}),
        ],
        "eps": [Interval(Real, 0, None, closed="neither")],
        "compute_inverse_components": ["boolean"],
        "random_state": ["random_state"],
    }

    @abstractmethod
    def __init__(
        self,
        n_components="auto",
        *,
        eps=0.1,
        compute_inverse_components=False,
        random_state=None,
    ):
        self.n_components = n_components
        self.eps = eps
        self.compute_inverse_components = compute_inverse_components
        self.random_state = random_state

    @abstractmethod
    def _make_random_matrix(self, n_components, n_features):
        """Generate the random projection matrix."""

    def _compute_inverse_components(self):
        """Compute the pseudo-inverse of the (densified) components."""
        components = self.components_
        if sp.issparse(components):
            components = components.toarray()
        return linalg.pinv(components, check_finite=False)

    @_fit_context(prefer_skip_nested_validation=True)
    def fit(self, X, y=None):
        """Generate a sparse random projection matrix."""
        X = validate_data(
            self, X, accept_sparse=["csr", "csc"], dtype=[np.float64, np.float32]
        )

        n_samples, n_features = X.shape

        if self.n_components == "auto":              # ① 自动模式：用 JL 引理算维度
            self.n_components_ = johnson_lindenstrauss_min_dim(
                n_samples=n_samples, eps=self.eps
            )

            if self.n_components_ <= 0:              # ② 理论维度非正报错
                raise ValueError(
                    "eps=%f and n_samples=%d lead to a target dimension of "
                    "%d which is invalid" % (self.eps, n_samples, self.n_components_)
                )

            elif self.n_components_ > n_features:    # ③ 理论维度超原维度报错
                raise ValueError(
                    "eps=%f and n_samples=%d lead to a target dimension of "
                    "%d which is larger than the original space with "
                    "n_features=%d"
                    % (self.eps, n_samples, self.n_components_, n_features)
                )
        else:                                         # ④ 手动模式：用户指定维度
            if self.n_components > n_features:
                warnings.warn(
                    "The number of components is higher than the number of"
                    " features: n_features < n_components (%s < %s)."
                    "The dimensionality of the problem will not be reduced."
                    % (n_features, self.n_components),
                    DataDimensionalityWarning,
                )

            self.n_components_ = self.n_components

        # Generate a projection matrix of size [n_components, n_features]
        self.components_ = self._make_random_matrix(  # ⑤ 调用子类实现生成矩阵
            self.n_components_, n_features
        ).astype(X.dtype, copy=False)

        if self.compute_inverse_components:          # ⑥ 可选：计算伪逆用于逆变换
            self.inverse_components_ = self._compute_inverse_components()

        # Required by ClassNamePrefixFeaturesOutMixin.get_feature_names_out.
        self._n_features_out = self.n_components

        return self

    def inverse_transform(self, X):
        """Project data back to its original space."""
        check_is_fitted(self)

        X = check_array(X, dtype=[np.float64, np.float32], accept_sparse=("csr", "csc"))

        if self.compute_inverse_components:
            return X @ self.inverse_components_.T

        inverse_components = self._compute_inverse_components()
        return X @ inverse_components.T

    def __sklearn_tags__(self):
        tags = super().__sklearn_tags__()
        tags.transformer_tags.preserves_dtype = ["float64", "float32"]
        tags.input_tags.sparse = True
        return tags

概述：BaseRandomProjection 统一了随机投影的 fit/transform/inverse_transform 流程。fit 方法根据 n_components 参数决定目标维度：auto 模式调用 JL 引理计算并校验合法性，手动模式仅在超维时警告。核心步骤是调用抽象方法 _make_random_matrix 生成投影矩阵，并通过 .astype(X.dtype, copy=False) 与输入数据类型对齐，避免 transform 时隐式转换。可选计算伪逆支持 inverse_transform 近似重构原始数据。

架构图：BaseRandomProjection.fit 流程

graph TD A[fit(X, y=None)] --> B[validate_data 校验输入形状与类型] B --> C{n_components == 'auto'?} C -- 是 --> D[johnson_lindenstrauss_min_dim 计算 n_components_] D --> E[校验 n_components_ > 0] E --> F[校验 n_components_ <= n_features] C -- 否 --> G[使用用户指定 n_components] G --> H{n_components > n_features?} H -- 是 --> I[发出 DataDimensionalityWarning] H -- 否 --> J[直接赋值 n_components_] F --> K[_make_random_matrix 生成 components_] J --> K K --> L[.astype(X.dtype, copy=False) 类型对齐] L --> M{compute_inverse_components?} M -- 是 --> N[_compute_inverse_components 计算伪逆] M -- 否 --> O[跳过] N --> P[设置 _n_features_out] O --> P P --> Q[返回 self]

21.4.4 GaussianRandomProjection：稠密投影的标准实现

源码路径：sklearn/random_projection.py - GaussianRandomProjection.__init__（第 336-358 行）

    def __init__(
        self,
        n_components="auto",
        *,
        eps=0.1,
        compute_inverse_components=False,
        random_state=None,
    ):
        super().__init__(
            n_components=n_components,
            eps=eps,
            compute_inverse_components=compute_inverse_components,
            random_state=random_state,
        )

概述：构造函数仅将参数透传给基类 BaseRandomProjection，不包含额外逻辑。这种设计遵循模板方法模式，将维度计算、矩阵生成、类型对齐等通用逻辑下沉到基类，子类仅需实现 _make_random_matrix。

源码路径：sklearn/random_projection.py - GaussianRandomProjection._make_random_matrix（第 360-375 行）

    def _make_random_matrix(self, n_components, n_features):
        """Generate the random projection matrix."""
        random_state = check_random_state(self.random_state)  # ① 解析随机状态
        return _gaussian_random_matrix(                      # ② 委托工厂函数生成稠密高斯矩阵
            n_components, n_features, random_state=random_state
        )

概述：该方法解析随机状态后，直接委托给 _gaussian_random_matrix 生成形状为 (n_components, n_features) 的稠密高斯随机矩阵。矩阵元素服从 N(0, 1/n_components)，满足 JL 引理要求。

源码路径：sklearn/random_projection.py - GaussianRandomProjection.transform（第 377-394 行）

    def transform(self, X):
        """Project the data by using matrix product with the random matrix."""
        check_is_fitted(self)
        X = validate_data(
            self,
            X,
            accept_sparse=["csr", "csc"],
            reset=False,
            dtype=[np.float64, np.float32],
        )

        return X @ self.components_.T  # ③ 矩阵乘法：利用 NumPy BLAS 加速，输出稠密数组

概述：transform 仅执行矩阵乘法 X @ components_.T，利用 NumPy 底层 BLAS 库加速稠密矩阵运算。输入支持稀疏 CSR/CSC 格式，但输出始终为稠密数组，适合中小规模稠密数据场景。

架构图：GaussianRandomProjection 核心流程

graph LR A[GaussianRandomProjection] --> B[继承 BaseRandomProjection] B --> C[__init__: 透传参数给基类] B --> D[_make_random_matrix: 调用 _gaussian_random_matrix] D --> E[生成稠密 ndarray N(0, 1/n_components)] B --> F[fit: 复用基类逻辑] B --> G[transform: X @ components_.T] G --> H[利用 BLAS 加速稠密矩阵乘法] H --> I[输出稠密 ndarray]

21.4.5 SparseRandomProjection：稀疏投影的工程智慧

源码路径：sklearn/random_projection.py - SparseRandomProjection.__init__（第 396-420 行）

    def __init__(
        self,
        n_components="auto",
        *,
        density="auto",
        eps=0.1,
        dense_output=False,
        compute_inverse_components=False,
        random_state=None,
    ):
        super().__init__(
            n_components=n_components,
            eps=eps,
            compute_inverse_components=compute_inverse_components,
            random_state=random_state,
        )

        self.dense_output = dense_output
        self.density = density

概述：构造函数新增 density（稀疏度）和 dense_output（输出格式控制）两个参数，其余透传基类。density='auto' 时将在 _make_random_matrix 中按 1/sqrt(n_features) 计算。

源码路径：sklearn/random_projection.py - SparseRandomProjection._make_random_matrix（第 422-436 行）

    def _make_random_matrix(self, n_components, n_features):
        """Generate the random projection matrix"""
        random_state = check_random_state(self.random_state)
        self.density_ = _check_density(self.density, n_features)
        return _sparse_random_matrix(
            n_components, n_features, density=self.density_, random_state=random_state
        )

概述：解析随机状态，调用 _check_density 确定具体稀疏度并保存为 density_ 属性，随后委托 _sparse_random_matrix 生成 CSR 格式稀疏矩阵。

源码路径：sklearn/random_projection.py - SparseRandomProjection.transform（第 438-456 行）

    def transform(self, X):
        """Project the data by using matrix product with the random matrix."""
        check_is_fitted(self)
        X = validate_data(
            self,
            X,
            accept_sparse=["csr", "csc"],
            reset=False,
            dtype=[np.float64, np.float32],
        )

        return safe_sparse_dot(X, self.components_.T, dense_output=self.dense_output)

概述：使用 safe_sparse_dot 执行稀疏感知矩阵乘法。该函数自动处理“稀疏×稀疏”、“稀疏×稠密”、“稠密×稀疏”三种组合，并根据 dense_output 参数决定返回稠密还是稀疏结果。当输入稀疏且 dense_output=False 时输出保持稀疏，节省内存；投影维度极小（非零元极少）时设为 True 反而更高效。

核心差异对比：

下表对比了两种投影的关键特性差异，帮助根据数据特征选择合适的投影器：

| 特性 | GaussianRandomProjection | SparseRandomProjection |

|------|--------------------------|------------------------|

| 矩阵类型 | 稠密 ndarray | 稀疏 csr_matrix |

| 生成函数 | _gaussian_random_matrix | _sparse_random_matrix |

| 乘法方式 | X @ components_.T (BLAS) | safe_sparse_dot (稀疏感知) |

| 密度参数 | 无 | density='auto' (默认 1/sqrt(n_features)) |

| 输出控制 | 总是稠密 | dense_output 控制稀疏/稠密输出 |

| 内存占用 | O(n_components × n_features) | O(nnz) ≈ O(density × n_components × n_features) |

| 适用场景 | 中小规模稠密数据 | 高维稀疏数据、超大规模稠密数据 |

架构图：SparseRandomProjection 核心流程

graph LR A[SparseRandomProjection] --> B[继承 BaseRandomProjection] B --> C[__init__: 新增 density/dense_output 参数] B --> D[_make_random_matrix: 调用 _sparse_random_matrix] D --> E[_check_density 确定稀疏度] E --> F[生成 CSR 稀疏矩阵 Achlioptas 分布] B --> G[fit: 复用基类逻辑] B --> H[transform: safe_sparse_dot] H --> I{dense_output?} I -- True --> J[返回稠密 ndarray] I -- False --> K[输入稀疏则返回稀疏 CSC/CSR]

21.4.6 随机投影数据流图

graph TD A[输入数据 X (n_samples, n_features)] --> B{fit 阶段} B --> C[n_components='auto'?] C -- 是 --> D[调用 johnson_lindenstrauss_min_dim 计算 n_components_] C -- 否 --> E[使用用户指定 n_components] D --> F[校验 n_components_ <= n_features] E --> F F --> G[调用 _make_random_matrix 生成 components_] G --> H[可选: 计算 inverse_components_] H --> I[fit 完成] I --> J[transform 阶段] J --> K[validate_data 校验输入] K --> L{矩阵乘法} L -- Gaussian --> M[X @ components_.T] L -- Sparse --> N[safe_sparse_dot(X, components_.T)] M --> O[输出 X_new (n_samples, n_components)] N --> O

21.5 Dummy 估计器：建立模型评估的“基线坐标系”

21.5.1 核心概念：为什么需要基线？

DummyClassifier 和 DummyRegressor 通过极简规则（如众数、均值、分位数或常量）生成预测，不依赖任何特征。支持多种 strategy：分类器含 'prior'、'stratified'、'uniform'、'constant'；回归器含 'mean'、'median'、'quantile'、'constant'。fit 过程中仅学习一个常量统计量（如训练标签的众数或均值），predict 忽略输入 X 直接返回该常量。score 方法允许使用任意评估函数（如 accuracy、f1、r2），为真实模型提供客观的下限参考。即使在特征全噪声或缺失的极端情况下，Dummy 估计器仍能提供稳定的基线表现。

21.5.2 DummyClassifier：分类基线的多策略实现

源码路径：sklearn/dummy.py - DummyClassifier.__init__（第 70-78 行）

    def __init__(self, *, strategy="prior", random_state=None, constant=None):
        self.strategy = strategy
        self.random_state = random_state
        self.constant = constant

概述：构造函数接收三个仅关键字参数：strategy 指定基线策略（默认 "prior"），random_state 控制随机策略的可复现性，constant 仅在 "constant" 策略下使用。参数约束由 _parameter_constraints 定义。

源码路径：sklearn/dummy.py - DummyClassifier.fit（第 80-170 行）

    @_fit_context(prefer_skip_nested_validation=True)
    def fit(self, X, y, sample_weight=None):
        """Fit the baseline classifier."""
        validate_data(self, X, skip_check_array=True)  # ① 仅校验 X 形状，不检查数值

        self._strategy = self.strategy

        if self._strategy == "uniform" and sp.issparse(y):  # ② uniform 策略不支持稀疏 y，转稠密并警告
            y = y.toarray()
            warnings.warn(...)

        self.sparse_output_ = sp.issparse(y)  # ③ 记录输出格式：稀疏输入则稀疏输出

        if not self.sparse_output_:
            y = np.asarray(y)
            y = np.atleast_1d(y)

        if y.ndim == 1:
            y = np.reshape(y, (-1, 1))        # ④ 统一为二维 (n_samples, n_outputs) 支持多输出

        self.n_outputs_ = y.shape[1]

        check_consistent_length(X, y)

        if sample_weight is not None:
            sample_weight = _check_sample_weight(sample_weight, X)

        if self._strategy == "constant":        # ⑤ constant 策略校验常量合法性
            if self.constant is None:
                raise ValueError("Constant target value has to be specified ...")
            else:
                constant = np.reshape(np.atleast_1d(self.constant), (-1, 1))
                if constant.shape[0] != self.n_outputs_:
                    raise ValueError("Constant target value should have shape (%d, 1)." % self.n_outputs_)

        # ⑥ 核心：计算类别分布（支持样本权重）
        (self.classes_, self.n_classes_, self.class_prior_) = class_distribution(
            y, sample_weight
        )

        if self._strategy == "constant":        # ⑦ 校验常量是否在训练类别中
            for k in range(self.n_outputs_):
                if not any(constant[k][0] == c for c in self.classes_[k]):
                    raise ValueError(
                        "The constant target value must be present in "
                        "the training data. You provided constant={}. "
                        "Possible values are: {}.".format(
                            self.constant, self.classes_[k].tolist()
                        )
                    )

        if self.n_outputs_ == 1:                # ⑧ 单输出时扁平化属性
            self.n_classes_ = self.n_classes_[0]
            self.classes_ = self.classes_[0]
            self.class_prior_ = self.class_prior_[0]

        return self

概述：fit 方法仅校验 X 形状，核心工作是调用 class_distribution 计算加权类别分布（classes_、n_classes_、class_prior_）。class_prior_ 即经验先验概率，支持样本权重。根据策略校验常量合法性，单输出时扁平化属性便于后续使用。

架构图：DummyClassifier.fit 核心逻辑

graph TD A[fit(X, y, sample_weight)] --> B[validate_data 仅校验 X 形状] B --> C{strategy == 'uniform' 且 y 稀疏?} C -- 是 --> D[y.toarray() 转稠密并警告] C -- 否 --> E[记录 sparse_output_] E --> F[统一 y 为二维 (n_samples, n_outputs)] F --> G[check_consistent_length] G --> H{strategy == 'constant'?} H -- 是 --> I[校验 constant 非空且形状匹配] H -- 否 --> J[class_distribution 计算类别分布] I --> J J --> K[保存 classes_, n_classes_, class_prior_] K --> L{单输出?} L -- 是 --> M[扁平化属性为一维] L -- 否 --> N[保持列表结构] M --> O[返回 self] N --> O

源码路径：sklearn/dummy.py - DummyClassifier.predict（第 172-240 行）

    def predict(self, X):
        """Perform classification on test vectors X."""
        check_is_fitted(self)

        n_samples = _num_samples(X)
        rs = check_random_state(self.random_state)

        n_classes_ = self.n_classes_
        classes_ = self.classes_
        class_prior_ = self.class_prior_
        constant = self.constant
        if self.n_outputs_ == 1:
            # Get same type even for self.n_outputs_ == 1
            n_classes_ = [n_classes_]
            classes_ = [classes_]
            class_prior_ = [class_prior_]
            constant = [constant]
        # Compute probability only once
        if self._strategy == "stratified":
            proba = self.predict_proba(X)
            if self.n_outputs_ == 1:
                proba = [proba]

        if self.sparse_output_:                 # ① 稀疏输出路径
            class_prob = None
            if self._strategy in ("most_frequent", "prior"):
                classes_ = [np.array([cp.argmax()]) for cp in class_prior_]

            elif self._strategy == "stratified":
                class_prob = class_prior_

            elif self._strategy == "uniform":
                raise ValueError("Sparse target prediction is not supported with the uniform strategy")

            elif self._strategy == "constant":
                classes_ = [np.array([c]) for c in constant]

            y = _random_choice_csc(n_samples, classes_, class_prob, self.random_state)
        else:                                   # ② 稠密输出路径
            if self._strategy in ("most_frequent", "prior"):
                y = np.tile(
                    [
                        classes_[k][class_prior_[k].argmax()]
                        for k in range(self.n_outputs_)
                    ],
                    [n_samples, 1],
                )

            elif self._strategy == "stratified":
                y = np.vstack(
                    [
                        classes_[k][proba[k].argmax(axis=1)]
                        for k in range(self.n_outputs_)
                    ]
                ).T

            elif self._strategy == "uniform":
                ret = [
                    classes_[k][rs.randint(n_classes_[k], size=n_samples)]
                    for k in range(self.n_outputs_)
                ]
                y = np.vstack(ret).T

            elif self._strategy == "constant":
                y = np.tile(self.constant, (n_samples, 1))

            if self.n_outputs_ == 1:
                y = np.ravel(y)

        return y