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

Sklearn 源码解析(书)v1.0(四)

RANSACRegressor 使用 _dynamic_max_trials 根据当前内点比例动态调整 max_trials。相比固定迭代次数,动态策略在内点比例高时大幅减少计算,在内点比例低时自动增加尝试次数以保证置信度。stop_probability 控制的是 “至少采样到一次全内点子集的概率”,这是 RANSAC 理论保证的核心。max_skips 限制无效采样的累积,防止因 is_data_valid/is_model_valid 过于严格导致无限循环。

TheilSen 全组合与随机采样的精度/计算权衡

TheilSenRegressorC(n_samples, n_subsamples) <= max_subpopulation 时使用 全组合,保证空间中位数的统计精度;超出时退化为 随机子采样,以 max_subpopulation 控制计算与内存上限。n_subsamples 默认取 n_features + 1(含截距),这是达到 最大破坏点 的最小子集大小。_spatial_median 的 Weiszfeld 迭代复杂度为 O(n_subpopulation * n_features),通常远小于组合数。

SparseCoefMixin 稀疏/稠密转换的内存/性能权衡

SparseCoefMixin 提供 densify/sparsify 双向转换。稀疏格式(CSR)在 高维稀疏系数(如文本分类)下节省内存并加速 X @ coef;稠密格式在 中低维或系数非稀疏 时因缓存友好且无间接寻址而更快。用户可根据下游任务(预测、持久化、解释性分析)灵活切换,但需注意频繁转换的开销。

10.11 动手练习

  • 追踪 LinearRegression 的三种求解路径

    • 阅读 LinearRegression.fit() 中的 if self.positiveelif sp.issparse(X)else 分支。

    • 回答:

      • 为什么稀疏 X 不能使用 positive=True

      • lsqratol / btoltol 参数有什么关系?

      • 为什么稀疏 X 的 copy_X 恒为 False

  • 逆向 _preprocess_datasample_weight 缩放

    • 阅读 _preprocess_data_rescale_data

    • 回答:

      • 稀疏 X 为什么只能记录 X_offset 而不直接中心化?

      • sample_weight_sqrt 在哪些求解器中会被使用?

      • 为什么 _preprocess_data 返回的 X_scale 恒为全 1?

  • 解析 LinearModelLoss 与贝叶斯回归的参数更新

    • 浏览 LinearModelLoss_bayes.py 中的 BayesianRidge.fit_update_coef__log_marginal_likelihood

    • 回答:

      • LinearModelLoss 如何处理 multiclass 系数布局?

      • BayesianRidge 何时需要 full_matrices=True 的完整 SVD?

      • ARDRegressionn_samples < n_features 时使用了什么矩阵求逆技巧?

  • 对比四种鲁棒回归器的核心机制

    • 阅读 _huber.py_quantile.py_ransac.py_theil_sen.py 中的 fit 实现。

    • 回答:

      • HuberRegressorbounds 为什么对最后一个参数(sigma)设置了下界?

      • QuantileRegressor 为何要对 alpha 乘以 sum(sample_weight)

      • RANSACstop_probabilitymax_trials 如何交互决定迭代上限?

      • TheilSenRegressor 何时需要采用随机子采样而非全组合?

10.12 本章小结

这一章我们深入探索了 线性模型的内部引擎

  • 通过 LinearRegression.fit 的三条求解路径,掌握了 稠密、稀疏、正系数 三种情形的数值实现与容错设计。

  • _preprocess_data_rescale_data加权中心化 提供统一框架,特别注意了稀疏数据的 不破坏稀疏性 策略。

  • make_dataset 为坐标下降法提供统一的数据抽象与稀疏截距衰减机制。

  • _pre_fit_check_precomputed_gram_matrix 展示了 Gram 矩阵预计算 的细粒度校验与重复利用机制。

  • LinearClassifierMixinSparseCoefMixin线性分类稀疏系数 提供了统一且高效的 API。

  • LinearModelLoss 将各种基损失统一为 loss / gradient / Hessian 接口,核心计算通过 sandwich_dot 实现高效矩阵-对角-矩阵乘积。

  • BayesianRidgeARDRegression 通过 SVD / Woodbury 以及 证据最大化 完成贝叶斯推断,获得了 后验均值、协方差超参数 的自适应更新。

  • 四种鲁棒回归器(HuberRegressorQuantileRegressorRANSACRegressorTheilSenRegressor)分别在 损失函数形态线性规划随机抽样共识空间中位数 方面提供了对离群点的强大抵抗能力。

接下来,在第11章 支持向量机 —— 揭开“最大间隔分类器”的面纱 中,我们将从底层 libsvm/liblinear 的 C/C++ 内核到 Python 包装层,完整剖析 核函数、拉格朗日乘子、SMO 求解器 以及 稀疏/稠密数据的统一接口,帮助你掌握 SVM 的数学与实现细节。

第 11 章 —— 支持向量机 —— 揭开“最大间隔分类器”的面纱

11.1 学习目标

  • 难度:★★★☆☆(3/5)

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

  • 理解 Python 层 BaseLibSVM/BaseSVC 如何统一管理核函数、支持向量、决策函数与稀疏/稠密数据的分发

  • 掌握 libsvm C++ 内核中 SMO 求解器的工作集选择、收缩启发式与核缓存机制

  • 掌握 liblinear C++ 内核中坐标下降与 TRON 信任域 Newton 法的算法实现与求解器选择逻辑

  • 理解 Cython 桥接层如何实现稠密/稀疏数据的内存视图映射、CSR 直通与内存管理

  • 理解跨平台随机数生成器(Mersenne Twister + Lemire 后处理)如何保证 SVM 优化的确定性

  • 理解 l1_min_c 边界计算如何指导 L1 正则化 SVM 的参数搜索

11.2 生活类比

想象 SVM 模块是一个百万级物流中心的智能分拣系统:

  • BaseLibSVM = 分拣中心的中控台,负责核函数(运输方式)、稀疏/稠密(货物形态)、决策函数(配送路线)的统一调度

  • libsvm 的 SMO 求解器 = 精挑细选的质检员,每次只挑出两个“问题最严重”的包裹(工作集选择)进行调整

  • Cache 核缓存 = 常用工具就近摆放的货架,避免每次都去远处仓库取

  • 收缩启发式 = 将已经“确定不再变化”的包裹暂存到旁边,只处理活跃部分

  • liblinear 坐标下降 = 逐个调整每件货品的位置(按坐标轴逐个优化),每次只动一个维度

  • TRON 信任域 = 使用“大步流星”的策略,先用二次模型估计整批移动方向,如果估计不准就缩小步长重试

  • Cython 桥接层 = 会说两种语言的“翻译官”,将 Python 的 NumPy 包裹直接搬运到 C 层的货车,零拷贝直通

  • Mersenne Twister 随机数 = 全平台统一的“骰子”,确保 Windows 和 Linux 上优化路径完全一致

就像物流中心的每一个环节都需要精密的协作——从收货验货(数据验证)到分拣打包(SMO/坐标下降),再到装车配送(Cython 桥接预测),SVM 模块正是这样一套环环相扣的高效系统。这一类比将贯穿全章,帮助你在阅读源码时建立直观的心智模型。

11.3 源码地图

graph TD subgraph Python层 A[sklearn/svm/_base.py] --> A1[BaseLibSVM] A --> A2[BaseSVC] A --> A3[_fit_liblinear] A --> A4[_get_liblinear_solver_type] B[sklearn/svm/_classes.py] --> B1[LinearSVC/LinearSVR] B --> B2[SVC/NuSVC/SVR/NuSVR/OneClassSVM] C[sklearn/svm/_bounds.py] --> C1[l1_min_c] D[sklearn/svm/_newrand.pyx] --> D1[set_seed_wrap/bounded_rand_int_wrap] end subgraph Cython桥接层 E[sklearn/svm/_libsvm.pyx] --> E1[fit/predict/decision_function] F[sklearn/svm/_libsvm_sparse.pyx] --> F1[稀疏训练/预测] G[sklearn/svm/_liblinear.pyx] --> G1[train_wrap] end subgraph C辅助层 H[sklearn/svm/src/libsvm/libsvm_helper.c] I[sklearn/svm/src/libsvm/libsvm_sparse_helper.c] J[sklearn/svm/src/liblinear/liblinear_helper.c] end subgraph libsvm核心 K[sklearn/svm/src/libsvm/svm.cpp] --> K1[Cache/Kernel/Solver] K --> K2[SVC_Q/ONE_CLASS_Q/SVR_Q] K --> K3[solve_c_svc/solve_nu_svc/...] end subgraph liblinear核心 L[sklearn/svm/src/liblinear/linear.cpp] --> L1[l2r_lr_fun/l2r_l2_svc_fun] L --> L2[Solver_MCSVM_CS] L --> L3[solve_l2r_l1l2_svc/...] M[sklearn/svm/src/liblinear/tron.cpp] --> M1[TRON/tron/trcg] end subgraph 随机数 N[sklearn/svm/src/newrand/newrand.h] --> N1[mt19937/Lemire] end A --> E A --> F A --> G E --> H F --> I G --> J E --> K F --> K G --> L G --> M K --> N L --> N

11.4 Python 层基座:BaseLibSVM 的“总调度室”

为什么需要 BaseLibSVM 抽象基类?

正如物流中心需要一个中控台统一调度运输方式、货物形态和配送路线,BaseLibSVM 作为 SVM 家族的抽象基类,统一管理 C-SVC、Nu-SVC、SVR、OneClassSVM 共有的核函数、支持向量和决策函数逻辑。它通过 _sparse_kernels 列表和 _sparse 标志实现稠密/稀疏数据的自动分发,fit() 方法完成从数据验证到 libsvm 调用的完整流水线。这种设计避免了在每个具体 SVM 变体中重复实现核函数计算、预测分发等通用逻辑,体现了“中控台统一调度”的架构思想。

fit() 的“安全检查站”

就像货物入库前必须经过安检和验货,BaseLibSVM.fit() 首先进行严格的输入验证。它使用 validate_data 强制 float64、C 连续和 CSR 格式要求,对 precomputed 核拒绝稀疏输入,对 callable 核跳过数值验证。通过 _validate_targets 钩子委托子类处理标签和类别权重,体现了模板方法模式。数据验证阶段就像物流中心的“收货验货环节”,确保后续处理的数据质量。

gamma 参数的三态解析

gamma 参数支持三种模式,就像物流中心根据货物类型选择不同的运输方式:

  • 'scale':基于 X 的方差自动缩放(处理稀疏时用 E[X²]-E[X]² 计算方差),相当于根据货物密度动态调整运输策略

  • 'auto':使用 1/n_features,相当于按固定规格标准化运输

  • 浮点值:直接使用用户指定的 gamma,相当于完全自定义运输参数

概率参数弃用迁移

probability 参数在 1.9 版本被标记为 FutureWarning,建议改用 CalibratedClassifierCV(ensemble=False)。弃用字符串 'deprecated' 作为默认值,保留旧行为但引导新 API。这就像物流中心逐步淘汰旧式分拣设备,引导客户使用新的智能分拣系统,同时兼容过渡期的混合使用。

fit() 的“收尾工作”

拟合结束后保存 shape_fit_ 用于后续预测验证,二分类时翻转 intercept_dual_coef_ 符号(内部保存 _intercept__dual_coef_),检查对偶系数和截距的有限性防止大数值导致 NaN。这相当于分拣完成后的“出库复核和归档”,确保模型状态一致可用。

