scikit-learn：根据对偶系数计算截距

Question

我正在尝试使用对偶系数计算 SVM 分类器（线性核）的截距。我使用 SVC 类而不是 LinearSVC，因为后者没有 Dual_coef_ 属性。

这是代码：

from sklearn import svm
from sklearn.utils import shuffle
from sklearn.inspection import DecisionBoundaryDisplay
import matplotlib.pyplot as plt
import numpy as np
from numpy.testing import assert_array_almost_equal


seed = np.random.seed(123)

# Data
m = 1000

c2 = (-1, 1)
c1 = (2, 2.5)

r1 = 1.5
r2 = 3.5

rand_1_theta = np.random.uniform(0, 1, size=m)
rand_1_radius = np.random.uniform(0, 1, size=m)
ds1_x = (c1[0] + np.sqrt(rand_1_radius * r1) * np.cos(rand_1_theta * 2 * np.pi))
ds1_y = (c1[1] + np.sqrt(rand_1_radius * r1) * np.sin(rand_1_theta * 2 * np.pi))

rand_2_theta = np.random.uniform(0, 1, size=m)
rand_2_radius = np.random.uniform(0, 1, size=m)
ds2_x = (c2[0] + np.sqrt(rand_2_radius * r2) * np.cos(rand_2_theta * 2 * np.pi))
ds2_y = (c2[1] + np.sqrt(rand_2_radius * r2) * np.sin(rand_2_theta * 2 * np.pi))

X = np.stack((np.concatenate((ds1_x, ds2_x)), 
              np.concatenate((ds1_y, ds2_y))), 
              axis=1)
y = np.concatenate((np.zeros(m), np.ones(m)))
rgb = [(1,0,0) if val==0 else (0,0,1) for val in y]

# Preprocessing
X, y, rgb = shuffle(X, y, rgb, random_state=123)


# SVM
model = svm.SVC(probability=True, kernel='linear', random_state=123)
model.fit(X, y)


# Evaluation
support_vectors = model.support_vectors_
support_indexes = model.support_
coef = model.coef_
dual_coeff = model.dual_coef_
intercept = model.intercept_

primal_coeff = np.matmul(dual_coeff, support_vectors)


print(primal_coeff )
print(coef)

y_support = y[support_indexes]
all_bias = y_support - np.dot(support_vectors, primal_coeff.squeeze())

intercept_estimate_1 = np.mean(all_bias)
intercept_estimate_2 = np.min(all_bias)

print(intercept)
print(intercept_estimate_1)
print(intercept_estimate_2)

注意使用脚本末尾打印的平均值和最小值进行的截距估计。 据我所知，截距是所有截距估计的平均值。至少，大多数书籍都是这样介绍的。但 scikit-learn 似乎正在使用最小值。我不明白为什么，文档也没有解释它是如何完成的。

谁能解释一下？

Answer 1

截距既不是所有偏差的平均值也不是最小值。

相反，对于 libsvm（以及 sklearn）来说，它是 所有自由向量的 y* 梯度的平均值。您可以在 LIBSVM 论文，4.1.5.

中看到这一点

因为在预测中，只有支持向量很重要，而 libsvm 不会保存所有 coef，我认为如果不深入研究核心，你无法轻松检查此属性。

但是如果你真的对此感兴趣，你可以检查 libsvm 代码来计算偏差（=-rho）这里

double Solver::calculate_rho()
{
    double r;
    int nr_free = 0;
    double ub = INF, lb = -INF, sum_free = 0;
    for(int i=0;i<active_size;i++)
    {
        double yG = y[i]*G[i];

        if(is_upper_bound(i))
        {
            if(y[i]==-1)
                ub = min(ub,yG);
            else
                lb = max(lb,yG);
        }
        else if(is_lower_bound(i))
        {
            if(y[i]==+1)
                ub = min(ub,yG);
            else
                lb = max(lb,yG);
        }
        else
        {
            ++nr_free;
            sum_free += yG;
        }
    }

    if(nr_free>0)
        r = sum_free/nr_free;
    else
        r = (ub+lb)/2;

    return r;
}

