sklearn.svm.SVCNavigation : << 2.7b-SMO-From-Scratch | Index | 2.8-Theorie-PAC >>
Kernel : Python 3 · Durée estimée : ~30 min
Concept-phare : 2.7b a écrit la boucle SMO à la main (~150 lignes, 21 560 itérations, 1,19 s) et l’a validée contre
sklearn.svm.SVC. Ce notebook inverse la perspective : le solveur industriel devient le sujet.SVCest l’enveloppe Python de LIBSVM (Chang & Lin, 2011) — du C++ compilé, un working set au second ordre (Fan, Chen & Lin, 2005), du shrinking et un cache de noyau. On mesure ce que ce packaging achète : exactitude, temps, vecteurs supports — puis on ouvre le dual que SVC ne rend pas (dual_coef_), on chiffre le gap de dualité et les violations KKT côté LIBSVM avec la même machinerie que 2.7b, et on actionne les molettes que la boîte noire expose (tol).
Position dans la série : bloc B.5 de l’EPIC optimisation convexe — le pendant SOTA de 2.7b (bloc A.1). Les blocs B.6 (Lasso-SOTA) et B.7 (CVXPY) sont déjà couverts par 2.11b — Proximal Operators, qui confronte ISTA/FISTA à sklearn.linear_model et à cvxpy sur le même problème Lasso — ce notebook ne les double pas.
SVC ne retourne pas : alpha complet depuis dual_coef_ + support_, et vérifier la contrainte d’égalité Σ α_i y_i = 0.tol, et en exercice, GridSearchCV sur C × gamma.NumPy et scikit-learn suffisent ; le noyau RBF vient de sklearn.metrics.pairwise (le côté SOTA jusqu’au bout).
Référence. Chang, C.-C. & Lin, C.-J. (2011), LIBSVM: a library for support vector machines, ACM TIST 2(27). Fan, R.-E., Chen, P.-H. & Lin, C.-J. (2005), Working set selection using second order information for training support vector machines, JMLR 6.
numpy=2.2.6 scikit-learn=1.8.0Tout tourne sur CPU en quelques secondes. La comparaison avec 2.7b n’a de sens que sur le même dataset, la même graine et les mêmes hyperparamètres — c’est la section suivante.
Le protocole de 2.7b section 7 est recopié exactement : make_moons(200, noise=0.25, random_state=42), labels en {-1, +1}, split stratifié 70/30 (random_state=42), C = 1.0, gamma = 2.0, tol = 1e-3. La seule variable expérimentale est le solveur : LIBSVM ici, la SMO de Platt maison dans 2.7b.
X, y_pm = make_moons(n_samples=200, noise=0.25, random_state=42)
y = np.where(y_pm == 0, -1, 1).astype(float)
X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.30, random_state=42, stratify=y)
C, GAMMA, TOL = 1.0, 2.0, 1e-3
print(f"train : {len(y_tr)} points test : {len(y_te)} points (2.7b : identique)")
print(f"C={C} gamma={GAMMA} tol={TOL:.0e} (2.7b : identique)")train : 140 points test : 60 points (2.7b : identique)
C=1.0 gamma=2.0 tol=1e-03 (2.7b : identique)Mêmes points, même split, mêmes hyperparamètres qu’en 2.7b : à la fin du notebook, toute différence de résultat ne pourra être attribuée qu’au solveur.
Un fit LIBSVM sur 140 points dure une fraction de milliseconde : un perf_counter autour d’un appel unique affiche 0.00 s (c’est ce que montre la sortie de 2.7b section 7). Pour comparer des temps de cet ordre, on répète le fit à froid (nouvel objet SVC à chaque itération, pas de cache d’objets Python entre les mesures) et on rapporte la moyenne.
clf = SVC(C=C, kernel="rbf", gamma=GAMMA, tol=TOL, max_iter=20000)
clf.fit(X_tr, y_tr)
acc_skl = accuracy_score(y_te, clf.predict(X_te))
n_sv = int(clf.n_support_.sum())
N_REPETS = 50
t0 = time.perf_counter()
for _ in range(N_REPETS):
SVC(C=C, kernel="rbf", gamma=GAMMA, tol=TOL, max_iter=20000).fit(X_tr, y_tr)
t_skl = (time.perf_counter() - t0) / N_REPETS
print(f"sklearn SVC (LIBSVM) : accuracy test = {acc_skl:.3f}, |SV| = {n_sv}, "
f"temps moyen = {t_skl * 1e3:.2f} ms (moyenne de {N_REPETS} fits a froid)")
print("2.7b SMO from scratch (commite) : accuracy test = 0.967, |SV| = 48, temps = 1.19 s, n_iter = 21560")sklearn SVC (LIBSVM) : accuracy test = 0.967, |SV| = 48, temps moyen = 0.47 ms (moyenne de 50 fits a froid)
2.7b SMO from scratch (commite) : accuracy test = 0.967, |SV| = 48, temps = 1.19 s, n_iter = 21560SVC ne rend pasSVC n’expose pas le vecteur alpha complet — il rend dual_coef_ (les y_i · α_i des seuls supports, concaténés) et support_ (leurs indices dans le train). Le dual complet se reconstruit : alpha nul hors supports, |dual_coef_| sur les supports. Avec lui, tous les instruments de 2.7b section 10 se pointent sur LIBSVM :
||w||² reconstruit depuis le dual : Σ_ij (y_i α_i)(y_j α_j) K(x_i, x_j) ;D(α) = Σ α_i − ½ ||w||² ;P = ½ ||w||² + C Σ ξ_i, les ξ lus sur decision_function ;P − D et violations KKT résiduelles.# --- reconstruction du dual complet --------------------------------------
alpha = np.zeros(len(y_tr))
alpha[clf.support_] = np.abs(clf.dual_coef_[0])
print(f"alpha reconstruit : {int((alpha > 1e-12).sum())} non nuls "
f"(n_support_ = {n_sv}) somme = {alpha.sum():.4f}")
print(f"contrainte d'egalite sum(alpha_i * y_i) = {float((alpha * y_tr).sum()):.2e}")
# --- objectifs primal et dual sur la MEME solution -------------------------
ay = clf.dual_coef_[0]
K_ss = pairwise.rbf_kernel(clf.support_vectors_, clf.support_vectors_, gamma=GAMMA)
w2 = float(ay @ K_ss @ ay)
D_libsvm = float(alpha.sum() - 0.5 * w2)
f_tr = clf.decision_function(X_tr)
xi = np.maximum(0.0, 1.0 - y_tr * f_tr)
P_libsvm = float(0.5 * w2 + C * xi.sum())
gap = P_libsvm - D_libsvm
# --- violations KKT residuelles (meme forme que 2.7b section 10) -----------
viol_low = np.where(alpha < C - 1e-12, np.maximum(0.0, 1.0 - y_tr * f_tr), 0.0)
viol_high = np.where(alpha > 1e-12, np.maximum(0.0, y_tr * f_tr - 1.0), 0.0)
V = float(viol_low.sum() + viol_high.sum())
print(f"||w||^2 depuis le dual = {w2:.4f}")
print(f"objectif dual D(alpha) = {D_libsvm:.4f}")
print(f"objectif primal P(w, b, xi) = {P_libsvm:.4f}")
print(f"gap de dualite P - D = {gap:.2e} (>= 0 : dualite faible)")
print(f"violation KKT residuelle sum V_i = {V:.2e}")
print(f" pire violation ponctuelle = {max(viol_low.max(), viol_high.max()):.2e}")
print(f"repartition des alpha : =0 -> {(alpha <= 1e-6).sum()}, "
f"0<alpha<C -> {((alpha > 1e-6) & (alpha < C - 1e-6)).sum()}, "
f"=C -> {(alpha >= C - 1e-6).sum()}")
print("2.7b SMO (commite) : D = 30.3960, gap = 3.3e-03, sum V = 6.3e-03")alpha reconstruit : 48 non nuls (n_support_ = 48) somme = 40.1958
contrainte d'egalite sum(alpha_i * y_i) = 2.22e-16
||w||^2 depuis le dual = 19.5995
objectif dual D(alpha) = 30.3961
objectif primal P(w, b, xi) = 30.3981
gap de dualite P - D = 1.99e-03 (>= 0 : dualite faible)
violation KKT residuelle sum V_i = 4.12e-03
pire violation ponctuelle = 5.22e-04
repartition des alpha : =0 -> 92, 0<alpha<C -> 16, =C -> 32
2.7b SMO (commite) : D = 30.3960, gap = 3.3e-03, sum V = 6.3e-03tol prèsΣ α_i y_i ≈ 10⁻¹⁶) : la reconstruction depuis dual_coef_ est fidèle, et c’est l’exercice 2 qui la fait vérifier à fond.3.3 × 10⁻³) — même ordre de grandeur, epsilon plus petit. Aucun des deux solveurs n’atteint l’optimum exact : les deux s’arrêtent à tol = 10⁻³ près, chacun avec son critère. La complémentarité KKT (points hors marge à α = 0, points à α = C hors ou sur la marge, points dans la marge à 0 < α < C) se lit sur la répartition des alpha comme en 2.7b.D(α) de LIBSVM en contre-épreuve (écart 4 × 10⁻⁵ avec sa propre solution). Ici on le mesure depuis l’intérieur : c’est la même mathématique, conditionnée par un arrêt approché assumé.Les colonnes 2.7b citent les valeurs committées de son exécution (même machine, même dataset) ; les colonnes LIBSVM sont mesurées ci-dessus. Les lignes de code comptent le code spécifique au solveur : rbf_kernel + sous-problème 2D + cascade de Platt + smo_fit pour 2.7b, contre l’appel SVC(...).fit(...) ici.
SMO_D, SMO_GAP, SMO_ACC, SMO_SV, SMO_T, SMO_ITER, SMO_LOC = 30.3960, 3.3e-3, 0.967, 48, 1.19, 21560, 150
print("=" * 96)
print(f"{'solveur':<26}{'D(alpha)':>10}{'gap P-D':>10}{'accuracy':>10}"
f"{'|SV|':>6}{'temps':>10}{'n_iter':>9}{'LOC':>6}")
print("-" * 96)
print(f"{'LIBSVM via SVC (mesure)':<26}{D_libsvm:>10.4f}{gap:>10.1e}{acc_skl:>10.3f}"
f"{n_sv:>6}{t_skl * 1e3:>8.2f}ms{'n/d':>9}{2:>6}")
print(f"{'SMO Platt maison (2.7b)':<26}{SMO_D:>10.4f}{SMO_GAP:>10.1e}{SMO_ACC:>10.3f}"
f"{SMO_SV:>6}{SMO_T:>8.2f}s{SMO_ITER:>9}{SMO_LOC:>6}")
print("=" * 96)
print("n_iter LIBSVM : non expose par l'API sklearn (boite noire) -- 2.7b maison le mesure.")================================================================================================
solveur D(alpha) gap P-D accuracy |SV| temps n_iter LOC
------------------------------------------------------------------------------------------------
LIBSVM via SVC (mesure) 30.3961 2.0e-03 0.967 48 0.47ms n/d 2
SMO Platt maison (2.7b) 30.3960 3.3e-03 0.967 48 1.19s 21560 150
================================================================================================
n_iter LIBSVM : non expose par l'API sklearn (boite noire) -- 2.7b maison le mesure.predict/decision_function/GridSearchCV/pipeline).alpha, le mécanisme du working set — 2.7b rend tout cela visible parce que tout y est écrit à la main. Diagnostiquer un LIBSVM qui ne converge pas (max_iter atteint, warnings de shrinking) exige de savoir ce qui vit dessous.tol est le contrat de convergence que LIBSVM signe avec vous : le solveur s’arrête dès que la violation KKT la plus forte passe sous ce seuil. Plus il est serré, plus le dual monte (vers le même optimum), plus le temps monte — légèrement : sur 140 points, tout reste sub-milliseconde.
print(f"{'tol':>8}{'temps moyen':>14}{'accuracy':>10}{'|SV|':>6}{'D(alpha)':>10}{'gap P-D':>10}")
for tol_i in (1e-2, 1e-3, 1e-4, 1e-5):
clf_i = SVC(C=C, kernel="rbf", gamma=GAMMA, tol=tol_i, max_iter=200000)
t0 = time.perf_counter()
for _ in range(N_REPETS):
SVC(C=C, kernel="rbf", gamma=GAMMA, tol=tol_i, max_iter=200000).fit(X_tr, y_tr)
t_i = (time.perf_counter() - t0) / N_REPETS
clf_i.fit(X_tr, y_tr)
acc_i = accuracy_score(y_te, clf_i.predict(X_te))
a_i = np.zeros(len(y_tr))
a_i[clf_i.support_] = np.abs(clf_i.dual_coef_[0])
ay_i = clf_i.dual_coef_[0]
K_i = pairwise.rbf_kernel(clf_i.support_vectors_, clf_i.support_vectors_, gamma=GAMMA)
w2_i = float(ay_i @ K_i @ ay_i)
D_i = float(a_i.sum() - 0.5 * w2_i)
f_i = clf_i.decision_function(X_tr)
P_i = float(0.5 * w2_i + C * np.maximum(0.0, 1.0 - y_tr * f_i).sum())
print(f"{tol_i:>8.0e}{t_i * 1e3:>11.2f}ms{acc_i:>10.3f}"
f"{int(clf_i.n_support_.sum()):>6}{D_i:>10.4f}{P_i - D_i:>10.1e}") tol temps moyen accuracy |SV| D(alpha) gap P-D
1e-02 0.47ms 0.967 49 30.3953 1.6e-02
1e-03 0.44ms 0.967 48 30.3961 2.0e-03
1e-04 0.46ms 0.967 48 30.3961 1.7e-04
1e-05 0.49ms 0.967 48 30.3961 2.4e-05L’objectif dual grimpe doucement quand tol se resserre et le gap rétrécit — l’optimum ne change pas de nature, seule la précision du contrat change. Sur ce dataset, l’accuracy est stable sur toute la grille et les supports se stabilisent dès tol = 10⁻³ (49 supports à 10⁻², 48 ensuite) : c’est l’indice qu’un tol lâche suffit quand le problème est bien séparé, et qu’un tol serré est de l’argent gaspillé — sauf sur les problèmes dégénérés où les supports instables font la différence. Le prix reste ici sub-milliseconde ; c’est sur les gros datasets qu’il se paie réellement.
Trois mécanismes séparent LIBSVM de la boucle de Platt de 2.7b — aucun ne change l’optimum, tous changent le chemin :
η), qui réduit le nombre d’itérations d’un ordre de grandeur sur les grands problèmes. C’est exactement le working set alternatif que l’exercice 3 de 2.7b fait explorer.α_i qui n’ont pas bougé depuis longtemps sont provisoirement retirés du problème actif (ils seront presque tous à 0 ou à C à l’optimum) — le sous-problème rétrécit pendant la résolution, et est re-vérifié à la fin.K en entier — impossible au-delà de quelques dizaines de milliers de points.La contrepartie du packaging : n_iter invisible, trajectoire des α invisible, paramètres internes (cache_size, shrinking) réduits à des molettes opaques.
Trois exercices progressifs : le premier fait varier C sur le solveur SOTA et relie le budget de marge au nombre de supports ; le deuxième vérifie en profondeur la reconstruction du dual de la section 3 ; le troisième est le geste de production complet — la recherche d’hyperparamètres par validation croisée.
# Exercice 1 : l'effet de C sur le nombre de vecteurs supports, cote LIBSVM
# TODO etudiant : completez sweep_C(valeurs) qui, pour chaque C de valeurs,
# entraine un SVC(C=C, kernel='rbf', gamma=2.0, tol=1e-3) sur (X_tr, y_tr),
# et retourne la liste des couples (n_supports, accuracy_test). Faites varier
# C dans [0.01, 0.1, 1, 10, 100] et expliquez la MONOTONIE du nombre de
# supports : pourquoi un C petit fabrique beaucoup de supports a la borne C,
# et un C grand n'en fabrique presque plus ?
def sweep_C(valeurs):
# Indice : clf.n_support_.sum() pour le compte ; accuracy_score(y_te, ...).
# Etape 1 : boucle sur valeurs, fit a froid.
# Etape 2 : collecter (n_supports, accuracy) de chaque fit.
result = None # TODO etudiant
return result
print("Exercice 1 a completer -- sweep_C([0.01, 1, 100]) retourne", sweep_C([0.01, 1, 100]))Exercice 1 a completer -- sweep_C([0.01, 1, 100]) retourne NoneLa section 3 reconstruit alpha depuis dual_coef_ et constate Σ α_i y_i ≈ 10⁻¹⁶. Une vérification plus forte compare les scores de décision recalculés depuis le dual reconstruit à ceux de clf.decision_function.
# Exercice 2 : le dual reconstruit reproduit-il decision_function ?
# TODO etudiant : ecrivez ecart_reconstruction(clf, X) qui recalcule le score
# de decision f(x) = sum_i alpha_i y_i K(x_i, x) + b sur CHAQUE ligne de X
# depuis (alpha reconstruit, y_tr, X_tr, clf.intercept_), en utilisant
# pairwise.rbf_kernel(X, X_tr, gamma), et retourne l'ecart max
# |f_recalcule - clf.decision_function(X)|. Resultat attendu : ~1e-15.
def ecart_reconstruction(clf, X):
# Indice : les alpha non supports valent 0 -- la somme se reduit aux
# supports, mais LAISSER la somme complete sur X_tr est plus simple et
# numeriquement equivalent. Le biais est clf.intercept_[0].
# Etape 1 : reconstruire alpha (section 3).
# Etape 2 : f = (alpha * y_tr) @ rbf_kernel(X, X_tr).T + b.
# Etape 3 : ecart max avec clf.decision_function(X).
result = None # TODO etudiant
return result
print("Exercice 2 a completer -- ecart_reconstruction() retourne",
ecart_reconstruction(clf, X_te))Exercice 2 a completer -- ecart_reconstruction() retourne NoneGridSearchCV sur C × gammaLe SOTA ne s’arrête pas au solveur : l’écosystème enchaîne la validation croisée sur la grille d’hyperparamètres, avec le scoring de votre choix.
# Exercice 3 : grille C x gamma par validation croisee
# TODO etudiant : avec sklearn.model_selection.GridSearchCV, cherchez la
# meilleure combinaison C dans [0.1, 1, 10] x gamma dans [0.5, 2, 8] en
# 5-fold sur (X_tr, y_tr), scoring='accuracy'. Affichez best_params_,
# best_score_, et comparez l'accuracy test du meilleur modele a celle
# mesuree ici (0.967) avec C=1, gamma=2. La grille trouve-t-elle mieux ?
def meilleure_grille():
# Indice : GridSearchCV(SVC(kernel='rbf'), {'C': [...], 'gamma': [...]},
# cv=5).fit(X_tr, y_tr) ; le meilleur estimateur est deja refitte,
# accessible via .best_estimator_.
# Etape 1 : poser la grille et fitter.
# Etape 2 : retourner (best_params, best_score_cv, accuracy_test).
result = None # TODO etudiant
return result
print("Exercice 3 a completer -- meilleure_grille() retourne", meilleure_grille())Exercice 3 a completer -- meilleure_grille() retourne NoneLe geste SOTA vaut trois ordres de grandeur de temps gagnés, mais ce notebook a surtout montré que la boîte noire s’ouvre : le dual complet se reconstruit depuis dual_coef_, le gap de dualité et les violations KKT se mesurent sur LIBSVM exactement comme sur la SMO maison, et la tolérance d’arrêt est un contrat explicite — pas une promesse d’optimalité exacte. Les deux solveurs rendent le même classifieur au tol près : c’est la définition opérationnelle d’un problème bien posé.
Couverture de l’EPIC optimisation convexe (#16061) : ce notebook couvre le bloc B.5 (SVM-SOTA). Les blocs B.6 (Lasso-SOTA : sklearn.linear_model coordinate descent, CV du λ) et B.7 (CVXPY déclaratif, solvers ECOS/SCS) sont couverts par 2.11b — Proximal Operators from scratch, qui confronte ISTA/FISTA aux implémentations scikit-learn et à cvxpy sur le même problème Lasso. Les blocs A.1 (SMO) et A.3 (proximal) sont livrés par 2.7b et 2.11b ; le bloc A.2 (ADMM) reste ouvert.
Pour aller plus loin : 2.7b pour la boucle SMO complète et ses diagnostics visuels ; 2.7 pour le cadre SVM/kernel ; 2.11b pour l’autre grande famille de l’optimisation convexe non différentiable.