源码路径:sklearn/svm/_base.py - BaseLibSVM.fit()(117-250行)

    @_fit_context(prefer_skip_nested_validation=True)
    def fit(self, X, y, sample_weight=None):
        """Fit the SVM model according to the given training data.

        Parameters
        ----------
        X : {array-like, sparse matrix} of shape (n_samples, n_features) \
                or (n_samples, n_samples)
            Training vectors, where `n_samples` is the number of samples
            and `n_features` is the number of features.
            For kernel="precomputed", the expected shape of X is
            (n_samples, n_samples).

        y : array-like of shape (n_samples,)
            Target values (class labels in classification, real numbers in
            regression).

        sample_weight : array-like of shape (n_samples,), default=None
            Per-sample weights. Rescale C per sample. Higher weights
            force the classifier to put more emphasis on these points.

        Returns
        -------
        self : object
            Fitted estimator.

        Notes
        -----
        If X and y are not C-ordered and contiguous arrays of np.float64 and
        X is not a scipy.sparse.csr_matrix, X and/or y may be copied.

        If X is a dense array, then the other methods will not support sparse
        matrices as input.
        """
        rnd = check_random_state(self.random_state)  # 创建随机数生成器,用于后续生成随机种子

        sparse = sp.issparse(X)  # 检测输入是否为稀疏矩阵
        if sparse and self.kernel == "precomputed":
            raise TypeError("Sparse precomputed kernels are not supported.")  # precomputed 核拒绝稀疏输入
        self._sparse = sparse and not callable(self.kernel)  # 设置稀疏标志:稀疏输入且非 callable 核时为 True

        if callable(self.kernel):
            check_consistent_length(X, y)  # callable 核仅检查长度一致性
        else:
            X, y = validate_data(  # 强制验证:float64、C 连续、接受 CSR、拒绝大稀疏
                self,
                X,
                y,
                dtype=np.float64,
                order="C",
                accept_sparse="csr",
                accept_large_sparse=False,
            )

        y = self._validate_targets(y)  # 委托子类验证目标(分类/回归/单类)

        sample_weight = np.asarray(  # 标准化样本权重为 float64 数组
            [] if sample_weight is None else sample_weight, dtype=np.float64
        )
        solver_type = LIBSVM_IMPL.index(self._impl)  # 将 _impl 字符串映射为 libsvm 整数编码

        # TODO(1.11): remove probability
        self._effective_probability = self.probability  # 概率参数有效值
        if self._impl in ["c_svc", "nu_svc"]:
            if self._impl == "nu_scv":
                est_dep = "NuSVC"
            else:
                est_dep = "SVC"
            if self.probability != "deprecated":  # 1.9 版本起弃用 probability 参数
                warnings.warn(
                    f"The `probability` parameter was deprecated in 1.9 and "
                    f"will be removed in version 1.11. "
                    f"Use `CalibratedClassifierCV({est_dep}(), ensemble=False)` "
                    f"instead of `{est_dep}(probability=True)`",
                    FutureWarning,
                )
            else:
                self._effective_probability = False

        # input validation
        n_samples = _num_samples(X)
        if solver_type != 2 and n_samples != y.shape[0]:  # one_class 除外需检查样本数匹配
            raise ValueError(
                "X and y have incompatible shapes.\n"
                + "X has %s samples, but y has %s." % (n_samples, y.shape[0])
            )

        if self.kernel == "precomputed" and n_samples != X.shape[1]:  # precomputed 核必须是方阵
            raise ValueError(
                "Precomputed matrix must be a square matrix."
                " Input is a {}x{} matrix.".format(X.shape[0], X.shape[1])
            )

        if sample_weight.shape[0] > 0 and sample_weight.shape[0] != n_samples:  # 样本权重维度检查
            raise ValueError(
                "sample_weight and X have incompatible shapes: "
                "%r vs %r\n"
                "Note: Sparse matrices cannot be indexed w/"
                "boolean masks (use `indices=True` in CV)."
                % (sample_weight.shape, X.shape)
            )

        kernel = "precomputed" if callable(self.kernel) else self.kernel  # 统一核名称

        if kernel == "precomputed":
            # unused but needs to be a float for cython code that ignores
            # it anyway
            self._gamma = 0.0  # precomputed 核不需要 gamma
        elif isinstance(self.gamma, str):
            if self.gamma == "scale":
                # var = E[X^2] - E[X]^2 if sparse
                X_var = (X.multiply(X)).mean() - (X.mean()) ** 2 if sparse else X.var()  # 稀疏用 E[X²]-E[X]²,稠密用 var()
                self._gamma = 1.0 / (X.shape[1] * X_var) if X_var != 0 else 1.0  # gamma = 1/(n_features * var)
            elif self.gamma == "auto":
                self._gamma = 1.0 / X.shape[1]  # gamma = 1/n_features
        elif isinstance(self.gamma, Real):
            self._gamma = self.gamma  # 用户指定浮点值直接使用

        fit = self._sparse_fit if self._sparse else self._dense_fit  # 根据 _sparse 分发训练路径
        if self.verbose:
            print("[LibSVM]", end="")

        seed = rnd.randint(np.iinfo("i").max)  # 生成随机种子传给底层
        fit(X, y, sample_weight, solver_type, kernel, random_seed=seed)
        # see comment on the other call to np.iinfo in this file

        self.shape_fit_ = X.shape if hasattr(X, "shape") else (n_samples,)  # 保存训练形状供预测验证

        # In binary case, we need to flip the sign of coef, intercept and
        # decision function. Use self._intercept_ and self._dual_coef_
        # internally.
        self._intercept_ = self.intercept_.copy()  # 保存内部副本
        self._dual_coef_ = self.dual_coef_
        if self._impl in ["c_svc", "nu_svc"] and len(self.classes_) == 2:  # 二分类时翻转符号
            self.intercept_ *= -1
            self.dual_coef_ = -self.dual_coef_

        dual_coef = self._dual_coef_.data if self._sparse else self._dual_coef_  # 稀疏取 data 数组
        intercept_finiteness = np.isfinite(self._intercept_).all()  # 检查截距有限性
        dual_coef_finiteness = np.isfinite(dual_coef).all()  # 检查对偶系数有限性
        if not (intercept_finiteness and dual_coef_finiteness):
            raise ValueError(
                "The dual coefficients or intercepts are not finite."
                " The input data may contain large values and need to be"
                " preprocessed."
            )

        # Since, in the case of SVC and NuSVC, the number of models optimized by
        # libSVM could be greater than one (depending on the input), `n_iter_`
        # stores an ndarray.
        # For the other sub-classes (SVR, NuSVR, and OneClassSVM), the number of
        # models optimized by libSVM is always one, so `n_iter_` stores an
        # integer.
        if self._impl in ["c_svc", "nu_svc"]:
            self.n_iter_ = self._num_iter  # 多分类存数组
        else:
            self.n_iter_ = self._num_iter.item()  # 回归/单类存标量

        return self

代码解析:

这段代码实现了 SVM 拟合的完整流水线:

  1. 随机状态初始化(第 1 行):创建随机数生成器用于后续的随机种子生成

  2. 稀疏性检测与分发标记(第 4-7 行):检测输入是否为稀疏矩阵,precomputed 核拒绝稀疏输入,设置 _sparse 标志

  3. 数据验证(第 9-17 行):callable 核仅检查长度一致性,否则强制 float64、C 连续、CSR 格式

  4. 目标验证委托(第 19 行):调用子类重写的 _validate_targets 处理标签

  5. 样本权重标准化(第 21-23 行):转为 float64 数组,空则置为空数组

  6. 求解器类型映射(第 25 行):将 _impl 字符串映射为 libsvm 整数编码

  7. 概率参数弃用处理(第 27-39 行):1.9 版本起弃用 probability,引导使用 CalibratedClassifierCV

  8. 形状兼容性检查(第 41-54 行):验证 X/y 样本数、precomputed 核方阵、sample_weight 维度

  9. gamma 三态解析(第 56-70 行):scale/auto/float 三种模式,稀疏时用 E[X²]-E[X]² 计算方差

  10. 分发训练(第 72-76 行):根据 _sparse 选择稠密/稀疏训练路径,生成随机种子调用底层

  11. 后处理(第 78-103 行):保存训练形状、二分类符号翻转、有限性检查、n_iter_ 赋值


源码路径:sklearn/svm/_base.py - BaseLibSVM.__init__()(70-116行)

    def __init__(
        self,
        kernel,
        degree,
        gamma,
        coef0,
        tol,
        C,
        nu,
        epsilon,
        shrinking,
        probability,
        cache_size,
        class_weight,
        verbose,
        max_iter,
        random_state,
    ):
        if self._impl not in LIBSVM_IMPL:  # 校验 _impl 是否在合法列表中
            raise ValueError(
                "impl should be one of %s, %s was given" % (LIBSVM_IMPL, self._impl)
            )

        self.kernel = kernel
        self.degree = degree
        self.gamma = gamma
        self.coef0 = coef0
        self.tol = tol
        self.C = C
        self.nu = nu
        self.epsilon = epsilon
        self.shrinking = shrinking
        self.probability = probability
        self.cache_size = cache_size
        self.class_weight = class_weight
        self.verbose = verbose
        self.max_iter = max_iter
        self.random_state = random_state

代码解析:

构造函数完成两件核心工作:

  1. 实现类型校验(第 3-6 行):检查 _impl 属性是否在合法列表 LIBSVM_IMPL 中(c_svc/nu_svc/one_class/epsilon_svr/nu_svr),这是子类必须设置的类属性

  2. 超参数存储(第 8-23 行):将所有 SVM 超参数直接存储为实例属性,供 fit() 和预测方法使用。这种设计让子类只需设置 _impl 即可复用基类的完整参数管理


源码路径:sklearn/svm/_base.py - BaseLibSVM.__sklearn_tags__()(118-123行)

    def __sklearn_tags__(self):
        tags = super().__sklearn_tags__()
        # Used by cross_val_score.
        tags.input_tags.pairwise = self.kernel == "precomputed"  # 接受预计算核矩阵
        tags.input_tags.sparse = self.kernel != "precomputed"    # 非 precomputed 时接受稀疏输入
        return tags

代码解析:

此方法声明估计器的数据能力标签,供 cross_val_score 等工具检查数据兼容性:

  • pairwise=True 表示接受预计算核矩阵(n_samples × n_samples)

  • sparse=True 表示接受稀疏输入(非 precomputed 核时)

这就像物流中心在门口挂牌说明“接受标准货箱”和“接受预打包托盘”两种入库方式。


源码路径:sklearn/svm/_base.py - BaseLibSVM._validate_targets()(252-260行)

    def _validate_targets(self, y):
        """Validation of y and class_weight.

        Default implementation for SVR and one-class; overridden in BaseSVC.
        """
        return column_or_1d(y, warn=True).astype(np.float64, copy=False)

代码解析:

基类默认实现用于 SVR 和 OneClassSVM:将 y 转为一维 float64 数组。BaseSVC 会重写此方法处理分类标签编码和类别权重计算。这体现了模板方法模式——基类定义流程骨架,子类填充特定验证逻辑。


设计中的取舍:BaseLibSVM 的统一调度 vs 子类特化

BaseLibSVM 采用模板方法模式,将通用的 fit/predict/decision_function 流程固化在基类,通过 _validate_targets_sparse 标志、_impl 类型等钩子让子类注入差异化逻辑。这种设计的优势是:

  • 代码复用:核函数计算、稀疏/稠密分发、预测验证等通用逻辑只写一次

  • 一致性保证:所有 SVM 变体共享相同的验证流程、异常处理、后处理逻辑

  • 扩展性:新增 SVM 变体(如加权 SVR)只需设置 _impl 和重写 _validate_targets

劣势是基类膨胀,fit() 方法承担过多职责(验证、分发、参数解析、后处理)。替代方案是使用组合模式将验证、核计算、求解器调用拆分为独立组件,但会增加调用栈深度和数据拷贝开销。scikit-learn 选择继承+模板方法,在保持性能的前提下实现了足够的灵活性。

11.5 Python 层预测引擎:决策函数的“双通道输出”

稠密/稀疏预测的分野

正如物流中心针对标准货箱和散件货物设有不同的分拣线,SVM 预测引擎维护双通道:

  • _dense_fit 调用 libsvm.fit 直接使用内存视图传递数据

  • _sparse_fit 先排序 CSR 索引,再调用 libsvm_sparse.libsvm_sparse_train

  • 稀疏拟合后重建 dual_coef_ 为 CSR 矩阵,支持多分类的块状布局

预测时的模型重建

每次 predict/decision_function 都从当前属性重建 C 层 svm_model 结构,set_model 负责将 NumPy 数组转换为 libsvm 的 svm_node 表示。precomputed 核在预测时通过 support_ 索引直接查表。这就像每次发货前都要根据最新库存重新生成拣货单。

二分类的符号翻转陷阱

libsvm 内部以 +1/-1 编码标签,sklearn 需要翻回用户原始标签顺序。fit() 末尾对二分类的 intercept_dual_coef_ 取反,保证语义一致。这是 libsvm 内部编码与 sklearn 外部语义的“翻译层”。

预测前的安全校验

_validate_for_predict 检查模型是否已拟合,检查 n_support_.sum()support_vectors_.shape[0] 的一致性(防止 CVE-2020-28975),对稀疏输入进行 sort_indices() 保证索引有序。这是发货前的“安检复核”。

源码路径:sklearn/svm/_base.py - BaseLibSVM.predict()(262-280行)

    def predict(self, X):
        """Perform regression on samples in X.

        For a one-class model, +1 (inlier) or -1 (outlier) is returned.

        Parameters
        ----------
        X : {array-like, sparse matrix} of shape (n_samples, n_features)
            For kernel="precomputed", the expected shape of X is
            (n_samples_test, n_samples_train).

        Returns
        -------
        y_pred : ndarray of shape (n_samples,)
            The predicted values.
        """
        X = self._validate_for_predict(X)  # 预测前安全校验(含 is_fitted、稀疏/稠密一致性、precomputed 形状)
        predict = self._sparse_predict if self._sparse else self._dense_predict  # 双通道分发
        return predict(X)

代码解析:

预测统一入口:先调用 _validate_for_predict 安全检查,再根据 _sparse 标志分发到稠密/稀疏预测方法。这种双通道分发机制确保了同一模型能同时处理稠密和稀疏测试数据。


源码路径:sklearn/svm/_base.py - BaseLibSVM._dense_fit()(315-356行)

    def _dense_fit(self, X, y, sample_weight, solver_type, kernel, random_seed):
        if callable(self.kernel):  # callable 核:存储训练数据引用,计算核矩阵替换 X
            # you must store a reference to X to compute the kernel in predict
            # TODO: add keyword copy to copy on demand
            self.__Xfit = X
            X = self._compute_kernel(X)

            if X.shape[0] != X.shape[1]:
                raise ValueError("X.shape[0] should be equal to X.shape[1]")

        libsvm.set_verbosity_wrap(self.verbose)  # 设置 C 层输出详细程度

        # we don't pass **self.get_params() to allow subclasses to
        # add other parameters to __init__
        (
            self.support_,
            self.support_vectors_,
            self._n_support,
            self.dual_coef_,
            self.intercept_,
            self._probA,
            self._probB,
            self.fit_status_,
            self._num_iter,
        ) = libsvm.fit(  # 调用 Cython 绑定的 libsvm.fit
            X,
            y,
            svm_type=solver_type,
            sample_weight=sample_weight,
            class_weight=getattr(self, "class_weight_", np.empty(0)),
            kernel=kernel,
            C=self.C,
            nu=self.nu,
            probability=self._effective_probability,
            degree=self.degree,
            shrinking=self.shrinking,
            tol=self.tol,
            cache_size=self.cache_size,
            coef0=self.coef0,
            gamma=self._gamma,
            epsilon=self.epsilon,
            max_iter=self.max_iter,
            random_seed=random_seed,
        )

        self._warn_from_fit_status()  # 检查收敛状态,必要时发出警告

代码解析:

稠密训练路径核心逻辑:

  1. Callable 核处理(第 2-8 行):存储训练数据引用 __Xfit,计算核矩阵替换 X,验证方阵

  2. Verbosity 设置(第 10 行):传递 verbose 标志控制 C 层输出

  3. 调用 libsvm.fit(第 13-38 行):传递所有超参数,接收支持向量、对偶系数、截距等返回值

  4. 收敛警告(第 40 行):检查 fit_status_ 发出 ConvergenceWarning

注意:callable 核时必须保存训练数据用于预测时计算核矩阵,这就是 __Xfit 的作用。