Answer 2

我想补充已接受的答案。我最近有一篇post正是关于这个问题的。 StackOverflow 似乎不能很好地处理数学，所以我没有在这里粘贴帖子源代码。因此，请参阅我的原始帖子以进行澄清。

我的主张

我的主张是，根据 LIBSVM 论文 4.1.5，截距 IS 所有截距估计的平均值：

在 LIBSVM 中，为了数值稳定性，我们对所有这些值进行平均。

其中“所有这些值”指的是截距估计。在我的帖子中，我表明 负 y 倍梯度（在 @hychou 的回答中）只是表示截距的另一种形式。这部分你是对的。

你错在哪里

基本上，

• 截距估计不是取自所有支持向量（@hychou的答案在这部分是正确的）。仅考虑自由 alpha，并删除有界的 alpha，特别是上界的 alpha。您的代码没有删除那些 alpha 和支持向量。
• SVM 公式中的 y 取值 +1 和 -1，但在您的代码中，

  y_support

  取 1 和 0。是的，当输入

  SVC

  时，y 取 1 和 0；但是当按照数学公式手动计算截距时，+1 和 -1 是你的朋友。

我重写了你的代码，它在我的计算机上完美运行，至少使用你提供的种子。 我用

###

突出显示了您应该注意的地方：

from sklearn import svm
from sklearn.utils import shuffle
from sklearn.inspection import DecisionBoundaryDisplay
import matplotlib.pyplot as plt
import numpy as np
from numpy.testing import assert_array_almost_equal


seed = np.random.seed(123)

# Data
m = 1000

c2 = (-1, 1)
c1 = (2, 2.5)

r1 = 1.5
r2 = 3.5

rand_1_theta = np.random.uniform(0, 1, size=m)
rand_1_radius = np.random.uniform(0, 1, size=m)
ds1_x = (c1[0] + np.sqrt(rand_1_radius * r1) * np.cos(rand_1_theta * 2 * np.pi))
ds1_y = (c1[1] + np.sqrt(rand_1_radius * r1) * np.sin(rand_1_theta * 2 * np.pi))

rand_2_theta = np.random.uniform(0, 1, size=m)
rand_2_radius = np.random.uniform(0, 1, size=m)
ds2_x = (c2[0] + np.sqrt(rand_2_radius * r2) * np.cos(rand_2_theta * 2 * np.pi))
ds2_y = (c2[1] + np.sqrt(rand_2_radius * r2) * np.sin(rand_2_theta * 2 * np.pi))

X = np.stack((np.concatenate((ds1_x, ds2_x)), 
              np.concatenate((ds1_y, ds2_y))), 
              axis=1)
y = np.concatenate((np.zeros(m), np.ones(m)))
rgb = [(1,0,0) if val==0 else (0,0,1) for val in y]

# Preprocessing
X, y, rgb = shuffle(X, y, rgb, random_state=123)


# SVM
model = svm.SVC(probability=True, kernel='linear', random_state=123)
model.fit(X, y)


# Evaluation
support_vectors = model.support_vectors_
support_indexes = model.support_
coef = model.coef_
dual_coeff = model.dual_coef_
intercept = model.intercept_

primal_coeff = np.matmul(dual_coeff, support_vectors)


print(primal_coeff )
print(coef)

### BEGIN NOTICE HERE ###
y_support = y[support_indexes] * 2 - 1
S = np.ravel(np.abs(dual_coeff)) < 1
all_bias = y_support[S] - np.dot(support_vectors[S], primal_coeff.squeeze())
### END NOTICE HERE ###

intercept_estimate_1 = np.mean(all_bias)
intercept_estimate_2 = np.min(all_bias)

print(intercept)
print(intercept_estimate_1)
print(intercept_estimate_2)