源码路径:sklearn/svm/_base.py - BaseLibSVM._sparse_fit()(358-407行)

    def _sparse_fit(self, X, y, sample_weight, solver_type, kernel, random_seed):
        X.data = np.asarray(X.data, dtype=np.float64, order="C")  # 强制 float64 data
        X.sort_indices()  # 排序 indices(libsvm 要求)

        kernel_type = self._sparse_kernels.index(kernel)  # 字符串核名转索引

        libsvm_sparse.set_verbosity_wrap(self.verbose)

        (
            self.support_,
            self.support_vectors_,
            dual_coef_data,
            self.intercept_,
            self._n_support,
            self._probA,
            self._probB,
            self.fit_status_,
            self._num_iter,
        ) = libsvm_sparse.libsvm_sparse_train(  # 调用稀疏训练绑定
            X.shape[1],
            X.data,
            X.indices,
            X.indptr,
            y,
            solver_type,
            kernel_type,
            self.degree,
            self._gamma,
            self.coef0,
            self.tol,
            self.C,
            getattr(self, "class_weight_", np.empty(0)),
            sample_weight,
            self.nu,
            self.cache_size,
            self.epsilon,
            int(self.shrinking),
            int(self._effective_probability),
            self.max_iter,
            random_seed,
        )

        self._warn_from_fit_status()

        if hasattr(self, "classes_"):  # 分类器:n_class = n_classes - 1 (OvO)
            n_class = len(self.classes_) - 1
        else:  # regression
            n_class = 1
        n_SV = self.support_vectors_.shape[0]

        dual_coef_indices = np.tile(np.arange(n_SV), n_class)  # 构建 CSR 索引
        if not n_SV:
            self.dual_coef_ = sp.csr_matrix([])
        else:
            dual_coef_indptr = np.arange(
                0, dual_coef_indices.size + 1, dual_coef_indices.size / n_class
            )
            self.dual_coef_ = sp.csr_matrix(  # 重建 CSR 矩阵,形状 (n_class, n_SV)
                (dual_coef_data, dual_coef_indices, dual_coef_indptr), (n_class, n_SV)
            )

代码解析:

稀疏训练路径特有处理:

  1. CSR 标准化(第 2-3 行):强制 float64 data、排序 indices(libsvm 要求)

  2. 核类型映射(第 5 行):字符串核名转 _sparse_kernels 索引

  3. 调用 libsvm_sparse_train(第 7-25 行):传递 CSR 三元组 (data, indices, indptr) 和超参数

  4. dual_coef_ 重建(第 29-42 行):将返回的扁平 dual_coef_data 重构为 CSR 矩阵,形状 (n_class, n_SV),支持多分类 OvO 的块状布局


源码路径:sklearn/svm/_base.py - BaseLibSVM._dense_predict()(282-315行)

    def _dense_predict(self, X):
        X = self._compute_kernel(X)  # callable 核时计算核矩阵
        if X.ndim == 1:
            X = check_array(X, order="C", accept_large_sparse=False)

        kernel = self.kernel
        if callable(self.kernel):
            kernel = "precomputed"
            if X.shape[1] != self.shape_fit_[0]:  # 验证测试样本数匹配训练样本数
                raise ValueError(
                    "X.shape[1] = %d should be equal to %d, "
                    "the number of samples at training time"
                    % (X.shape[1], self.shape_fit_[0])
                )

        svm_type = LIBSVM_IMPL.index(self._impl)

        return libsvm.predict(  # 调用 Cython 绑定预测
            X,
            self.support_,
            self.support_vectors_,
            self._n_support,
            self._dual_coef_,
            self._intercept_,
            self._probA,
            self._probB,
            svm_type=svm_type,
            kernel=kernel,
            degree=self.degree,
            coef0=self.coef0,
            gamma=self._gamma,
            cache_size=self.cache_size,
        )

代码解析:

稠密预测路径:

  1. 核计算(第 2 行):callable 核时计算测试集与训练集的核矩阵

  2. 输入验证(第 3-4 行):1D 输入转为 2D C 连续数组

  3. Callable 核特殊处理(第 6-11 行):验证测试样本数与训练样本数匹配

  4. 调用 libsvm.predict(第 13-26 行):传递模型参数(支持向量、对偶系数、截距等)和核参数


源码路径:sklearn/svm/_base.py - BaseLibSVM._sparse_predict()(316-350行)

    def _sparse_predict(self, X):
        # Precondition: X is a csr_matrix of dtype np.float64.
        kernel = self.kernel
        if callable(kernel):
            kernel = "precomputed"

        kernel_type = self._sparse_kernels.index(kernel)

        C = 0.0  # C is not useful here

        return libsvm_sparse.libsvm_sparse_predict(  # 直接传递 CSR 三元组
            X.data,
            X.indices,
            X.indptr,
            self.support_vectors_.data,
            self.support_vectors_.indices,
            self.support_vectors_.indptr,
            self._dual_coef_.data,
            self._intercept_,
            LIBSVM_IMPL.index(self._impl),
            kernel_type,
            self.degree,
            self._gamma,
            self.coef0,
            self.tol,
            C,
            getattr(self, "class_weight_", np.empty(0)),
            self.nu,
            self.epsilon,
            self.shrinking,
            self._effective_probability,
            self._n_support,
            self._probA,
            self._probB,
        )

代码解析:

稀疏预测路径直接传递 CSR 三元组:

  1. 测试数据 CSR:X.data, X.indices, X.indptr

  2. 支持向量 CSR:support_vectors_.data, .indices, .indptr

  3. 对偶系数数据:_dual_coef_.data(CSR 矩阵的 data 数组)

  4. 其余超参数透传,C 设为 0(预测时不需要)


源码路径:sklearn/svm/_base.py - BaseLibSVM._compute_kernel()(352-360行)

    def _compute_kernel(self, X):
        """Return the data transformed by a callable kernel"""
        if callable(self.kernel):
            # in the case of precomputed kernel given as a function, we
            # have to compute explicitly the kernel matrix
            kernel = self.kernel(X, self.__Xfit)  # 调用用户提供的核函数
            if sp.issparse(kernel):
                kernel = kernel.toarray()
            X = np.asarray(kernel, dtype=np.float64, order="C")
        return X

代码解析:

Callable 核的显式计算:调用用户提供的核函数 kernel(X, X_fit),稀疏结果转密集,确保 float64 C 连续。这是预测时动态计算核矩阵的关键步骤。


源码路径:sklearn/svm/_base.py - BaseLibSVM._dense_decision_function()(420-455行)

    def _dense_decision_function(self, X):
        X = check_array(X, dtype=np.float64, order="C", accept_large_sparse=False)

        kernel = self.kernel
        if callable(kernel):
            kernel = "precomputed"

        return libsvm.decision_function(  # 调用 Cython 绑定计算决策值
            X,
            self.support_,
            self.support_vectors_,
            self._n_support,
            self._dual_coef_,
            self._intercept_,
            self._probA,
            self._probB,
            svm_type=LIBSVM_IMPL.index(self._impl),
            kernel=kernel,
            degree=self.degree,
            cache_size=self.cache_size,
            coef0=self.coef0,
            gamma=self._gamma,
        )

代码解析:

稠密决策函数:验证输入格式,callable 核标记为 precomputed,调用 libsvm.decision_function 计算决策值。注意不需要传递 tol、C、nu 等训练参数,预测只需核参数和模型参数。


源码路径:sklearn/svm/_base.py - BaseLibSVM._sparse_decision_function()(457-496行)

    def _sparse_decision_function(self, X):
        X.data = np.asarray(X.data, dtype=np.float64, order="C")

        kernel = self.kernel
        if hasattr(kernel, "__call__"):
            kernel = "precomputed"

        kernel_type = self._sparse_kernels.index(kernel)

        return libsvm_sparse.libsvm_sparse_decision_function(  # 稀疏决策函数绑定
            X.data,
            X.indices,
            X.indptr,
            self.support_vectors_.data,
            self.support_vectors_.indices,
            self.support_vectors_.indptr,
            self._dual_coef_.data,
            self._intercept_,
            LIBSVM_IMPL.index(self._impl),
            kernel_type,
            self.degree,
            self._gamma,
            self.coef0,
            self.tol,
            self.C,
            getattr(self, "class_weight_", np.empty(0)),
            self.nu,
            self.epsilon,
            self.shrinking,
            self._effective_probability,
            self._n_support,
            self._probA,
            self._probB,
        )

代码解析:

稀疏决策函数:与 _sparse_predict 类似传递 CSR 三元组,但额外传递了 tol、C、nu 等参数(某些变体的决策函数计算需要)。


源码路径:sklearn/svm/_base.py - BaseLibSVM._decision_function()(398-418行)

    def _decision_function(self, X):
        """Evaluates the decision function for the samples in X.

        Parameters
        ----------
        X : array-like of shape (n_samples, n_features)

        Returns
        -------
        X : array-like of shape (n_samples, n_class * (n_class-1) / 2)
            Returns the decision function of the sample for each class
            in the model.
        """
        # NOTE: _validate_for_predict contains check for is_fitted
        # hence must be placed before any other attributes are used.
        X = self._validate_for_predict(X)  # 必须最先调用(含 is_fitted 检查)
        X = self._compute_kernel(X)  # callable 核时计算核矩阵

        if self._sparse:
            dec_func = self._sparse_decision_function(X)  # 双通道分发
        else:
            dec_func = self._dense_decision_function(X)

        # In binary case, we need to flip the sign of coef, intercept and
        # decision function.
        if self._impl in ["c_svc", "nu_svc"] and len(self.classes_) == 2:  # 二分类符号翻转
            return -dec_func.ravel()

        return dec_func

代码解析:

决策函数统一入口,体现双通道分发和二分类符号翻转:

  1. 安全验证(第 9 行):_validate_for_predict 必须最先调用(含 is_fitted 检查)

  2. 核计算(第 10 行):callable 核时计算核矩阵

  3. 双通道分发(第 12-15 行):根据 _sparse 选择稀疏/稠密路径

  4. 二分类符号翻转(第 17-19 行):C-SVC/Nu-SVC 二分类时取反并展平,映射回用户原始标签顺序


源码路径:sklearn/svm/_base.py - BaseLibSVM._validate_for_predict()(498-535行)

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

        if not callable(self.kernel):
            X = validate_data(  # 非 callable 核:验证数据格式
                self,
                X,
                accept_sparse="csr",
                dtype=np.float64,
                order="C",
                accept_large_sparse=False,
                reset=False,
            )

        if self._sparse and not sp.issparse(X):  # 模型稀疏但输入稠密:转 CSR
            X = sp.csr_matrix(X)
        if self._sparse:
            X.sort_indices()  # 确保索引有序

        if sp.issparse(X) and not self._sparse and not callable(self.kernel):  # 模型稠密但输入稀疏:报错
            raise ValueError(
                "cannot use sparse input in %r trained on dense data"
                % type(self).__name__
            )

        if self.kernel == "precomputed":  # precomputed 核形状检查
            if X.shape[1] != self.shape_fit_[0]:
                raise ValueError(
                    "X.shape[1] = %d should be equal to %d, "
                    "the number of samples at training time"
                    % (X.shape[1], self.shape_fit_[0])
                )
        # Fixes https://nvd.nist.gov/vuln/detail/CVE-2020-28975
        # Check that _n_support is consistent with support_vectors
        sv = self.support_vectors_
        if not self._sparse and sv.size > 0 and self.n_support_.sum() != sv.shape[0]:  # 内部一致性检查
            raise ValueError(
                f"The internal representation of {self.__class__.__name__} was altered"
            )
        return X

代码解析:

预测前综合安全校验:

  1. 已拟合检查(第 2 行):check_is_fitted

  2. 非 callable 核的数据验证(第 4-10 行):float64、C 连续、接受 CSR

  3. 稀疏/稠密一致性(第 12-17 行):模型稀疏但输入稠密则转 CSR;模型稠密但输入稀疏报错

  4. Precomputed 核形状检查(第 19-24 行):测试样本数必须等于训练样本数

  5. 内部一致性检查(第 26-30 行):n_support_.sum() == support_vectors_.shape[0] 防止模型被篡改(CVE-2020-28975 修复)


设计中的取舍:双通道预测 vs 统一接口

BaseLibSVM 维护 _dense_*_sparse_* 两套完整的预测路径,通过 _sparse 布尔标志在 Python 层分发。这种设计的优势:

  • 零拷贝直通:稠密用内存视图(float64_t[:, ::1]),稀疏用 CSR 三元组直通,避免格式转换开销

  • 算法特化:稠密路径利用 BLAS dot 加速核计算,稀疏路径利用索引稀疏性跳过零元素

  • 内存效率:稀疏模型存储 CSR 格式的 support_vectors_ 和 dual_coef_,预测时直接传递

劣势是代码重复:预测、决策函数、概率预测各需维护两套实现。替代方案是统一转为稠密或统一转为稀疏,但会牺牲大规模稀疏数据的性能优势。scikit-learn 选择“以空间换时间”,接受代码重复换取最优性能。

11.6 BaseSVC 分类器层:多分类的“合纵连横”

one-vs-one 的内部机制

libsvm 对多分类天然使用 OvO 策略,训练 n*(n-1)/2 个二分类器,就像物流中心为每一对城市建立直达专线。decision_function_shape='ovr' 通过 _ovr_decision_function 从 OvO 投票转换为 OvR 格式,break_ties=True 时用 argmax 打破平局,代价是额外的决策函数计算。

_one_vs_one_coef 的线性核捷径

线性核下可直接从 dual_coef 和 support_vectors 重构原始权重,嵌套循环遍历所有类别对,用 safe_sparse_dot 累加 alpha*x。这避免了显式存储完整的 OvO 系数矩阵,节省内存。

predict_proba 的条件暴露

使用 @available_if(_check_proba) 装饰器动态控制方法可见性。_check_proba 检查 probability 参数和 SVM 类型(仅 SVC/NuSVC 支持)。概率模型通过 5 折交叉验证训练 Platt scaling 参数 probA/probB。

标签编码与类别权重

_validate_targets 使用 np.unique 求出类别并按排序编码,compute_class_weight 支持 'balanced' 和自定义字典,标签编码结果通过 classes_.take() 映射回原始标签。

源码路径:sklearn/svm/_base.py - _one_vs_one_coef()(50-83行)

def _one_vs_one_coef(dual_coef, n_support, support_vectors):
    """Generate primal coefficients from dual coefficients
    for the one-vs-one multi class LibSVM in the case
    of a linear kernel."""

    # get 1vs1 weights for all n*(n-1) classifiers.
    # this is somewhat messy.
    # shape of dual_coef_ is nSV * (n_classes -1)
    # see docs for details
    n_class = dual_coef.shape[0] + 1  # 类别数 = dual_coef 行数 + 1

    # XXX we could do preallocation of coef but
    # would have to take care in the sparse case
    coef = []
    sv_locs = np.cumsum(np.hstack([[0], n_support]))  # 支持向量分段索引累积和
    for class1 in range(n_class):
        # SVs for class1:
        sv1 = support_vectors[sv_locs[class1] : sv_locs[class1 + 1], :]
        for class2 in range(class1 + 1, n_class):
            # SVs for class1:
            sv2 = support_vectors[sv_locs[class2] : sv_locs[class2 + 1], :]

            # dual coef for class1 SVs:
            alpha1 = dual_coef[class2 - 1, sv_locs[class1] : sv_locs[class1 + 1]]
            # dual coef for class2 SVs:
            alpha2 = dual_coef[class1, sv_locs[class2] : sv_locs[class2 + 1]]
            # build weight for class1 vs class2

            coef.append(safe_sparse_dot(alpha1, sv1) + safe_sparse_dot(alpha2, sv2))  # 稀疏安全点积累加
    return coef

代码解析:

线性核下 OvO 对偶系数转原始权重的核心算法:

  1. 类别数推导(第 9 行):dual_coef.shape[0] + 1(OvO 有 n_class-1 行)

  2. 支持向量分段索引(第 13 行):sv_locs 累积和定位各类别 SV 范围

  3. 双层循环遍历类别对(第 14-15 行):class1 < class2,共 n*(n-1)/2 对

  4. 提取各类别 SV 及对应 alpha(第 17-22 行):注意 dual_coef 行索引的偏移逻辑

  5. 稀疏安全点积累加(第 24 行):safe_sparse_dot(alpha1, sv1) + safe_sparse_dot(alpha2, sv2)

  6. 返回权重列表(第 25 行):长度 n*(n-1)/2,每个元素形状 (n_features,)

这种实现避免了显式构建巨大的 OvO 系数矩阵,直接计算每个类别对的原始权重向量。


源码路径:sklearn/svm/_base.py - BaseSVC._validate_targets()(570-600行)

    def _validate_targets(self, y):
        y_ = column_or_1d(y, warn=True)
        check_classification_targets(y)
        cls, y = np.unique(y_, return_inverse=True)  # 排序类别 + 逆映射索引
        self.class_weight_ = compute_class_weight(self.class_weight, classes=cls, y=y_)  # 计算类别权重
        if len(cls) < 2:
            raise ValueError(
                "The number of classes has to be greater than one; got %d class"
                % len(cls)
            )

        self.classes_ = cls  # 保存排序后的类别列表

        return np.asarray(y, dtype=np.float64, order="C")  # 返回编码后标签

代码解析:

分类器专用的目标验证:

  1. 一维化与检查(第 2-3 行):转为 1D,验证是分类目标

  2. 标签编码(第 4 行):np.unique 返回排序后的类别 cls 和逆映射索引 y(0,1,2...)

  3. 类别权重计算(第 5 行):支持 'balanced' 和 dict,基于原始标签 y_ 计算

  4. 类别数检查(第 6-10 行):至少 2 类

  5. 保存类别列表(第 12 行):classes_ 供预测时映射回原始标签

  6. 返回编码后标签(第 14 行):float64 C 连续数组供 libsvm 使用


源码路径:sklearn/svm/_base.py - BaseSVC.decision_function()(602-640行)

    def decision_function(self, X):
        """Evaluate the decision function for the samples in X.

        Parameters
        ----------
        X : array-like of shape (n_samples, n_features)
            The input samples.

        Returns
        -------
        X : ndarray of shape (n_samples, n_classes * (n_classes-1) / 2)
            Returns the decision function of the sample for each class
            in the model.
            If decision_function_shape='ovr', the shape is (n_samples,
            n_classes).

        Notes
        -----
        If decision_function_shape='ovo', the function values are proportional
        to the distance of the samples X to the separating hyperplane. If the
        exact distances are required, divide the function values by the norm of
        the weight vector (``coef_``). See also `this question
        <https://stats.stackexchange.com/questions/14876/
        interpreting-distance-from-hyperplane-in-svm>`_ for further details.
        If decision_function_shape='ovr', the decision function is a monotonic
        transformation of ovo decision function.
        """
        dec = self._decision_function(X)  # 获取原始 OvO 决策值
        if self.decision_function_shape == "ovr" and len(self.classes_) > 2:  # OvR 转换
            return _ovr_decision_function(dec < 0, -dec, len(self.classes_))
        return dec

代码解析:

决策函数形状控制:

  1. 获取原始 OvO 决策值(第 24 行):调用基类 _decision_function,形状 (n_samples, n_class*(n_class-1)/2)

  2. OvR 转换(第 25-26 行):若 decision_function_shape='ovr' 且多类,调用 _ovr_decision_function 将 OvO 投票转为 OvR 置信度

  3. 直接返回 OvO(第 27 行):否则返回原始 OvO 决策值

_ovr_decision_function(dec < 0, -dec, n_classes) 利用决策值符号进行投票统计,再取反得到类似 OvR 的置信度分数。


源码路径:sklearn/svm/_base.py - BaseSVC.predict()(642-680行)

    def predict(self, X):
        """Perform classification on samples in X.

        For a one-class model, +1 or -1 is returned.

        Parameters
        ----------
        X : {array-like, sparse matrix} of shape (n_samples, n_features) or \
                (n_samples_test, n_samples_train)
            For kernel="precomputed", the expected shape of X is
            (n_samples_test, n_samples_train).

        Returns
        -------
        y_pred : ndarray of shape (n_samples,)
            Class labels for samples in X.
        """
        check_is_fitted(self)
        if self.break_ties and self.decision_function_shape == "ovo":
            raise ValueError(
                "break_ties must be False when decision_function_shape is 'ovo'"
            )

        if (
            self.break_ties
            and self.decision_function_shape == "ovr"
            and len(self.classes_) > 2
        ):
            y = np.argmax(self.decision_function(X), axis=1)  # OvR argmax 打破平局
        else:
            y = super().predict(X)  # 默认路径:调用 BaseLibSVM.predict(OvO 投票)
        return self.classes_.take(np.asarray(y, dtype=np.intp))  # 映射回原始标签

代码解析:

分类预测入口,处理 break_ties 逻辑:

  1. 拟合检查(第 13 行)

  2. 参数冲突检查(第 14-17 行):OvO 模式下不可用 break_ties

  3. break_ties=True 时的 OvR 投票(第 19-23 行):多类 OvR 时用 argmax 打破平局

  4. 默认路径(第 24 行):调用父类 predict()(即 BaseLibSVM.predict)得到 OvO 投票结果

  5. 标签映射(第 25 行):classes_.take() 将数值索引映射回原始标签


源码路径:sklearn/svm/_base.py - BaseSVC._check_proba()(682-695行)

    def _check_proba(self):
        if self.probability == "deprecated" or not self.probability:
            raise AttributeError(
                "predict_proba is not available when probability=False"
            )
        if self._impl not in ("c_svc", "nu_svc"):
            raise AttributeError("predict_proba only implemented for SVC and NuSVC")
        return True

代码解析:

概率预测可用性检查,供 @available_if 装饰器使用:

  1. probability 参数检查"deprecated" 字符串或 False 时不可用

  2. SVM 类型检查:仅 C-SVC 和 Nu-SVC 支持概率输出

  3. 返回 True:表示方法可用,否则抛出 AttributeError 导致方法不可见


源码路径:sklearn/svm/_base.py - BaseSVC.predict_proba()(697-735行)

    @available_if(_check_proba)
    def predict_proba(self, X):
        """Compute probabilities of possible outcomes for samples in X.

        The model needs to have probability information computed at training
        time: fit with attribute `probability` set to True.

        Parameters
        ----------
        X : array-like of shape (n_samples, n_features)
            For kernel="precomputed", the expected shape of X is
            (n_samples_test, n_samples_train).

        Returns
        -------
        T : ndarray of shape (n_samples, n_classes)
            Returns the probability of the sample for each class in
            the model. The columns correspond to the classes in sorted
            order, as they appear in the attribute :term:`classes_`.

        Notes
        -----
        The probability model is created using cross validation, so
        the results can be slightly different than those obtained by
        predict. Also, it will produce meaningless results on very small
        datasets.
        """
        X = self._validate_for_predict(X)
        if self.probA_.size == 0 or self.probB_.size == 0:
            raise NotFittedError(
                "predict_proba is not available when fitted with probability=False"
            )
        pred_proba = (
            self._sparse_predict_proba if self._sparse else self._dense_predict_proba
        )
        return pred_proba(X)

代码解析:

概率预测入口:

  1. 动态可见性@available_if(_check_proba) 仅在概率可用时暴露方法

  2. 安全验证(第 13 行):_validate_for_predict

  3. 概率模型存在性检查(第 14-18 行):probA_/probB_ 非空验证

  4. 双通道分发(第 19-21 行):根据 _sparse 选择稀疏/稠密概率预测


设计中的取舍:OvO 训练 vs OvR 推理

LibSVM 原始实现默认使用一对一(OvO)策略进行多分类训练,训练 n*(n-1)/2 个二分类器。scikit-learn 保留了这一默认行为,但通过 decision_function_shape='ovr' 提供了一个后处理步骤,将 OvO 的决策值转换为 OvR 格式。

  • 优势:OvO 训练在二分类问题上计算量较小,且每个子问题规模较小,有利于并行化和数值稳定性;决策函数的后处理(如 _ovr_decision_function)仅在预测时进行,开销可控。

  • 劣势:相比直接训练 OvR 模型,OvO 需要存储和计算更多的中间模型,在内存和预测延迟上略有不利。

  • 决策:scikit-learn 选择保留 LibSVM 的原始 OvO 训练流程,因为它已经经过高度优化(工作集选择、核缓存、收缩启发式),而将 OvR 转换作为轻量级的后处理步骤,避免重新实现整个训练算法,同时为用户提供熟悉的 OvR 接口。

11.7 libsvm 核心:SMO 求解器与核缓存

11.7.1 Cache 类:LRU 核缓存机制

classDiagram class Cache { +int l +long int size +head_t *head +head_t lru_head +Cache(int l, long int size) +~Cache() +int get_data(int index, Qfloat** data, int len) +void swap_index(int i, int j) -void lru_delete(head_t* h) -void lru_insert(head_t* h) } class head_t { +head_t *prev +head_t *next +Qfloat *data +int len } Cache --> head_t : 包含

Cache 类实现 LRU(最近最少使用)缓存策略,管理核矩阵列的缓存。每个样本对应一个 head_t 节点,缓存该样本与所有样本的核函数值(即 Q 矩阵的一列)。LRU 双向链表维护访问顺序,最近使用的移到尾部,淘汰时从头部移除。

源码路径:sklearn/svm/src/libsvm/svm.cpp - Cache 类(约 100-200 行)

class Cache
{
public:
    Cache(int l_, long int size_) : l(l_), size(size_)
    {
        head = (head_t *)calloc(l, sizeof(head_t));  // 初始化为 0
        size /= sizeof(Qfloat);                      // 换算为 Qfloat 元素数量
        size -= l * sizeof(head_t) / sizeof(Qfloat); // 减去 head_t 数组占用的空间
        size = max(size, 2 * (long int) l);          // 缓存至少要能容纳 2 列
        lru_head.next = lru_head.prev = &lru_head;   // 空的循环双向链表
    }

    ~Cache()
    {
        for (head_t *h = lru_head.next; h != &lru_head; h = h->next)
            free(h->data);
        free(head);
    }

    // 请求第 index 列的前 len 个元素
    // 返回需要填充的起始位置 p(p >= len 表示全命中)
    int get_data(const int index, Qfloat **data, int len)
    {
        head_t *h = &head[index];
        if (h->len) lru_delete(h);  // 若已在链表中,先删除
        int more = len - h->len;    // 需要新增的元素数

        if (more > 0)
        {
            // 缓存空间不足,淘汰 LRU 节点
            while (size < more)
            {
                head_t *old = lru_head.next;  // 头部为最久未使用
                lru_delete(old);
                free(old->data);
                size += old->len;
                old->data = 0;
                old->len = 0;
            }

            // 分配新空间
            h->data = (Qfloat *)realloc(h->data, sizeof(Qfloat) * len);
            size -= more;
            swap(h->len, len);
        }

        lru_insert(h);  // 插入到尾部(最近使用)
        *data = h->data;
        return len;
    }

    // 索引交换时同步缓存内容
    void swap_index(int i, int j)
    {
        if (i == j) return;

        if (head[i].len) lru_delete(&head[i]);
        if (head[j].len) lru_delete(&head[j]);
        swap(head[i].data, head[j].data);
        swap(head[i].len, head[j].len);
        if (head[i].len) lru_insert(&head[i]);
        if (head[j].len) lru_insert(&head[j]);

        if (i > j) swap(i, j);
        for (head_t *h = lru_head.next; h != &lru_head; h = h->next)
        {
            if (h->len > i)
            {
                if (h->len > j)
                    swap(h->data[i], h->data[j]);  // 交换缓存中的两列数据
                else
                {
                    // 列 j 超出该缓存条目长度,只能失效整条
                    lru_delete(h);
                    free(h->data);
                    size += h->len;
                    h->data = 0;
                    h->len = 0;
                }
            }
        }
    }

private:
    void lru_delete(head_t *h)  // 从 LRU 链表删除
    {
        h->prev->next = h->next;
        h->next->prev = h->prev;
    }

    void lru_insert(head_t *h)  // 插入到 LRU 链表尾部
    {
        h->next = &lru_head;
        h->prev = lru_head.prev;
        h->prev->next = h;
        h->next->prev = h;
    }
};

代码解析:

  1. 构造函数(第 4-11 行):按字节换算缓存容量,减去 head 数组开销,保证至少能存 2 列(避免频繁淘汰)

  2. get_data()(第 19-48 行):请求列数据,缓存命中则直接返回;未命中则按 LRU 淘汰旧数据,分配新内存,更新链表

  3. swap_index()(第 50-75 行):SMO 中变量交换时调用,同步交换缓存内容,失效长度不足的缓存条目

  4. LRU 维护lru_delete/lru_insert 维护双向循环链表,头部为最久未用,尾部为最近使用

11.7.2 Kernel 类:五种核函数实现

Kernel 类封装核函数计算,继承自 QMatrix,提供 get_Q() 按需计算 Q 矩阵列(配合 Cache),k_function() 单次核计算(预测用)。

源码路径:sklearn/svm/src/libsvm/svm.cpp - Kernel 类及核函数实现(约 200-400 行)

class Kernel : public QMatrix {
public:
#ifdef _DENSE_REP
    Kernel(int l, PREFIX(node) * x, const svm_parameter& param, BlasFunctions *blas_functions);
#else
    Kernel(int l, PREFIX(node) * const * x, const svm_parameter& param, BlasFunctions *blas_functions);
#endif
    virtual ~Kernel();

    static double k_function(const PREFIX(node) *x, const PREFIX(node) *y,
                 const svm_parameter& param, BlasFunctions *blas_functions);
    virtual Qfloat *get_Q(int column, int len) const = 0;
    virtual double *get_QD() const = 0;
    virtual void swap_index(int i, int j) const
    {
        swap(x[i], x[j]);
        if (x_square) swap(x_square[i], x_square[j]);
    }
protected:
    double (Kernel::*kernel_function)(int i, int j) const;  // 成员函数指针

private:
#ifdef _DENSE_REP
    PREFIX(node) *x;  // 稠密:连续内存
#else
    const PREFIX(node) **x;  // 稀疏:指针数组
#endif
    double *x_square;  // RBF 核用:预计算 ||x||^2
    BlasFunctions *m_blas;  // BLAS 函数指针

    const int kernel_type;
    const int degree;
    const double gamma;
    const double coef0;

    static double dot(const PREFIX(node) *px, const PREFIX(node) *py, BlasFunctions *blas_functions);
    // ... 五种核实现:
    double kernel_linear(int i, int j) const { return dot(x[i], x[j], m_blas); }
    double kernel_poly(int i, int j) const { return powi(gamma*dot(x[i],x[j],m_blas)+coef0, degree); }
    double kernel_rbf(int i, int j) const { return exp(-gamma*(x_square[i]+x_square[j]-2*dot(x[i],x[j],m_blas))); }
    double kernel_sigmoid(int i, int j) const { return tanh(gamma*dot(x[i],x[j],m_blas)+coef0); }
    double kernel_precomputed(int i, int j) const { /* 直接查表 */ }
};

代码解析:

  1. 构造函数(第 10-35 行):根据 kernel_type 设置 kernel_function 函数指针,RBF 核预计算 x_square

  2. dot() 静态方法(第 40-55 行):稠密用 BLAS dot,稀疏用双指针遍历索引

  3. 五种核实现(第 57-70 行):直接调用 dot 组合,RBF 利用预计算范数避免重复计算

  4. k_function() 静态方法(第 72-110 行):预测时单次核计算,稠密/稀疏分支,RBF 手动展开计算差值平方和

11.7.3 Solver 类:SMO 主循环与工作集选择

SMO(Sequential Minimal Optimization)每次选取两个变量求解二次规划子问题,通过工作集选择启发式加速收敛。

源码路径:sklearn/svm/src/libsvm/svm.cpp - Solver::Solve()select_working_set()(约 500-700 行)

void Solver::Solve(int l, const QMatrix& Q, const double *p_, const schar *y_,
           double *alpha_, const double *C_, double eps,
           SolutionInfo* si, int shrinking, int max_iter)
{
    // ... 初始化 alpha_status, active_set, G, G_bar ...

    int iter = 0;
    int counter = min(l, 1000) + 1;

    while (1)
    {
        if ((max_iter != -1) && (iter >= max_iter)) {  // 迭代上限检查
            si->solve_timed_out = true;
            break;
        }

        // 显示进度 & 收缩
        if (--counter == 0)
        {
            counter = min(l, 1000);
            if (shrinking) do_shrinking();  // 定期执行收缩
            info(".");
        }

        int i, j;
        if (select_working_set(i, j) != 0)  // 选择工作集
        {
            reconstruct_gradient();  // 重建梯度
            active_size = l;
            info("*");
            if (select_working_set(i, j) != 0)
                break;  // 最优解
            else
                counter = 1;
        }

        ++iter;

        // 更新 alpha[i], alpha[j] - 核心 SMO 步骤
        const Qfloat *Q_i = Q.get_Q(i, active_size);
        const Qfloat *Q_j = Q.get_Q(j, active_size);
        double C_i = get_C(i), C_j = get_C(j);
        double old_alpha_i = alpha[i], old_alpha_j = alpha[j];

        if (y[i] != y[j])  // 标签不同
        {
            double quad_coef = QD[i] + QD[j] + 2 * Q_i[j];
            if (quad_coef <= 0) quad_coef = TAU;
            double delta = (-G[i] - G[j]) / quad_coef;
            double diff = alpha[i] - alpha[j];
            alpha[i] += delta; alpha[j] += delta;
            // ... 边界裁剪逻辑 ...
        }
        else  // 标签相同
        {
            double quad_coef = QD[i] + QD[j] - 2 * Q_i[j];
            if (quad_coef <= 0) quad_coef = TAU;
            double delta = (G[i] - G[j]) / quad_coef;
            double sum = alpha[i] + alpha[j];
            alpha[i] -= delta; alpha[j] += delta;
            // ... 边界裁剪逻辑 ...
        }

        // 更新梯度 G
        double delta_alpha_i = alpha[i] - old_alpha_i;
        double delta_alpha_j = alpha[j] - old_alpha_j;
        for (int k = 0; k < active_size; k++)
            G[k] += Q_i[k] * delta_alpha_i + Q_j[k] * delta_alpha_j;

        // 更新 alpha_status 和 G_bar
        // ... 上界/下界变化时更新 G_bar ...
    }
    // ... 计算 rho, 目标值, 保存结果 ...
}

// 工作集选择:两阶段启发式
int Solver::select_working_set(int &out_i, int &out_j)
{
    // 第一阶段:找违反 KKT 最大的 i
    double Gmax = -INF;
    int Gmax_idx = -1;
    for (int t = 0; t < active_size; t++)
        if (y[t] == +1) {
            if (!is_upper_bound(t) && -G[t] >= Gmax) { Gmax = -G[t]; Gmax_idx = t; }
        } else {
            if (!is_lower_bound(t) && G[t] >= Gmax) { Gmax = G[t]; Gmax_idx = t; }
        }

    int i = Gmax_idx;
    const Qfloat *Q_i = (i != -1) ? Q->get_Q(i, active_size) : NULL;

    // 第二阶段:找目标下降最快的 j
    double Gmax2 = -INF;
    int Gmin_idx = -1;
    double obj_diff_min = INF;

    for (int j = 0; j < active_size; j++)
    {
        if (y[j] == +1) {
            if (!is_lower_bound(j)) {
                double grad_diff = Gmax + G[j];
                if (G[j] >= Gmax2) Gmax2 = G[j];
                if (grad_diff > 0) {
                    double quad_coef = QD[i] + QD[j] - 2.0 * y[i] * Q_i[j];
                    double obj_diff = (quad_coef > 0) ? -(grad_diff*grad_diff)/quad_coef : -(grad_diff*grad_diff)/TAU;
                    if (obj_diff <= obj_diff_min) { Gmin_idx = j; obj_diff_min = obj_diff; }
                }
            }
        } else {
            if (!is_upper_bound(j)) {
                double grad_diff = Gmax - G[j];
                if (-G[j] >= Gmax2) Gmax2 = -G[j];
                if (grad_diff > 0) {
                    double quad_coef = QD[i] + QD[j] + 2.0 * y[i] * Q_i[j];
                    double obj_diff = (quad_coef > 0) ? -(grad_diff*grad_diff)/quad_coef : -(grad_diff*grad_diff)/TAU;
                    if (obj_diff <= obj_diff_min) { Gmin_idx = j; obj_diff_min = obj_diff; }
                }
            }
        }
    }

    if (Gmax + Gmax2 < eps || Gmin_idx == -1) return 1;  // 收敛
    out_i = Gmax_idx; out_j = Gmin_idx;
    return 0;
}

代码解析:

  1. 主循环Solve):每次迭代选工作集、更新两个 alpha、更新梯度、定期收缩、检查收敛

  2. 工作集选择select_working_set):

    • 第一阶段:遍历活跃集,找 -y_i * G_i 最大的 i(最违反 KKT 条件)

    • 第二阶段:固定 i,遍历候选 j,计算目标函数下降量 obj_diff,选下降最快的 j

    • 二次项系数 quad_coefQD[i]+QD[j]±2*y[i]*Q_i[j],正定性保证用 TAU 下界

  3. 收敛判断Gmax + Gmax2 < eps 即 KKT 违反程度之和小于容差

11.7.4 收缩启发式:加速训练的关键

flowchart TD A[开始 do_shrinking] --> B[计算 Gmax1, Gmax2] B --> C{unshrink == false 且 Gmax1+Gmax2 <= eps*10?} C -->|是| D[unshrink = true, 重建梯度, active_size = l] C -->|否| E[遍历活跃集] E --> F{be_shrunk(i, Gmax1, Gmax2)?} F -->|是| G[active_size--, 与末尾交换] F -->|否| H[继续下一个] G --> H H --> I{遍历完成?} I -->|否| E I -->|是| J[结束]

收缩启发式将已收敛的变量(alpha 在边界且梯度满足条件)从活跃集移出,只优化活跃变量。

源码路径:sklearn/svm/src/libsvm/svm.cpp - Solver::do_shrinking()be_shrunk()(约 700-800 行)

bool Solver::be_shrunk(int i, double Gmax1, double Gmax2)
{
    if (is_upper_bound(i)) {
        if (y[i] == +1) return (-G[i] > Gmax1);
        else return (-G[i] > Gmax2);
    }
    else if (is_lower_bound(i)) {
        if (y[i] == +1) return (G[i] > Gmax2);
        else return (G[i] > Gmax1);
    }
    else return false;  // FREE 变量不收缩
}

void Solver::do_shrinking()
{
    int i;
    double Gmax1 = -INF, Gmax2 = -INF;

    // 找最大违反对
    for (i = 0; i < active_size; i++) {
        if (y[i] == +1) {
            if (!is_upper_bound(i) && -G[i] >= Gmax1) Gmax1 = -G[i];
            if (!is_lower_bound(i) && G[i] >= Gmax2) Gmax2 = G[i];
        } else {
            if (!is_upper_bound(i) && -G[i] >= Gmax2) Gmax2 = -G[i];
            if (!is_lower_bound(i) && G[i] >= Gmax1) Gmax1 = G[i];
        }
    }

    // 解除收缩条件
    if (unshrink == false && Gmax1 + Gmax2 <= eps * 10) {
        unshrink = true;
        reconstruct_gradient();
        active_size = l;
        info("*");
    }

    // 执行收缩
    for (i = 0; i < active_size; i++)
        if (be_shrunk(i, Gmax1, Gmax2)) {
            active_size--;
            while (active_size > i) {
                if (!be_shrunk(active_size, Gmax1, Gmax2)) {
                    swap_index(i, active_size);
                    break;
                }
                active_size--;
            }
        }
}

代码解析:

  1. be_shrunk():判断变量是否可收缩——上界变量看 -G,下界变量看 G,与标签配合

  2. do_shrinking()

    • 先计算活跃集内的 Gmax1, Gmax2(最大违反梯度)

    • 解除收缩:若 Gmax1+Gmax2 <= eps*10 且未解除过,重建完整梯度,恢复 active_size=l

    • 执行收缩:遍历活跃集,可收缩的与末尾交换并减小 active_size


设计中的取舍:SMO 工作集选择的两阶段启发式

Fan et al. (2005) 提出的两阶段工作集选择:

  • 第一阶段:选违反 KKT 最大的 imax -y_i*G_i

  • 第二阶段:固定 i,选目标下降最快的 jmin obj_diff

替代方案:

  • 全局最优对:遍历所有 i,j 对计算 obj_diff,复杂度 O(n²) 太高

  • 随机选择:收敛慢,不稳定

两阶段启发式在 O(n) 复杂度下近似全局最优,实践中效果极好。收缩启发式进一步将有效 n 缩小到活跃集大小,大幅加速大规模问题。

11.8 libsvm 训练流程:从问题构建到模型输出

flowchart TD A[train() 入口] --> B[remove_zero_weight 移除零权重样本] B --> C{svm_type?} C -->|C_SVC| D[solve_c_svc] C -->|NU_SVC| E[solve_nu_svc] C -->|ONE_CLASS| F[solve_one_class] C -->|EPSILON_SVR| G[solve_epsilon_svr] C -->|NU_SVR| H[solve_nu_svr] D --> I[svm_train_one] E --> I F --> I G --> I H --> I I --> J[Solver/Solver_NU.Solve] J --> K[构建 QMatrix: SVC_Q/ONE_CLASS_Q/SVR_Q] K --> L[SMO 主循环] L --> M[计算 rho, 目标值] M --> N[统计 nSV, nBSV] N --> O[构建 decision_function] O --> P[多分类: svm_group_classes 分组] P --> Q[训练 k*(k-1)/2 个二分类器] Q --> R[组装 model: SV, sv_coef, rho, probA/probB] R --> S[返回 model]

11.8.1 各问题求解入口

源码路径:sklearn/svm/src/libsvm/svm.cpp - solve_c_svc() 等(约 2000-2500 行)

static void solve_c_svc(
    const PREFIX(problem) *prob, const svm_parameter* param,
    double *alpha, Solver::SolutionInfo* si, double Cp, double Cn, BlasFunctions *blas_functions)
{
    int l = prob->l;
    double *minus_ones = new double[l];
    schar *y = new schar[l];
    double *C = new double[l];

    for (int i = 0; i < l; i++) {
        alpha[i] = 0;
        minus_ones[i] = -1;
        if (prob->y[i] > 0) {
            y[i] = +1;
            C[i] = prob->W[i] * Cp;  // 正类权重 * C
        } else {
            y[i] = -1;
            C[i] = prob->W[i] * Cn;  // 负类权重 * C
        }
    }

    Solver s;
    s.Solve(l, SVC_Q(*prob, *param, y, blas_functions), minus_ones, y,
        alpha, C, param->eps, si, param->shrinking, param->max_iter);

    for (int i = 0; i < l; i++)
        alpha[i] *= y[i];  // 恢复 alpha * y

    delete[] C; delete[] minus_ones; delete[] y;
}

代码解析:

  • 标签转换y > 0+1y <= 0-1

  • 类别权重C[i] = W[i] * Cp/Cn,样本权重 W[i] 乘以类别惩罚参数

  • Q 矩阵SVC_Q 继承 Kernelget_Q 返回 y_i*y_j*K(x_i,x_j)

  • 目标函数min 0.5*alpha^T Q alpha - sum(alpha),约束 y^T alpha = 0, 0 <= alpha_i <= C_i

11.8.2 多分类训练:OvO 组装

源码路径:sklearn/svm/src/libsvm/svm.cpp - PREFIX(train) 分类分支(约 2600-3000 行)

// 分类分支:训练 k*(k-1)/2 个二分类器
bool *nonzero = Malloc(bool, l);
for (i = 0; i < l; i++) nonzero[i] = false;

decision_function *f = Malloc(decision_function, nr_class*(nr_class-1)/2);
double *probA = NULL, *probB = NULL;
if (param->probability) {
    probA = Malloc(double, nr_class*(nr_class-1)/2);
    probB = Malloc(double, nr_class*(nr_class-1)/2);
}

int p = 0;
for (i = 0; i < nr_class; i++)
    for (int j = i+1; j < nr_class; j++) {
        // 构建子问题:类别 i vs j
        problem sub_prob;
        sub_prob.l = ci + cj;
        sub_prob.y = [+1...+1, -1...-1];
        sub_prob.W = 对应权重;

        if (param->probability)
            svm_binary_svc_probability(&sub_prob, param, weighted_C[i], weighted_C[j],
                probA[p], probB[p], status, blas_functions);

        f[p] = svm_train_one(&sub_prob, param, weighted_C[i], weighted_C[j], status, blas_functions);

        // 标记非零 alpha 对应的样本
        for (k = 0; k < ci; k++) if (!nonzero[si+k] && fabs(f[p].alpha[k]) > 0) nonzero[si+k] = true;
        for (k = 0; k < cj; k++) if (!nonzero[sj+k] && fabs(f[p].alpha[ci+k]) > 0) nonzero[sj+k] = true;
        ++p;
    }

// 组装 model->sv_coef: 交织存储各分类器的 alpha
for (i = 0; i < nr_class; i++)
    for (int j = i+1; j < nr_class; j++) {
        // classifier (i,j): i 的系数在 sv_coef[j-1], j 的系数在 sv_coef[i]
        int q = nz_start[i];
        for (k = 0; k < ci; k++) if (nonzero[si+k]) model->sv_coef[j-1][q++] = f[p].alpha[k];
        q = nz_start[j];
        for (k = 0; k < cj; k++) if (nonzero[sj+k]) model->sv_coef[i][q++] = f[p].alpha[ci+k];
    }

代码解析:

  1. 类别分组svm_group_classes 按标签排序分组,生成 start[], count[], perm[]

  2. 子问题构建:每对类别提取对应样本,标签设为 +1/-1

  3. 概率估计svm_binary_svc_probability 5折 CV 训练 Platt scaling 参数

  4. 支持向量合并nonzero[] 标记任一分类器中 alpha 非零的样本

  5. 系数交织存储sv_coef[j-1] 存类别 i 对 j 的系数,sv_coef[i] 存类别 j 对 i 的系数,按 nz_start 偏移排列


11.9 liblinear 核心:坐标下降与 TRON 信任域

Liblinear 专为线性核设计,无需核缓存,通过坐标下降或 TRON 求解原始/对偶问题。

11.9.1 坐标下降求解器族

classDiagram class SolverBase { <<abstract>> +Solve() } class Solver_CD_Dual { +solve_l2r_l1l2_svc() +solve_l2r_l1l2_svr() +solve_l2r_lr_dual() } class Solver_CD_Primal { +solve_l1r_l2_svc() +solve_l1r_lr() } class Solver_MCSVM_CS { +Solve() +solve_sub_problem() +be_shrunk() } SolverBase <|-- Solver_CD_Dual SolverBase <|-- Solver_CD_Primal SolverBase <|-- Solver_MCSVM_CS

11.9.1.1 对偶坐标下降:L2R_L2LOSS_SVC_DUAL 等

源码路径:sklearn/svm/src/liblinear/linear.cpp - solve_l2r_l1l2_svc()(约 1000-1300 行)

static int solve_l2r_l1l2_svc(
    const problem *prob, double *w, double eps,
    double Cp, double Cn, int solver_type, int max_iter)
{
    // ... 初始化 alpha=0, w=0, QD, index, diag, upper_bound ...
    // L1-loss: diag=0, upper_bound=C
    // L2-loss: diag=1/(2C), upper_bound=INF

    while (iter < max_iter) {
        PGmax_new = -INF; PGmin_new = INF;
        // 随机打乱 index
        for (i = 0; i < active_size; i++) {
            int j = i + bounded_rand_int(active_size - i);
            swap(index[i], index[j]);
        }

        for (s = 0; s < active_size; s++) {
            i = index[s];
            G = y[i] * (w^T x_i) - 1 + alpha[i] * diag[i];  // 梯度

            C = upper_bound[i];
            PG = 计算投影梯度;
            PGmax_new = max(PGmax_new, PG);
            PGmin_new = min(PGmin_new, PG);

            if (fabs(PG) > 1e-12) {
                alpha_old = alpha[i];
                alpha[i] = min(max(alpha[i] - G/QD[i], 0.0), C);  // 投影更新
                d = (alpha[i] - alpha_old) * y[i];
                w += d * x_i;  // 更新原始权重
            }
        }

        iter++;
        if (PGmax_new - PGmin_new <= eps) {
            if (active_size == l) break;
            else { active_size = l; info("*"); continue; }  // 重新激活
        }
        PGmax_old = PGmax_new; PGmin_old = PGmin_new;
    }
    return iter;
}

代码解析:

  1. 问题设定:对偶问题 min 0.5*alpha^T(Q+D)alpha - e^T alpha, s.t. 0<=alpha<=upper_bound

  2. 坐标更新:单变量二次函数最小化 alpha_new = clip(alpha - G/QD, 0, C)

  3. 原始权重同步w += d * x_i 增量更新,利用稀疏性仅遍历非零特征

  4. 收缩:投影梯度 PG 指示最优性,连续多轮 PGmax-PGmin <= eps 则收敛

11.9.1.2 原始坐标下降:L1R_L2LOSS_SVC

源码路径:sklearn/svm/src/liblinear/linear.cpp - solve_l1r_l2_svc()(约 1500-1800 行)

static int solve_l1r_l2_svc(
    problem *prob_col, double *w, double eps,
    double Cp, double Cn, int max_iter)
{
    // 转置为列格式:prob_col->x[j] 是第 j 个特征的所有样本值
    // b = 1 - y*w^T x (margin)

    while (iter < max_iter) {
        Gmax_new = 0; Gnorm1_new = 0;
        // 随机打乱特征索引
        for (j = 0; j < active_size; j++) { swap... }

        for (s = 0; s < active_size; s++) {
            j = index[s];
            G_loss = 0; H = 0;
            // 遍历第 j 个特征非零的样本
            for (x = prob_col->x[j]; x->index != -1; x++) {
                ind = x->index - 1;
                if (b[ind] > 0) {  // hinge loss 活跃
                    G_loss -= C[ind] * x->value * b[ind];
                    H += C[ind] * x->value * x->value;
                }
            }
            G_loss *= 2; G = G_loss; H = max(2*H, 1e-12);

            // 软阈值更新:w_new = sign(-G/H) * max(|G/H| - 1/H, 0)
            if (Gp < H*w[j]) d = -Gp/H;
            else if (Gn > H*w[j]) d = -Gn/H;
            else d = -w[j];

            // 回溯线搜索保证目标下降
            for (num_linesearch = 0; num_linesearch < max_num_linesearch; num_linesearch++) {
                cond = ...;  // 近似条件
                if (appxcond <= 0) break;
                d *= 0.5;
            }
            w[j] += d;
            // 更新 b
        }
        // 收敛判断 & 收缩
    }
}

代码解析:

  1. 列格式转置transpose() 将行格式转为列格式,便于按特征访问

  2. 软阈值操作:L1 正则导致坐标更新带软阈值,实现稀疏解

  3. 回溯线搜索:确保目标函数单调下降,处理非光滑 L1 项

11.9.2 TRON 信任域 Newton 法

用于 L2 正则化的原始问题(L2R_LR, L2R_L2LOSS_SVC, L2R_L2LOSS_SVR)。

flowchart TD A[TRON 主循环] --> B[计算梯度 g, 函数值 f] B --> C[trcg 求解信任域子问题] C --> D[得到步长 s] D --> E[计算 actred/prered] E --> F{actred/prered 比率} F -->|< eta0| G[缩小 delta] F -->|< eta1| H[微调 delta] F -->|< eta2| I[保持/扩大 delta] F -->|>= eta2| J[扩大 delta] G --> K[拒绝步长] H --> K I --> L[接受步长: w = w_new] J --> L K --> B L --> M[更新 g, f, gnorm] M --> N{gnorm <= eps*gnorm1?} N -->|是| O[收敛] N -->|否| B

源码路径:sklearn/svm/src/liblinear/tron.cpp - TRON::tron()trcg()(约 50-150 行)

int TRON::tron(double *w)
{
    double eta0 = 1e-4, eta1 = 0.25, eta2 = 0.75;
    double sigma1 = 0.25, sigma2 = 0.5, sigma3 = 4;

    int n = fun_obj->get_nr_variable();
    double *s = new double[n], *r = new double[n], *w_new = new double[n], *g = new double[n];

    for (i = 0; i < n; i++) w[i] = 0;
    f = fun_obj->fun(w);
    fun_obj->grad(w, g);
    delta = blas->nrm2(n, g, inc);  // 初始信任域半径 = 梯度范数
    double gnorm1 = delta;

    iter = 1;
    while (iter <= max_iter && search) {
        cg_iter = trcg(delta, g, s, r);  // 截断共轭梯度求解子问题

        memcpy(w_new, w, sizeof(double)*n);
        blas->axpy(n, 1.0, s, inc, w_new, inc);  // w_new = w + s

        gs = blas->dot(n, g, inc, s, inc);
        prered = -0.5*(gs - blas->dot(n, s, inc, r, inc));  // 预测下降
        fnew = fun_obj->fun(w_new);
        actred = f - fnew;  // 实际下降

        snorm = blas->nrm2(n, s, inc);
        if (iter == 1) delta = min(delta, snorm);  // 首步调整

        // 更新信任域半径
        if (fnew - f - gs <= 0) alpha = sigma3;
        else alpha = max(sigma1, -0.5*(gs/(fnew - f - gs)));

        if (actred < eta0*prered) delta = min(max(alpha, sigma1)*snorm, sigma2*delta);
        else if (actred < eta1*prered) delta = max(sigma1*delta, min(alpha*snorm, sigma2*delta));
        else if (actred < eta2*prered) delta = max(sigma1*delta, min(alpha*snorm, sigma3*delta));
        else delta = max(delta, min(alpha*snorm, sigma3*delta));

        info("iter %2d act %5.3e pre %5.3e delta %5.3e f %5.3e |g| %5.3e CG %3d\n", ...);

        if (actred > eta0*prered) {  // 接受步长
            iter++;
            memcpy(w, w_new, sizeof(double)*n);
            f = fnew;
            fun_obj->grad(w, g);
            gnorm = blas->nrm2(n, g, inc);
            if (gnorm <= eps*gnorm1) break;
        }
    }
    return --iter;
}

int TRON::trcg(double delta, double *g, double *s, double *r)
{
    // 共轭梯度法求解 min g^T s + 0.5 s^T H s, s.t. ||s|| <= delta
    for (i = 0; i < n; i++) { s[i]=0; r[i]=-g[i]; d[i]=r[i]; }
    cgtol = 0.1 * blas->nrm2(n, g, inc);

    while (1) {
        if (blas->nrm2(n, r, inc) <= cgtol) break;
        cg_iter++;
        fun_obj->Hv(d, Hd);  // Hessian-vector 积

        alpha = rTr / blas->dot(n, d, inc, Hd, inc);
        blas->axpy(n, alpha, d, inc, s, inc);  // s += alpha*d

        if (blas->nrm2(n, s, inc) > delta) {  // 超出信任域边界
            // 截断到边界
            blas->axpy(n, -alpha, d, inc, s, inc);
            // 计算截断 alpha
            blas->axpy(n, alpha, d, inc, s, inc);
            blas->axpy(n, -alpha, Hd, inc, r, inc);
            break;
        }
        blas->axpy(n, -alpha, Hd, inc, r, inc);
        rnewTrnew = blas->dot(n, r, inc, r, inc);
        beta = rnewTrnew / rTr;
        blas->scal(n, beta, d, inc);
        blas->axpy(n, 1.0, r, inc, d, inc);
        rTr = rnewTrnew;
    }
    return cg_iter;
}

代码解析:

  1. 信任域框架:用二次模型 m(s) = f + g^T s + 0.5 s^T H s 近似目标,约束 ||s|| <= delta

  2. trcg 子问题求解:截断共轭梯度法,若步长超出信任域则截断到边界

  3. Hessian-vector 积fun_obj->Hv(d, Hd) 避免显式构建 Hessian,利用问题结构(如 X^T D X + I)高效计算

  4. 信任域更新:根据实际下降/预测下降比率 actred/prered 动态调整 delta

11.9.3 求解器选择逻辑

源码路径:sklearn/svm/_base.py - _get_liblinear_solver_type()(约 800-870 行)

def _get_liblinear_solver_type(multi_class, penalty, loss, dual):
    """Find the liblinear magic number for the solver."""
    _solver_type_dict = {
        "logistic_regression": {"l1": {False: 6}, "l2": {False: 0, True: 7}},
        "hinge": {"l2": {True: 3}},
        "squared_hinge": {"l1": {False: 5}, "l2": {False: 2, True: 1}},
        "epsilon_insensitive": {"l2": {True: 13}},
        "squared_epsilon_insensitive": {"l2": {False: 11, True: 12}},
        "crammer_singer": 4,
    }

    if multi_class == "crammer_singer":
        return _solver_type_dict[multi_class]
    elif multi_class != "ovr":
        raise ValueError(...)

    _solver_pen = _solver_type_dict.get(loss, None)
    if _solver_pen is None:
        error_string = "loss='%s' is not supported" % loss
    else:
        _solver_dual = _solver_pen.get(penalty, None)
        if _solver_dual is None:
            error_string = "The combination of penalty='%s' and loss='%s' is not supported" % (penalty, loss)
        else:
            solver_num = _solver_dual.get(dual, None)
            if solver_num is None:
                error_string = "The combination of penalty='%s' and loss='%s' are not supported when dual=%s" % (penalty, loss, dual)
            else:
                return solver_num
    raise ValueError(...)

代码解析:

三层嵌套字典映射:loss -> penalty -> dual -> solver_id。覆盖所有合法组合:

  • L2R_LR (0) / L2R_LR_DUAL (7):Logistic 回归原始/对偶

  • L2R_L2LOSS_SVC (2) / L2R_L2LOSS_SVC_DUAL (1):L2 损失 SVM 原始/对偶

  • L2R_L1LOSS_SVC_DUAL (3):L1 损失 SVM 对偶

  • L1R_L2LOSS_SVC (5) / L1R_LR (6):L1 正则化原始坐标下降

  • MCSVM_CS (4):Crammer-Singer 多分类

  • L2R_L2LOSS_SVR (11/12/13):SVR 变体


设计中的取舍:坐标下降 vs TRON

| 维度 | 坐标下降 | TRON 信任域 |

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

| 适用问题 | L1 正则、对偶 L1/L2 损失 | L2 正则原始问题 (LR, L2-SVC, L2-SVR) |

| 收敛速度 | 线性,大规模稀疏极快 | 超线性/二次,中小规模稠密快 |

| 内存 | 仅需存储 w, alpha, 梯度 | 需存 Hessian-vector 积临时数组 |

| 稀疏性利用 | 天然支持,仅遍历非零 | 需稠密梯度/Hv 计算 |

| 超参数敏感度 | 低,eps 控制 | 中等,信任域参数影响收敛 |

选择逻辑:_get_liblinear_solver_type 根据 (loss, penalty, dual) 自动派发。L1 正则必须用坐标下降(原始);L2 正则对偶用坐标下降;L2 正则原始用 TRON。

11.10 Cython 桥接层:零拷贝内存视图映射

11.10.1 稠密绑定:_libsvm.pyx

flowchart LR A[Python: np.ndarray float64_t[:, ::1]] --> B[Cython: const float64_t[:, ::1] X] B --> C[libsvm_helper.c: dense_to_libsvm] C --> D[libsvm: svm_node* 连续内存] D --> E[svm_train] E --> F[svm_model*] F --> G[libsvm_helper.c: set_model/copy_*] G --> H[Python: ndarray]

源码路径:sklearn/svm/_libsvm.pyx - fit()(约 50-200 行)

def fit(
    const float64_t[:, ::1] X,  # 内存视图:C 连续 float64 二维数组
    const float64_t[::1] Y,
    int svm_type=0,
    kernel='rbf',
    ...
):
    cdef svm_parameter param
    cdef svm_problem problem
    cdef svm_model *model
    cdef intp_t SV_len

    # ... 参数设置 ...

    set_problem(  # C 辅助函数:填充 svm_problem
        &problem,
        <char*> &X[0, 0],  # 直接传递内存视图指针,零拷贝
        <char*> &Y[0],
        <char*> &sample_weight[0],
        <intp_t*> X.shape,
        kernel_index,
    )

    cdef BlasFunctions blas_functions
    blas_functions.dot = _dot[double]  # 注入 BLAS dot 函数指针

    cdef int fit_status = 0
    with nogil:  # 释放 GIL,C++ 线程并行
        model = svm_train(&problem, &param, &fit_status, &blas_functions)

    # ... 复制结果到 NumPy 数组 ...
    return (support.base, support_vectors.base, ...)

代码解析:

  1. 内存视图const float64_t[:, ::1] 要求 C 连续、float64,直接映射 NumPy 数组缓冲区

  2. 指针传递&X[0, 0] 取首元素地址传给 C,无数据拷贝

  3. nogil 上下文:训练期间释放 GIL,允许多线程并行(OpenMP/BLAS)

  4. BLAS 注入blas_functions.dot = _dot[double] 将 SciPy BLAS dot 传给 C++ 核计算

11.10.2 稀疏绑定:_libsvm_sparse.pyx

源码路径:sklearn/svm/_libsvm_sparse.pyx - libsvm_sparse_train()(约 30-150 行)

def libsvm_sparse_train(
    int n_features,
    const float64_t[::1] values,
    const int32_t[::1] indices,
    const int32_t[::1] indptr,
    const float64_t[::1] Y,
    int svm_type, int kernel_type, ...
):
    # CSR 三元组直接传递
    problem = csr_set_problem(
        <char *> &values[0],
        <intp_t *> indices.shape,
        <char *> &indices[0],
        <intp_t *> indptr.shape,
        <char *> &indptr[0],
        <char *> &Y[0],
        <char *> &sample_weight[0],
        kernel_type,
    )
    # ... 训练 ...
    # 返回:support, support_vectors_(csr_matrix), sv_coef_data, intercept, ...

代码解析:

稀疏绑定直接传递 CSR 三元组 (data, indices, indptr)

  1. CSR 直通:避免转为稠密或 libsvm 内部链表格式的中间拷贝

  2. csr_set_problem:C 层将 CSR 转为 libsvm 的 svm_csr_node** 链表表示

  3. 返回 CSR 模型support_vectors_ 重建为 scipy.sparse.csr_matrixdual_coef_ 也是 CSR

11.10.3 liblinear 统一绑定:_liblinear.pyx

源码路径:sklearn/svm/_liblinear.pyx - train_wrap()(约 30-150 行)

def train_wrap(
    object X,  # 接受稠密 ndarray 或 scipy.sparse.csr_matrix
    const float64_t[::1] Y,
    bint is_sparse,
    int solver_type,
    ...
):
    if is_sparse:
        X_indices = X.indices
        X_indptr = X.indptr
        if X_has_type_float64:
            X_data_64 = X.data
            X_data_bytes_ptr = <char *> &X_data_64[0]
        else:
            X_data_32 = X.data
            X_data_bytes_ptr = <char *> &X_data_32[0]
        problem = csr_set_problem(...)  # CSR -> liblinear 稀疏格式
    else:
        X_as_1d_array = X.reshape(-1)
        problem = set_problem(...)  # 稠密 -> liblinear 稀疏格式

    # ... 训练 ...
    # 返回 Fortran 顺序权重矩阵 w (nr_class, nr_feature)
    return w.base, n_iter.base

代码解析:

  1. 统一入口:Python 层判断稀疏/稠密,统一转为 liblinear 的 feature_node** 稀疏格式

  2. float32/float64 兼容:自动处理两种精度

  3. 偏置处理bias > 0 时在特征末尾追加常数项

  4. Fortran 顺序输出w = np.empty((nr_class, nr_feature), order='F') 匹配 liblinear 列主序布局


11.11 C 辅助层:内存管理与数据转换

11.11.1 libsvm_helper.c:稠密模型重建

源码路径:sklearn/svm/src/libsvm/libsvm_helper.c - set_model()copy_*(约 50-200 行)

struct svm_model *set_model(struct svm_parameter *param, int nr_class,
                            char *SV, Py_ssize_t *SV_dims,
                            char *support, Py_ssize_t *support_dims,
                            Py_ssize_t *sv_coef_strides,
                            char *sv_coef, char *rho, char *nSV,
                            char *probA, char *probB)
{
    struct svm_model *model;
    double *dsv_coef = (double *) sv_coef;
    int i, m = nr_class * (nr_class-1)/2;

    if ((model = malloc(sizeof(struct svm_model))) == NULL) goto model_error;
    // ... 分配 nSV, label, sv_coef, rho ...

    if (param->kernel_type == PRECOMPUTED) {
        // precomputed: SV 只存索引
        model->SV = malloc(model->l * sizeof(struct svm_node));
        for (i=0; i<model->l; ++i) {
            model->SV[i].ind = ((int *) support)[i];
            model->SV[i].values = NULL;
        }
    } else {
        model->SV = dense_to_libsvm((double *) SV, SV_dims);  // 稠密连续内存映射
    }

    // sv_coef 是二维指针数组,需摊平拷贝
    for (i=0; i < model->nr_class-1; i++) {
        model->sv_coef[i] = dsv_coef + i*(model->l);
    }

    // rho 取反并去除 -0.0
    for (i=0; i<m; ++i) {
        t = model->rho[i];
        *ddata = (t != 0) ? -t : 0;
    }
    // ... probA/probB ...
    return model;
}

代码解析:

  1. 模型重建:预测时从 Python 属性重建 C svm_model 结构

  2. 稠密直通dense_to_libsvm 直接映射连续内存,sv_node.values 指向原数组

  3. sv_coef 摊平:Python 传入连续数组,C 层按行切片赋值给 sv_coef[i] 指针

  4. 错误回滚goto 链式错误处理,分配失败时依次释放已分配内存

11.11.2 libsvm_sparse_helper.c:稀疏 CSR 转换

源码路径:sklearn/svm/src/libsvm/libsvm_sparse_helper.c - csr_to_libsvm()csr_copy_SV()(约 30-150 行)

struct svm_csr_node **csr_to_libsvm(double *values, int* indices, int* indptr, int n_samples)
{
    struct svm_csr_node **sparse, *temp;
    for (i=0; i<n_samples; ++i) {
        n = indptr[i+1] - indptr[i];
        temp = malloc((n+1) * sizeof(struct svm_csr_node));
        for (j=0; j<n; ++j) {
            temp[j].value = values[k];
            temp[j].index = indices[k] + 1;  // libsvm 1-based 索引
            ++k;
        }
        temp[n].index = -1;  // 哨兵
        sparse[i] = temp;
    }
    return sparse;
}

int csr_copy_SV(char *data, Py_ssize_t *n_indices, char *indices, Py_ssize_t *n_indptr, char *indptr,
                struct svm_csr_model *model, int n_features)
{
    int i, j, k=0, index;
    double *dvalues = (double *) data;
    int *iindices = (int *) indices;
    int *iindptr = (int *) indptr;
    iindptr[0] = 0;
    for (i=0; i<model->l; ++i) {
        index = model->SV[i][0].index;
        for(j=0; index >=0; ++j) {
            iindices[k] = index - 1;  // 转回 0-based
            dvalues[k] = model->SV[i][j].value;
            index = model->SV[i][j+1].index;
            ++k;
        }
        iindptr[i+1] = k;
    }
    return 0;
}

代码解析:

  1. CSR -> 链表csr_to_libsvm 逐行分配 svm_csr_node 数组,索引 +1 适配 libsvm

  2. 链表 -> CSRcsr_copy_SV 遍历链表重建 CSR 三元组,索引 -1 恢复 0-based

  3. 内存管理:每行单独 malloc,释放时需逐行 freefree 指针数组

11.11.3 liblinear_helper.c:稠密转稀疏格式

源码路径:sklearn/svm/src/liblinear/liblinear_helper.c - dense_to_sparse()csr_to_sparse()(约 30-120 行)

static struct feature_node **dense_to_sparse(char *x, int double_precision,
        int n_samples, int n_features, int n_nonzero, double bias)
{
    struct feature_node **sparse;
    struct feature_node *T;
    n_nonzero += (have_bias+1) * n_samples;  // 预计算总非零数含偏置和哨兵
    T = malloc(n_nonzero * sizeof(struct feature_node));  # 一次性分配大块内存

    for (i=0; i<n_samples; ++i) {
        sparse[i] = T;
        for (j=1; j<=n_features; ++j) {
            if (val != 0) { T->value = val; T->index = j; ++T; }
        }
        if (have_bias) { T->value = bias; T->index = j; ++T; }
        T->index = -1; ++T;  // 哨兵
    }
    return sparse;
}

代码解析:

  1. 单次分配策略:预计算总非零数,一次性 malloc 大块内存,各样本指针指入其中,减少碎片和分配开销

  2. 偏置处理bias > 0 时在特征末尾追加 bias 作为合成特征

  3. 索引转换:Python 0-based → liblinear 1-based


11.12 跨平台随机数生成器:确定性优化路径

flowchart LR A[Python: random_state] --> B[_newrand.pyx: set_seed_wrap] B --> C[newrand.h: set_seed] C --> D[std::mt19937 mt_rand] D --> E[bounded_rand_int: Lemire 后处理] E --> F[libsvm/liblinear: SMO 工作集打乱/坐标下降打乱] F --> G[跨平台确定性结果]

11.12.1 Mersenne Twister + Lemire 后处理

源码路径:sklearn/svm/src/newrand/newrand.h(全文)

// Scikit-Learn-specific random number generator replacing `rand()` originally
// used in LibSVM / LibLinear, to ensure the same behaviour on windows-linux,
// with increased speed

// (1) Init a `mt_rand` object
std::mt19937 mt_rand(std::mt19937::default_seed);

// (2) public `set_seed()` function that should be used instead of `srand()` to set a new seed.
void set_seed(unsigned custom_seed) {
    mt_rand.seed(custom_seed);
}

// (3) New internal `bounded_rand_int` function, used instead of rand() everywhere.
inline uint32_t bounded_rand_int(uint32_t range) {
    // "LibSVM / LibLinear Original way" - make a 31bit positive
    // random number and use modulo to make it fit in the range
    // return abs( (int)mt_rand()) % range;

    // "Better way": tweaked Lemire post-processor
    // from http://www.pcg-random.org/posts/bounded-rands.html
    uint32_t x = mt_rand();
    uint64_t m = uint64_t(x) * uint64_t(range);
    uint32_t l = uint32_t(m);
    if (l < range) {
        uint32_t t = -range;
        if (t >= range) {
            t -= range;
            if (t >= range)
                t %= range;
        }
        while (l < t) {
            x = mt_rand();
            m = uint64_t(x) * uint64_t(range);
            l = uint32_t(m);
        }
    }
    return m >> 32;  // 高 32 位为无偏结果
}

代码解析:

  1. mt19937:C++11 标准 Mersenne Twister,周期 2^19937-1,跨平台行为一致

  2. Lemire 后处理m = x * range,高 32 位为结果,低 32 位 < range 时拒绝采样修正偏差

  3. 替代 rand() % range:后者有模偏差且 Windows/Linux rand() 实现不同,导致优化路径不一致

11.12.2 Python 侧包装

源码路径:sklearn/svm/_newrand.pyx(全文)

cdef extern from "newrand.h":
    void set_seed(unsigned int)
    unsigned int bounded_rand_int(unsigned int)

def set_seed_wrap(unsigned int custom_seed):
    set_seed(custom_seed)

def bounded_rand_int_wrap(unsigned int range_):
    return bounded_rand_int(range_)

代码解析:

Cython 直接暴露 C++ 函数,Python 测试可验证确定性:

  • set_seed_wrap(0); bounded_rand_int_wrap(100) → 固定值 54

  • KS 检验验证均匀分布


11.13 l1_min_c:L1 正则化参数下界计算

11.13.1 数学原理与实现

对于 L1 正则化线性 SVM/Logistic,当 C < l1_min_c 时最优解全为零。l1_min_c 给出非零解的最小 C,指导参数搜索下界。

源码路径:sklearn/svm/_bounds.py - l1_min_c()(约 20-80 行)

@validate_params(...)
def l1_min_c(X, y, *, loss="squared_hinge", fit_intercept=True, intercept_scaling=1.0):
    X = check_array(X, accept_sparse="csc")
    check_consistent_length(X, y)

    Y = LabelBinarizer(neg_label=-1).fit_transform(y).T  # (n_classes, n_samples), 标签 ±1
    # maximum absolute value over classes and features
    den = np.max(np.abs(safe_sparse_dot(Y, X)))  # ||Y^T X||_max

    if fit_intercept:
        bias = np.full((np.size(y), 1), intercept_scaling, dtype=...)
        den = max(den, abs(np.dot(Y, bias)).max())  # 考虑偏置项

    if den == 0.0:
        raise ValueError("Ill-posed l1_min_c calculation...")

    if loss == "squared_hinge":
        return 0.5 / den
    else:  # loss == 'log'
        return 2.0 / den

代码解析:

  1. 标签二值化LabelBinarizer(neg_label=-1) 将多类标签转为 ±1 矩阵 Y (n_classes, n_samples)

  2. 对偶梯度上界den = max_{c,j} |sum_i Y_{c,i} X_{i,j}|||Y^T X||_max

  3. 偏置项fit_intercept=True 时追加 intercept_scaling 列,更新 den

  4. 理论系数

    • squared_hinge0.5/den 来自对偶问题 KKT 条件 max |grad| <= 1/(2C)

    • log2.0/den 来自 Logistic 对偶梯度界 max |grad| <= 1/C

11.13.2 测试验证

源码路径:sklearn/svm/tests/test_bounds.py - test_l1_min_c()(约 20-50 行)

def check_l1_min_c(X, y, loss, fit_intercept=True, intercept_scaling=1.0):
    min_c = l1_min_c(X, y, loss=loss, fit_intercept=fit_intercept, intercept_scaling=intercept_scaling)

    clf = {"log": LogisticRegression(l1_ratio=1, solver="liblinear"),
           "squared_hinge": LinearSVC(loss="squared_hinge", penalty="l1", dual=False)}[loss]
    clf.fit_intercept = fit_intercept
    clf.intercept_scaling = intercept_scaling

    clf.C = min_c
    clf.fit(X, y)
    assert (np.asarray(clf.coef_) == 0).all()  # C=min_c 时全零解
    assert (np.asarray(clf.intercept_) == 0).all()

    clf.C = min_c * 1.01
    clf.fit(X, y)
    assert (np.asarray(clf.coef_) != 0).any() or (np.asarray(clf.intercept_) != 0).any()  # 略大则非零

代码解析:

测试验证边界正确性:C = min_c 得到零解,C = 1.01*min_c 得到非零解。


11.14 OneClassSVM/SVR/NuSVR 变体实现

11.14.1 OneClassSVM:无标签训练

源码路径:sklearn/svm/_classes.py - OneClassSVM.fit()(约 600-620 行)

class OneClassSVM(OutlierMixin, BaseLibSVM):
    _impl = "one_class"

    def fit(self, X, y=None, sample_weight=None):
        super().fit(X, np.ones(_num_samples(X)), sample_weight=sample_weight)  # y 全 1
        self.offset_ = -self._intercept_  # offset = -rho
        return self

    def decision_function(self, X):
        dec = self._decision_function(X).ravel()
        return dec

    def score_samples(self, X):
        return self.decision_function(X) + self.offset_  # 原始评分 = 决策值 + offset

    def predict(self, X):
        y = super().predict(X)
        return np.asarray(y, dtype=np.intp)  # +1 内点, -1 异常点

代码解析:

  1. 构造标签y = np.ones(n_samples) 全部视为正类

  2. libsvm 内部solve_one_class 设置标签 ones,上界 C = W[i]sum alpha = nu * sum C

  3. offset 定义decision_function = score_samples - offset_offset_ = -intercept_ = rho

11.14.2 SVR/NuSVR:回归变体

源码路径:sklearn/svm/src/libsvm/svm.cpp - solve_epsilon_svr()solve_nu_svr()(约 2100-2300 行)

// epsilon-SVR: 变量加倍 (alpha, alpha*)
static void solve_epsilon_svr(...) {
    int l = prob->l;
    double *alpha2 = new double[2*l];
    double *linear_term = new double[2*l];
    schar *y = new schar[2*l];
    double *C = new double[2*l];

    for (i=0; i<l; i++) {
        alpha2[i] = 0; linear_term[i] = param->p - prob->y[i]; y[i] = 1; C[i] = prob->W[i]*param->C;
        alpha2[i+l] = 0; linear_term[i+l] = param->p + prob->y[i]; y[i+l] = -1; C[i+l] = prob->W[i]*param->C;
    }
    Solver s;
    s.Solve(2*l, SVR_Q(...), linear_term, y, alpha2, C, ...);
    for (i=0; i<l; i++) alpha[i] = alpha2[i] - alpha2[i+l];  // alpha - alpha*
}

// Nu-SVR: Nu 约束 + 变量加倍
static void solve_nu_svr(...) {
    // C[i] = C[i+l] = W[i]*C
    // sum alpha = nu * sum C / 2
    // linear_term: -y[i], +y[i]
    // Solver_NU 求解
}

代码解析:

  • Epsilon-SVR:引入松弛变量,对偶变量加倍 (alpha, alpha*),标签 +1/-1,线性项 ±(p - y)

  • Nu-SVRnu 替代 epsilon,约束 sum alpha = nu * sum C / 2,用 Solver_NU


11.15 LinearSVC/LinearSVR:liblinear 驱动的线性模型

11.15.1 LinearSVC

源码路径:sklearn/svm/_classes.py - LinearSVC.fit()(约 150-250 行)

class LinearSVC(LinearClassifierMixin, SparseCoefMixin, BaseEstimator):
    def __init__(self, penalty="l2", loss="squared_hinge", dual="auto", ...):
        self.dual = dual
        self.penalty = penalty
        self.loss = loss
        ...

    @_fit_context(prefer_skip_nested_validation=True)
    def fit(self, X, y, sample_weight=None):
        X, y = validate_data(self, X, y, accept_sparse="csr", dtype=np.float64, order="C")
        check_classification_targets(y)
        self.classes_ = np.unique(y)

        _dual = _validate_dual_parameter(self.dual, self.loss, self.penalty, self.multi_class, X)

        self.coef_, self.intercept_, n_iter_ = _fit_liblinear(
            X, y, self.C, self.fit_intercept, self.intercept_scaling,
            self.class_weight, self.penalty, _dual, self.verbose,
            self.max_iter, self.tol, self.random_state, self.multi_class,
            self.loss, sample_weight=sample_weight,
        )
        self.n_iter_ = n_iter_.max().item()

        if self.multi_class == "crammer_singer" and len(self.classes_) == 2:
            self.coef_ = (self.coef_[1] - self.coef_[0]).reshape(1, -1)  # 二分类合并系数
            if self.fit_intercept:
                self.intercept_ = np.array([self.intercept_[1] - self.intercept_[0]])

        return self

代码解析:

  1. dual='auto' 自动选择_validate_dual_parameter 基于 n_samples vs n_features 和求解器支持性决定

  2. _fit_liblinear:统一预处理(标签编码、类别权重、偏置扩展),调用 liblinear.train_wrap

  3. Crammer-Singer 二分类后处理:合并两类系数差值作为二分类权重

11.15.2 LinearSVR

源码路径:sklearn/svm/_classes.py - LinearSVR.fit()(约 350-450 行)

class LinearSVR(RegressorMixin, LinearModel):
    def fit(self, X, y, sample_weight=None):
        X, y = validate_data(self, X, y, accept_sparse="csr", dtype=np.float64, order="C")
        penalty = "l2"  # SVR only accepts l2 penalty

        _dual = _validate_dual_parameter(self.dual, self.loss, penalty, "ovr", X)

        self.coef_, self.intercept_, n_iter_ = _fit_liblinear(
            X, y, self.C, self.fit_intercept, self.intercept_scaling,
            None, penalty, _dual, self.verbose, self.max_iter, self.tol,
            self.random_state, loss=self.loss, epsilon=self.epsilon,
            sample_weight=sample_weight,
        )
        self.coef_ = self.coef_.ravel()  # 回归展平为 1D
        self.n_iter_ = n_iter_.max().item()
        return self

代码解析:

  • 固定 penalty="l2"(liblinear SVR 仅支持 L2)

  • 损失函数:epsilon_insensitive (L1) 或 squared_epsilon_insensitive (L2)

  • 回归无类别权重,class_weight=None


11.16 测试验证体系:确保数值一致性与正确性

11.16.1 稠密/稀疏等价性测试

源码路径:sklearn/svm/tests/test_sparse.py - check_svm_model_equal()(约 30-70 行)

def check_svm_model_equal(dense_svm, X_train, y_train, X_test):
    sparse_svm = base.clone(dense_svm)
    dense_svm.fit(X_train.toarray(), y_train)
    sparse_svm.fit(X_train, y_train)

    assert sparse.issparse(sparse_svm.support_vectors_)
    assert sparse.issparse(sparse_svm.dual_coef_)
    assert_allclose(dense_svm.support_vectors_, sparse_svm.support_vectors_.toarray())
    assert_allclose(dense_svm.dual_coef_, sparse_svm.dual_coef_.toarray())
    if dense_svm.kernel == "linear":
        assert sparse.issparse(sparse_svm.coef_)
        assert_array_almost_equal(dense_svm.coef_, sparse_svm.coef_.toarray())
    assert_allclose(dense_svm.support_, sparse_svm.support_)
    assert_allclose(dense_svm.predict(X_test_dense), sparse_svm.predict(X_test))
    assert_array_almost_equal(dense_svm.decision_function(X_test_dense), sparse_svm.decision_function(X_test))

代码解析:

核心验证:同一模型在稠密/稀疏输入下产生数值一致的支持向量、对偶系数、预测、决策函数。

11.16.2 未排序索引测试

源码路径:sklearn/svm/tests/test_sparse.py - test_unsorted_indices()(约 100-150 行)

def test_unsorted_indices(csr_container):
    # reverse each row's indices
    def scramble_indices(X):
        new_data, new_indices = [], []
        for i in range(1, len(X.indptr)):
            row_slice = slice(*X.indptr[i-1:i+1])
            new_data.extend(X.data[row_slice][::-1])
            new_indices.extend(X.indices[row_slice][::-1])
        return csr_container((new_data, new_indices, X.indptr), shape=X.shape)

    X_sparse_unsorted = scramble_indices(X_sparse)
    assert not X_sparse_unsorted.has_sorted_indices

    unsorted_svc = svm.SVC(kernel="linear", probability=True, random_state=0).fit(X_sparse_unsorted, y)
    assert_allclose(unsorted_svc.coef_.toarray(), sorted_svc.coef_.toarray())
    assert_allclose(unsorted_svc.predict_proba(X_test_unsorted), sorted_svc.predict_proba(X_test))

代码解析:

验证 _sparse_fit 中的 X.sort_indices() 正确处理未排序 CSR,保证结果与排序后一致。

11.16.3 收敛警告测试

源码路径:sklearn/svm/tests/test_svm.py - test_libsvm_convergence_warnings()(约 1000-1050 行)

def test_libsvm_convergence_warnings(global_random_seed):
    a = svm.SVC(kernel=lambda x, y: np.dot(x, y.T), random_state=global_random_seed, max_iter=2)
    warning_msg = r"Solver terminated early \(max_iter=2\).  Consider pre-processing ..."
    with pytest.warns(ConvergenceWarning, match=warning_msg):
        a.fit(np.array(X), Y)
    assert np.all(a.n_iter_ == 2)

代码解析:

验证 max_iter 触发收敛警告,且 n_iter_ 正确记录迭代次数。


11.17 本章小结

本章我们深入剖析了 scikit-learn 中 SVM 模块的实现原理,从 Python 层的统一接口到底层 C++ 优化内核,再到 Cython 桥接层和跨平台随机数生成器,揭示了这一经典算法是如何通过模块化设计和性能优化在保持算法正确性的同时实现高效计算的。

我们首先理解了 BaseLibSVM 如何作为总调度室统一管理不同 SVM 变体的共享逻辑,特别是它如何通过 gamma 参数的三态解析和概率参数的弃用机制兼容不同使用场景;接着探讨了预测引擎如何通过稠密/稀疏双通道实现高效预测,并解释了二分类符号翻转的必要性;然后分析了 BaseSVC 如何在 OvO 基础上通过决策函数形状转换和打破平局机制支持更符合直觉的多分类输出;随后深入 libsvm 和 liblinear 的 C++ 内核,掌握了 SMO 求解器的工作集选择与收缩启发式,以及坐标下降和 TRON 信任域 Newton 法在线性 SVM 中的应用;进一步考察了 Cython 桥接层如何通过内存视图直通和 CSR 直通实现零拷贝数据传递,以及底层辅助层如何处理内存管理和数据转换;最后我们了解了跨平台随机数生成器如何通过 Mersenne Twister 和 Lemire 后处理确保确定性,以及 l1_min_c 如何为 L1 正则化 SVM 提供参数搜索的下界指导。

11.17.1 概念总结表

| 概念 | 解释 |

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

| BaseLibSVM | SVM 家族的抽象基类,统一管理核函数、支持向量、决策函数与稀疏/稠密数据分发 |

| SMO 求解器 | libsvm 核心优化算法,每次选取两个变量(工作集)求解二次规划子问题 |

| 工作集选择 | 两阶段启发式:先选违反 KKT 最大的变量,再选目标函数下降最快的配对变量 |

| 收缩启发式 | 将已收敛变量(alpha 在边界且梯度满足条件)从活跃集移出,加速训练 |

| 核缓存 (Cache) | LRU 策略缓存核矩阵列,避免重复计算,支持索引交换时的缓存同步 |

| 坐标下降 | liblinear 核心算法,逐坐标循环优化,利用稀疏性实现极快收敛 |

| TRON 信任域 | 信任域 Newton 方法,用二次模型近似目标函数,通过 CG 求解子问题 |

| Cython 桥接层 | _libsvm.pyx/_libsvm_sparse.pyx/_liblinear.pyx 实现零拷贝内存视图映射 |

| 双通道预测 | _dense_*_sparse_* 方法分离,Python 层通过 _sparse 标志自动分发 |

| 二分类符号翻转 | libsvm 内部用 +1/-1 编码,sklearn 需翻回用户标签顺序,fit/决策函数双重保证 |

| OvO/OvR 转换 | 训练用 OvO,预测可选 OvR 输出,_ovr_decision_function 实现投票转置信度 |

| Platt Scaling | 5折 CV 训练 Sigmoid 参数 probA/probB,将决策值映射为概率 |

| Mersenne Twister + Lemire | 跨平台确定性随机数生成器,修复 Windows/Linux 优化路径不一致问题 |

| l1_min_c | L1 正则化 SVM 非零解的 C 下界,基于对偶问题 KKT 条件推导 |

感谢你读到了这里,恭喜你,你已经完成了 支持向量机 模块的学习。

posted @ 2026-09-04 08:54  绝不原创的飞龙  阅读(3)  评论(0)    收藏  举报