07 - Support Vector Machines: kernel trick e SVM non lineari

Hai due nuvole di punti a forma di mezzaluna, una dentro la curva dell'altra, e un classificatore lineare le separa con una retta che, sulle mezzelune dell'esempio più avanti, sbaglia 16 punti su 100 (accuratezza 0,840 sul test). Una Support Vector Machine (SVM) risolve il problema in due mosse: sceglie, fra le infinite rette che separano le classi, quella più lontana dai punti di entrambe; poi, con il kernel trick, calcola quella retta in uno spazio dove le mezzelune diventano separabili, senza mai costruire quello spazio.

Il prezzo sono due numeri da scegliere, C e gamma, e un tempo di addestramento che cresce almeno quadraticamente con i campioni. Questo articolo è il settimo della serie ML Fundamentals; ti servono la regressione logistica e la cross-validation, e dopo gli alberi e le foreste è un altro modello non lineare che incontri nella serie.

Ogni formula che leggi qui sotto è stata controllata con un numero: il margine di cinque punti calcolato a mano, l'uguaglianza fra obiettivo primale e duale, il kernel polinomiale confrontato con le feature esplicite, l'effetto di C e gamma, i tempi di addestramento misurati. Dove l'esperimento ha smentito una frase che circola nei manuali o in versioni precedenti di questo testo, vale l'esperimento e lo dico.

Il margine, calcolato a mano su cinque punti

Un iperpiano in due dimensioni è una retta, descritta da un vettore $w$ (la direzione perpendicolare alla retta) e da un numero $b$: i punti con $w \cdot x + b = 0$ stanno sulla retta, quelli con segno positivo da una parte, con segno negativo dall'altra. Con le etichette scritte come $y_i \in \{-1, +1\}$, un punto è classificato bene quando $y_i (w \cdot x_i + b) > 0$.

Il margine è la larghezza della fascia vuota attorno alla retta: se imponi $y_i (w \cdot x_i + b) \geq 1$ per ogni punto, la fascia ha larghezza $2 / \|w\|$, quindi massimizzare il margine equivale a minimizzare $\tfrac{1}{2}\|w\|^2$. I punti che toccano il bordo della fascia sono i support vector.

Prendi due punti negativi, $(0,0)$ e $(0,1)$, e tre positivi, $(2,0)$, $(2,1)$ e $(5,3)$. A occhio la retta migliore è $x_1 = 1$, cioè $w = (1, 0)$ e $b = -1$: il punto $(2,0)$ dà $1 \cdot 2 - 1 = 1$, il punto $(0,0)$ dà $-(0 - 1) = 1$, e il margine vale $2 / 1 = 2$.

Il listato chiede a libsvm (la libreria che sta sotto SVC) la stessa risposta. Ho messo C=1e6: un C enorme rende l'errore sul margine così costoso che il modello si comporta come una SVM a margine rigido, adatta a dati separabili come questi.

# Cinque punti nel piano: il margine si calcola a mano e poi si confronta con libsvm
import numpy as np
from sklearn.svm import SVC

X = np.array([[0, 0], [0, 1], [2, 0], [2, 1], [5, 3]], dtype=float)
y = np.array([-1, -1, 1, 1, 1])

clf = SVC(kernel="linear", C=1e6).fit(X, y)   # C enorme: quasi hard-margin
w = clf.coef_[0]
b = clf.intercept_[0]
print("w =", np.round(w, 6), " b =", round(b, 6))
print("margine 2/||w|| =", round(2 / np.linalg.norm(w), 6))
print("indici dei support vector:", clf.support_.tolist())
print("alpha * y dei support vector:", np.round(clf.dual_coef_[0], 6).tolist())

# w ricostruito dal duale: somma di alpha_i y_i x_i
w_duale = clf.dual_coef_[0] @ clf.support_vectors_
print("w dal duale =", np.round(w_duale, 6))
# vincolo del duale: somma di alpha_i y_i = 0
print("somma alpha_i y_i =", round(float(clf.dual_coef_[0].sum()), 9))
# y_i f(x_i) per tutti i punti: vale 1 sui support vector, di piu' altrove
print("y f(x) =", np.round(y * clf.decision_function(X), 4).tolist())

# Output:
# w = [1. 0.]  b = -1.0
# margine 2/||w|| = 2.0
# indici dei support vector: [1, 3]
# alpha * y dei support vector: [-0.5, 0.5]
# w dal duale = [1. 0.]
# somma alpha_i y_i = 0.0
# y f(x) = [1.0, 1.0, 1.0, 1.0, 4.0]

Il modello trova esattamente $w = (1, 0)$, $b = -1$ e margine 2, come il conto a mano. La riga y f(x) mostra che quattro punti su cinque valgono 1 (stanno sul bordo), mentre $(5,3)$ vale 4 ed è lontano. C'è però una sorpresa: i support vector restituiti sono solo due, gli indici 1 e 3, cioè $(0,1)$ e $(2,1)$.

Gli altri due punti sul bordo hanno coefficiente zero. Succede perché la soluzione del duale non è unica quando più punti stanno sul margine: due coefficienti $\alpha = 0.5$ bastano a produrre $w = 0.5\,(2,1) - 0.5\,(0,1) = (1, 0)$. Quindi «support vector» in scikit-learn significa «punto con coefficiente diverso da zero», che è un sottoinsieme dei punti che toccano il margine.

Soft margin, C e il problema duale

Con dati veri le classi si sovrappongono e nessuna retta separa tutto. Il soft margin permette a ogni punto di violare il margine di una quantità $\xi_i \geq 0$ (la slack) e la fa pagare: si minimizza $\tfrac{1}{2}\|w\|^2 + C \sum_i \xi_i$.

La slack ottima di un punto è la sua hinge loss, $\max(0,\, 1 - y_i (w \cdot x_i + b))$: zero se il punto è fuori dal margine, cresce linearmente se lo viola. C è il prezzo della violazione: C grande punisce l'errore e restringe il margine, C piccolo accetta più violazioni in cambio di un margine largo.

Lo stesso problema si può scrivere nei moltiplicatori di Lagrange $\alpha_i$ (il duale): si massimizza $\sum_i \alpha_i - \tfrac{1}{2}\sum_{i,j} \alpha_i \alpha_j y_i y_j\, x_i \cdot x_j$ con $0 \leq \alpha_i \leq C$ e $\sum_i \alpha_i y_i = 0$. I dati compaiono solo dentro i prodotti scalari $x_i \cdot x_j$, ed è la porta d'ingresso del kernel.

Siccome il problema è convesso, all'ottimo i due obiettivi hanno lo stesso valore. Il listato lo verifica su due lune con rumore, calcolando l'obiettivo primale dai pesi e dalle slack e quello duale dai coefficienti restituiti da dual_coef_, che contiene i prodotti $\alpha_i y_i$.

# Dualita' forte: sul soft margin l'obiettivo primale e quello duale coincidono all'ottimo
import numpy as np
from sklearn.datasets import make_moons
from sklearn.svm import SVC
from sklearn.metrics.pairwise import linear_kernel

X, y01 = make_moons(n_samples=200, noise=0.35, random_state=0)
y = 2 * y01 - 1          # etichette in {-1, +1}
C = 1.0
clf = SVC(kernel="linear", C=C, tol=1e-8).fit(X, y)

w, b = clf.coef_[0], clf.intercept_[0]
xi = np.maximum(0.0, 1.0 - y * (X @ w + b))               # slack = hinge loss
primale = 0.5 * w @ w + C * xi.sum()

a = clf.dual_coef_[0]                                     # alpha_i * y_i
Ksv = linear_kernel(clf.support_vectors_)
duale = np.abs(a).sum() - 0.5 * a @ Ksv @ a

print(f"support vector: {len(clf.support_)} su {len(X)}")
print(f"alpha = C (violano il margine): {int((np.abs(a) > C - 1e-6).sum())}")
print(f"0 < alpha < C (sul margine):    {int((np.abs(a) < C - 1e-6).sum())}")
print(f"obiettivo primale = {primale:.6f}")
print(f"obiettivo duale   = {duale:.6f}")
print(f"differenza        = {abs(primale - duale):.2e}")

# Output:
# support vector: 86 su 200
# alpha = C (violano il margine): 83
# 0 < alpha < C (sul margine):    3
# obiettivo primale = 82.499680
# obiettivo duale   = 82.499680
# differenza        = 1.16e-07

I due valori (82,499680) differiscono per 1e-07, cioè per la tolleranza del solver: il duale di scikit-learn e il primale che hai scritto a mano descrivono lo stesso ottimo. Il conteggio dei coefficienti mostra i tre regimi del duale:

• $\alpha_i = 0$: il punto è fuori dal margine e non entra nel modello (114 punti su 200);
• $0 < \alpha_i < C$: il punto sta esattamente sul bordo del margine (3 punti);
• $\alpha_i = C$: il punto viola il margine o è classificato male (83 punti).

Con rumore alto quasi tutti i support vector sono del terzo tipo: la soluzione è determinata da chi sta dentro il margine o dal lato sbagliato, non da chi sta sul bordo. Per la predizione bastano gli 86 support vector, $f(x) = \sum_i \alpha_i y_i\, x_i \cdot x + b$; gli altri 114 punti si possono buttare e il modello non cambia.

Il kernel trick, verificato con i numeri

Le mezzelune non si separano con una retta, ma si separano se aggiungi alle coordinate $(x_1, x_2)$ anche i loro quadrati e il prodotto. Il kernel è una funzione $K(x, z)$ che restituisce direttamente il prodotto scalare fra le versioni «arricchite» $\varphi(x) \cdot \varphi(z)$ senza calcolare $\varphi$.

Poiché il duale usa i dati solo attraverso $x_i \cdot x_j$, basta sostituirlo con $K(x_i, x_j)$ e la SVM lavora nello spazio arricchito. Il caso più piccolo che si può controllare: in due dimensioni $(x \cdot z)^2 = \varphi(x) \cdot \varphi(z)$ con $\varphi(x) = (x_1^2,\ \sqrt{2}\,x_1 x_2,\ x_2^2)$.

I valori fra parentesi sono i default di SVC nella documentazione di scikit-learn 1.9. Il kernel di default è RBF; gamma='scale' usa $1 / (n_{feature} \cdot \mathrm{Var}(X))$, dove $\mathrm{Var}(X)$ è la varianza di tutti i valori della matrice presi insieme.

Il listato che segue controlla quattro cose: l'identità del kernel polinomiale con le feature esplicite, la formula RBF scritta a mano, la semidefinitezza della matrice di Gram RBF (la condizione perché un kernel corrisponda a un vero prodotto scalare, per il teorema di Mercer) e un controesempio per il sigmoide.

# Il kernel trick controllato con i numeri
import numpy as np
from sklearn.datasets import make_moons
from sklearn.metrics.pairwise import rbf_kernel, polynomial_kernel

rng = np.random.default_rng(0)
X = rng.normal(size=(6, 2))

# 1) Kernel polinomiale omogeneo di grado 2 in 2 dimensioni:
#    (x . z)^2 = phi(x) . phi(z)  con  phi(x) = (x1^2, sqrt(2) x1 x2, x2^2)
def phi(x):
    return np.array([x[0] ** 2, np.sqrt(2) * x[0] * x[1], x[1] ** 2])

K_kernel = (X @ X.T) ** 2
F = np.array([phi(x) for x in X])
K_esplicito = F @ F.T
print("poly grado 2: max |K - phi.phi| =", f"{np.abs(K_kernel - K_esplicito).max():.2e}")
K_sk = polynomial_kernel(X, degree=2, gamma=1.0, coef0=0.0)
print("coincide con sklearn polynomial_kernel:", np.allclose(K_kernel, K_sk))

# 2) RBF: formula scritta a mano contro sklearn
gamma = 0.7
D2 = ((X[:, None, :] - X[None, :, :]) ** 2).sum(-1)
K_rbf = np.exp(-gamma * D2)
print("RBF a mano = sklearn rbf_kernel:", np.allclose(K_rbf, rbf_kernel(X, gamma=gamma)))

# 3) Mercer in pratica: la matrice di Gram RBF su 300 punti e' semidefinita positiva
P, _ = make_moons(n_samples=300, noise=0.3, random_state=0)
G = rbf_kernel(P, gamma=1.0)
lam = np.linalg.eigvalsh(G)
print(f"autovalore minimo della Gram RBF: {lam.min():.2e}  (>= 0 a meno dell'errore numerico)")

# 4) Il kernel sigmoide non garantisce la semidefinitezza: un controesempio numerico
S = np.tanh(2.0 * (X @ X.T) - 1.0)
print(f"autovalore minimo Gram sigmoide (gamma=2, coef0=-1): {np.linalg.eigvalsh(S).min():.3f}")

# Output:
# poly grado 2: max |K - phi.phi| = 8.88e-16
# coincide con sklearn polynomial_kernel: True
# RBF a mano = sklearn rbf_kernel: True
# autovalore minimo della Gram RBF: -1.17e-14  (>= 0 a meno dell'errore numerico)
# autovalore minimo Gram sigmoide (gamma=2, coef0=-1): -3.434

Il kernel polinomiale e il prodotto scalare delle feature esplicite coincidono a meno di 1e-15, cioè all'errore di arrotondamento. Per la RBF non puoi scrivere $\varphi$ (lo spazio ha dimensione infinita), ma la matrice di Gram su 300 punti ha autovalore minimo dell'ordine di -1e-14: nessun autovalore negativo oltre il rumore numerico, come garantisce la teoria per ogni $\gamma > 0$.

Il sigmoide invece, con questi parametri e questi sei punti, produce un autovalore di -3,43: la sua matrice non è semidefinita positiva e non corrisponde a nessun prodotto scalare. Non è un difetto di implementazione, è la ragione per cui in pratica si parte da lineare e RBF.

Che cosa fanno C e gamma al confine

Con il kernel RBF gamma regola la scala della similarità: $K(x, z)$ vale quasi 1 per punti vicini e quasi 0 per punti lontani, e gamma grande restringe il «vicino». Un gamma piccolo dà confini lisci, quasi lineari; un gamma grande dà un confine che si avvolge attorno ai singoli punti.

Per vederlo ho addestrato la SVM sullo stesso insieme di 200 punti di training (400 campioni, rumore 0,30, metà e metà) con diverse coppie di valori e ho misurato l'accuratezza sui 200 punti tenuti da parte. Il listato stampa anche il numero di support vector, che è il segnale più rapido di come sta andando il modello.

# Effetto di C e gamma sul confine, su due lune rumorose (dati sintetici, seme fisso)
import numpy as np
from sklearn.datasets import make_moons
from sklearn.model_selection import train_test_split
from sklearn.svm import SVC

X, y = make_moons(n_samples=400, noise=0.30, random_state=1)
Xtr, Xte, ytr, yte = train_test_split(X, y, test_size=0.5, random_state=1, stratify=y)

print("kernel   C       gamma   acc_train acc_test  n_SV  (su", len(Xtr), "train)")
righe = [("linear", 1, None), ("rbf", 1, 0.01), ("rbf", 1, 1.0), ("rbf", 1, 100.0),
         ("rbf", 0.01, 1.0), ("rbf", 1000, 1.0), ("rbf", 1000, 100.0)]
for kern, C, g in righe:
    kw = {} if g is None else {"gamma": g}
    m = SVC(kernel=kern, C=C, **kw).fit(Xtr, ytr)
    gs = "-" if g is None else g
    print(f"{kern:7s} {C:7g} {str(gs):7s} {m.score(Xtr, ytr):9.3f} {m.score(Xte, yte):8.3f} {m.n_support_.sum():5d}")

# Output:
# kernel   C       gamma   acc_train acc_test  n_SV  (su 200 train)
# linear        1 -           0.860    0.840    72
# rbf           1 0.01        0.830    0.825   144
# rbf           1 1.0         0.905    0.895    68
# rbf           1 100.0       0.950    0.830   185
# rbf        0.01 1.0         0.845    0.830   200
# rbf        1000 1.0         0.915    0.900    50
# rbf        1000 100.0       1.000    0.805   170

La lettura va fatta riga per riga, ricordando che con 200 punti di test un punto vale mezzo punto percentuale e differenze di 1-2 punti sono rumore. Con gamma=0.01 il confine è troppo rigido (0,830 sul training, 0,825 sul test, 144 support vector): è underfitting, a livello della retta (0,840, differenza dentro il rumore). Con gamma=1 e C=1 si ottiene 0,905 e 0,895 con 68 support vector, un buon equilibrio.

Con gamma=100 il training sale a 0,950 ma il test scende a 0,830: l'overfitting. L'ultima riga lo porta all'estremo, C=1000 e gamma=100: training perfetto, 1,000, e test a 0,805, il valore peggiore della tabella. Il numero di support vector è alto in entrambe le direzioni sbagliate: 185 con gamma troppo alto, 200 su 200 con C=0.01 (nessun punto riesce a stare fuori dal margine larghissimo).

Per questo una frazione alta di support vector è un sintomo, non una diagnosi: dice che qualcosa non va, non se sia underfitting o overfitting.

Lo scaling non è facoltativo

La distanza euclidea somma i quadrati delle differenze in ogni feature, quindi la feature misurata con i numeri più grandi comanda. Il dataset del tumore al seno incluso in scikit-learn (load_breast_cancer, 569 campioni e 30 feature) è un caso da manuale: la varianza delle 30 feature va da 7e-06 a 3e+05, oltre dieci ordini di grandezza.

Il confronto qui sotto usa la stessa SVM RBF con gli stessi 5 fold, con e senza StandardScaler. Lo scaler sta dentro la Pipeline, e non è un dettaglio: se lo applicassi a tutto il dataset prima della cross-validation, la media e la deviazione dei fold di test entrerebbero nel calcolo: una piccola fuga di informazione dal test al training.

# Senza scaling la distanza euclidea la decide la feature con l'unita' piu' grande
import numpy as np
from sklearn.datasets import load_breast_cancer
from sklearn.model_selection import StratifiedKFold, cross_val_score
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.svm import SVC

X, y = load_breast_cancer(return_X_y=True)
cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=0)

print("varianza per feature: min = %.2e, max = %.2e" % (X.var(axis=0).min(), X.var(axis=0).max()))
senza = cross_val_score(SVC(kernel="rbf"), X, y, cv=cv)
con = cross_val_score(make_pipeline(StandardScaler(), SVC(kernel="rbf")), X, y, cv=cv)
print(f"RBF senza scaling:  accuracy {senza.mean():.4f} (deviazione {senza.std():.4f})")
print(f"RBF con StandardScaler nella pipeline: accuracy {con.mean():.4f} (deviazione {con.std():.4f})")

# gamma='scale' = 1 / (n_features * X.var()): ricalcolato a mano
m = SVC(kernel="rbf").fit(X, y)
print("gamma usata da sklearn:", m._gamma, "  formula:", 1.0 / (X.shape[1] * X.var()))

# Output:
# varianza per feature: min = 6.99e-06, max = 3.24e+05
# RBF senza scaling:  accuracy 0.9210 (deviazione 0.0308)
# RBF con StandardScaler nella pipeline: accuracy 0.9771 (deviazione 0.0070)
# gamma usata da sklearn: 6.395533747973492e-07   formula: 6.395533747973492e-07

Senza scaling l'accuratezza è 0,921 con deviazione 0,031 fra i fold; con lo scaling 0,977 con deviazione 0,007. Il guadagno di oltre cinque punti non viene da un modello migliore, ma dal fatto che il modello senza scaling guarda quasi solo le poche feature con valori enormi.

L'ultima riga conferma la formula della documentazione: gamma='scale' vale $1 / (30 \cdot \mathrm{Var}(X))$, circa 6,4e-07 sui dati grezzi. Nota anche che _gamma è un attributo privato, usato qui solo per il controllo. La guida pratica degli autori di libsvm raccomanda lo scaling di ogni attributo in [0, 1] o [-1, +1], con gli stessi coefficienti su training e test.

Scegliere C e gamma con una griglia

Nessun valore di C e gamma è giusto in assoluto: dipendono dalla scala dei dati e dal rumore. Gli autori di libsvm consigliano una ricerca a griglia con cross-validation su valori che crescono in modo esponenziale, perché l'effetto di questi parametri è moltiplicativo.

Il listato divide il dataset in 398 righe di training e 171 di test, cerca nella griglia con 5 fold stratificati solo sul training e poi misura una volta sola sul test. Quel test non è mai stato visto dalla ricerca: è l'unico numero onesto, il punteggio della cross-validation è invece ottimistico perché è il massimo su sedici combinazioni.

# Scelta di C e gamma con una griglia logaritmica; lo scaling resta DENTRO la validazione
import numpy as np
from sklearn.datasets import load_breast_cancer
from sklearn.model_selection import GridSearchCV, StratifiedKFold, train_test_split
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.svm import SVC

X, y = load_breast_cancer(return_X_y=True)
Xtr, Xte, ytr, yte = train_test_split(X, y, test_size=0.3, random_state=0, stratify=y)

pipe = Pipeline([("scala", StandardScaler()), ("svc", SVC(kernel="rbf"))])
griglia = {"svc__C": [0.1, 1, 10, 100], "svc__gamma": [0.0003, 0.003, 0.03, 0.3]}
cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=0)
gs = GridSearchCV(pipe, griglia, cv=cv, scoring="accuracy").fit(Xtr, ytr)

print("migliori parametri:", gs.best_params_)
print(f"accuracy in CV (train 398 righe): {gs.best_score_:.4f}")
print(f"accuracy sul test tenuto da parte (171 righe): {gs.score(Xte, yte):.4f}")
svc = gs.best_estimator_.named_steps["svc"]
print(f"support vector: {svc.n_support_.sum()} su {len(Xtr)} = {svc.n_support_.sum() / len(Xtr):.1%}")

# tabella della CV: media per ogni coppia (righe = C, colonne = gamma)
medie = gs.cv_results_["mean_test_score"].reshape(4, 4)
print("\nC \\ gamma   " + "  ".join(f"{g:>7g}" for g in griglia["svc__gamma"]))
for C, riga in zip(griglia["svc__C"], medie):
    print(f"{C:>9g}   " + "  ".join(f"{v:7.3f}" for v in riga))

# Output:
# migliori parametri: {'svc__C': 100, 'svc__gamma': 0.003}
# accuracy in CV (train 398 righe): 0.9874
# accuracy sul test tenuto da parte (171 righe): 0.9591
# support vector: 35 su 398 = 8.8%
#
# C \ gamma    0.0003    0.003     0.03      0.3
#       0.1     0.628    0.897    0.950    0.628
#         1     0.900    0.965    0.985    0.937
#        10     0.965    0.980    0.980    0.937
#       100     0.985    0.987    0.977    0.937

La ricerca sceglie C=100 e gamma=0.003 con 0,9874 in cross-validation e 0,9591 sul test: la differenza di quasi tre punti dipende in parte dall'aver scelto il massimo fra molte prove e in parte dall'avere solo 171 righe di test (un errore in più o in meno sposta il valore di 0,6 punti). La tabella mostra una cresta diagonale: valori alti di C vanno con gamma basso e viceversa (0,985 con C=1 e gamma=0.03, 0,985 con C=100 e gamma=0.0003).

Le due manopole compensano in parte lo stesso effetto, quindi più coppie diverse danno un modello quasi equivalente. Il modello scelto usa 35 support vector su 398, l'8,8%. Ai bordi della tabella i valori scendono a 0,628, che è la frequenza della classe più numerosa nel training: un modello che predice sempre la stessa classe farebbe lo stesso.

Più di due classi

La SVM è nata per due classi. Con $K$ classi SVC addestra $K(K-1)/2$ classificatori, uno per ogni coppia, e la predizione è un voto a maggioranza (strategia uno contro uno). La documentazione di scikit-learn dice che il solver usa sempre questa strategia e che decision_function_shape cambia solo la forma dell'uscita di decision_function: 'ovo' restituisce un punteggio per coppia, 'ovr' (il default) lo trasforma in un punteggio per classe.

La strategia uno contro resto è un'altra cosa: $K$ modelli, ciascuno addestrato su tutti i campioni, che si ottiene con OneVsRestClassifier. Il listato usa le cifre di load_digits (1.797 immagini 8x8, 10 classi).

# SVC e' sempre uno-contro-uno; decision_function_shape cambia solo la FORMA dell'uscita
import numpy as np
from sklearn.datasets import load_digits
from sklearn.model_selection import train_test_split
from sklearn.multiclass import OneVsRestClassifier
from sklearn.svm import SVC

X, y = load_digits(return_X_y=True)
Xtr, Xte, ytr, yte = train_test_split(X / 16.0, y, test_size=0.3, random_state=0, stratify=y)

ovo = SVC(kernel="rbf", gamma=0.05, C=10, decision_function_shape="ovo").fit(Xtr, ytr)
ovr_forma = SVC(kernel="rbf", gamma=0.05, C=10, decision_function_shape="ovr").fit(Xtr, ytr)
print("classi K =", len(ovo.classes_), " coppie K(K-1)/2 =", len(ovo.classes_) * (len(ovo.classes_) - 1) // 2)
print("decision_function (ovo):", ovo.decision_function(Xte).shape)
print("decision_function (ovr):", ovr_forma.decision_function(Xte).shape)
print("stesse predizioni:", bool((ovo.predict(Xte) == ovr_forma.predict(Xte)).all()))
print(f"accuracy SVC (uno-contro-uno): {ovo.score(Xte, yte):.4f}")

# Uno-contro-resto esplicito: K modelli, ciascuno addestrato su TUTTI i campioni
ovr = OneVsRestClassifier(SVC(kernel="rbf", gamma=0.05, C=10)).fit(Xtr, ytr)
print("modelli addestrati da OneVsRestClassifier:", len(ovr.estimators_))
print(f"accuracy uno-contro-resto: {ovr.score(Xte, yte):.4f}")

# Output:
# classi K = 10  coppie K(K-1)/2 = 45
# decision_function (ovo): (540, 45)
# decision_function (ovr): (540, 10)
# stesse predizioni: True
# accuracy SVC (uno-contro-uno): 0.9926
# modelli addestrati da OneVsRestClassifier: 10
# accuracy uno-contro-resto: 0.9907

Con 10 classi i classificatori a coppie sono 45, e la decision_function ha 45 colonne con 'ovo' e 10 con 'ovr'. Le predizioni sono identiche in entrambi i casi: il parametro non cambia il modello. L'accuratezza è 0,9926 per l'uno-contro-uno e 0,9907 per l'uno-contro-resto, cioè una differenza di un campione su 540: nessuna delle due strategie vince qui.

La differenza di costo invece c'è: ogni classificatore a coppie vede solo i campioni di due classi, mentre ciascuno dei dieci uno-contro-resto vede tutti i campioni, quindi con molte classi la strategia a coppie fa molti modelli piccoli.

Quando la SVM non scala

La documentazione di SVC dice che il tempo di addestramento cresce almeno quadraticamente con il numero di campioni e che può diventare impraticabile oltre alcune decine di migliaia, e consiglia LinearSVC, SGDClassifier o una trasformazione Nystroem (che approssima il kernel RBF con un numero fisso di feature esplicite) prima di un modello lineare.

Il listato misura entrambe le strade su un problema sintetico a 20 feature, con campioni di training da 3.000 a 48.000 e 4.000 punti di test. Questi sono tempi: cambiano da una esecuzione all'altra, perché dipendono dalla macchina e da cos'altro sta girando, e per questo ho eseguito il listato cinque volte e riporto gli intervalli.

# Tempo di addestramento: SVC RBF contro Nystroem + LinearSVC (TEMPI: cambiano da una esecuzione all'altra)
import time
import numpy as np
from sklearn.datasets import make_classification
from sklearn.kernel_approximation import Nystroem
from sklearn.model_selection import train_test_split
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.svm import SVC, LinearSVC

X, y = make_classification(n_samples=52000, n_features=20, n_informative=8,
                           flip_y=0.05, class_sep=1.0, random_state=0)
Xtr, Xte, ytr, yte = train_test_split(X, y, test_size=4000, random_state=0, stratify=y)
sc = StandardScaler().fit(Xtr)
Xtr, Xte = sc.transform(Xtr), sc.transform(Xte)

print(f"{'n train':>8} {'SVC t(s)':>9} {'SVC acc':>8} {'Nys+Lin t(s)':>13} {'acc':>6}")
tempi = []
for n in [3000, 6000, 12000, 24000, 48000]:
    t0 = time.perf_counter()
    m = SVC(kernel="rbf", gamma="scale").fit(Xtr[:n], ytr[:n])
    t_svc = time.perf_counter() - t0
    t0 = time.perf_counter()
    a = make_pipeline(Nystroem(gamma=1 / (20 * 1.0), n_components=300, random_state=0),
                      LinearSVC(dual="auto", max_iter=5000)).fit(Xtr[:n], ytr[:n])
    t_nys = time.perf_counter() - t0
    tempi.append(t_svc)
    print(f"{n:8d} {t_svc:9.2f} {m.score(Xte, yte):8.3f} {t_nys:13.2f} {a.score(Xte, yte):6.3f}")
esp = np.polyfit(np.log([3000, 6000, 12000, 24000, 48000]), np.log(tempi), 1)[0]
print(f"esponente stimato t ~ n^p per SVC: p = {esp:.2f}")

# Output:
#  n train  SVC t(s)  SVC acc  Nys+Lin t(s)    acc
#     3000      0.10    0.947          0.14  0.934
#     6000      0.34    0.951          0.31  0.941
#    12000      1.33    0.955          0.62  0.943
#    24000      4.52    0.959          1.29  0.943
#    48000     24.53    0.963          2.71  0.947
# esponente stimato t ~ n^p per SVC: p = 1.95

L'accuratezza è deterministica, i tempi no. Su AMD Ryzen 5 5600H con 2 thread, in cinque esecuzioni (con altri processi attivi sulla macchina), SVC ha impiegato 0,08-0,17 s a 3.000 campioni e 25-46 s a 48.000: sedici volte i campioni, da circa 240 a 580 volte il tempo, con un esponente stimato fra 1,95 e 2,22, cioè circa quadratico. La documentazione dice che il tempo cresce almeno quadraticamente e che la complessità del solver QP sta fra quadratica e cubica; una stima su cinque punti ha un'incertezza di quell'ordine.

Nystroem con 300 componenti seguito da LinearSVC è più lento o pari a 3.000-6.000 campioni (0,10-0,34 s), poi vince: 2,6-3,3 s a 48.000, da nove a sedici volte più veloce, ma perdendo circa 1,6 punti di accuratezza (0,947 contro 0,963). Non è quindi vero che l'approssimazione «accelera sempre di 10-50 volte»: sotto le decine di migliaia di campioni la SVM esatta è competitiva, e il vantaggio compare dove la SVM esatta diventa lenta. Questo problema sintetico non è un benchmark generale, serve a mostrare la forma della curva.

Probabilità da una SVM: Platt scaling

Una SVM produce un punteggio con segno, la distanza dall'iperpiano nello spazio del kernel, non una probabilità. Il Platt scaling addestra una regressione logistica a una sola variabile sul punteggio, usando punteggi calcolati fuori campione, e lo trasforma in una probabilità. In SVC lo faceva il parametro probability=True, che con scikit-learn 1.9 è deprecato e sarà rimosso nella 1.11: la documentazione indica CalibratedClassifierCV(SVC(), ensemble=False).

La stessa documentazione avverte che probability=True rallenta fit perché usa internamente una cross-validation a 5 fold e che predict_proba può essere incoerente con predict. Il listato verifica l'avviso, calibra con il metodo raccomandato e confronta il Brier score con una sigmoide non adattata.

# Probabilita' da una SVM: Platt scaling con CalibratedClassifierCV (SVC(probability=True) e' deprecato in 1.9)
import warnings
import numpy as np
import sklearn
from sklearn.calibration import CalibratedClassifierCV
from sklearn.datasets import load_breast_cancer
from sklearn.metrics import brier_score_loss
from sklearn.model_selection import train_test_split
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.svm import SVC

print("scikit-learn", sklearn.__version__)
X, y = load_breast_cancer(return_X_y=True)
Xtr, Xte, ytr, yte = train_test_split(X, y, test_size=0.3, random_state=0, stratify=y)
base = make_pipeline(StandardScaler(), SVC(kernel="rbf", C=10, gamma=0.003))

# 1) il vecchio modo: che cosa dice la versione installata
with warnings.catch_warnings(record=True) as w:
    warnings.simplefilter("always")
    make_pipeline(StandardScaler(), SVC(probability=True)).fit(Xtr, ytr)
print("avviso:", str(w[0].message)[:110] if w else "nessuno")

# 2) il modo indicato dalla documentazione: Platt scaling (sigmoid) su score fuori campione
cal = CalibratedClassifierCV(base, method="sigmoid", cv=5, ensemble=False).fit(Xtr, ytr)
p = cal.predict_proba(Xte)[:, 1]
print(f"Brier score sul test (piu' basso = meglio): {brier_score_loss(yte, p):.4f}")

# 3) l'alternativa ingenua: sigmoide dello score grezzo, senza adattarla ai dati
score = base.fit(Xtr, ytr).decision_function(Xte)
p_ing = 1 / (1 + np.exp(-score))
print(f"Brier con 1/(1+exp(-score)) non adattata:    {brier_score_loss(yte, p_ing):.4f}")
print(f"accordo fra predict e soglia 0.5 su predict_proba: {((p > 0.5) == base.predict(Xte)).mean():.3f}")

# Output:
# scikit-learn 1.9.0
# avviso: The `probability` parameter was deprecated in 1.9 and will be removed in version 1.11. Use `CalibratedClassifi
# Brier score sul test (piu' basso = meglio): 0.0315
# Brier con 1/(1+exp(-score)) non adattata:    0.0421
# accordo fra predict e soglia 0.5 su predict_proba: 0.994

Il Brier score è l'errore quadratico medio fra la probabilità prevista e l'esito 0/1, e più basso è meglio. La calibrazione adattata ai dati dà 0,0315, la sigmoide ingenua $1/(1+e^{-s})$ applicata al punteggio grezzo dà 0,0421: la forma della trasformazione va appresa, non assunta.

Nota anche che la soglia 0,5 sulla probabilità calibrata e predict concordano nel 99,4% dei casi, non nel 100%: un campione su 171 cambia classe, ed è esattamente l'incoerenza che la documentazione segnala. Se ti serve solo un punteggio di confidenza, e non una probabilità, usa decision_function e salta la calibrazione.

Esercizio svolto

Calcola a mano la matrice del kernel RBF con $\gamma = 0.5$ sui punti $a=(0,0)$, $b=(1,1)$ e $c=(3,0)$, e dì quale coppia è la più simile secondo il kernel. Il procedimento: calcola le distanze al quadrato, moltiplicale per $-\gamma$, esponenzia. Poi controlla con il listato.

# Esercizio: matrice del kernel RBF con gamma = 0.5 su tre punti
import numpy as np

P = np.array([[0.0, 0.0], [1.0, 1.0], [3.0, 0.0]])
gamma = 0.5
d2 = ((P[:, None, :] - P[None, :, :]) ** 2).sum(-1)    # distanze euclidee al quadrato
K = np.exp(-gamma * d2)
print("distanze^2:\n", d2)
print("K:\n", np.round(K, 4))

# Output:
# distanze^2:
#  [[0. 2. 9.]
#  [2. 0. 5.]
#  [9. 5. 0.]]
# K:
#  [[1.     0.3679 0.0111]
#  [0.3679 1.     0.0821]
#  [0.0111 0.0821 1.    ]]

Che cosa sai fare adesso

Dopo questo articolo sei in grado di fare cinque cose, ciascuna con un controllo numerico a portata di mano:

• calcolare a mano il margine di un insieme di punti e verificarlo con SVC;
• spiegare perché obiettivo primale e duale coincidono e che cosa vuol dire che un punto è un support vector;
• controllare che un kernel sia davvero un prodotto scalare;
• scegliere C e gamma con una griglia logaritmica dentro una Pipeline con lo scaler;
• riconoscere underfitting e overfitting dalla differenza fra training e test e sapere quando abbandonare la SVM esatta.

Le cinque capacità si tengono: il margine e il duale spiegano che cosa fa il modello, il kernel spiega come diventa non lineare, e le ultime due voci dicono come usarlo senza farsi ingannare dai numeri di una sola esecuzione o da una griglia che sceglie il massimo fra troppe prove.

Non abbiamo trattato la regressione a vettori di supporto, i kernel su stringhe o grafi, e la derivazione di SMO: la guida di scikit-learn è il punto di partenza, e la serie Machine Learning Engineer ha un articolo sul margine massimo e il trucco del kernel. Nel prossimo articolo passiamo dal supervisionato al non supervisionato con il K-Means e il clustering gerarchico; la valutazione dei modelli, con le metriche che qui abbiamo usato solo come accuratezza, è nell'articolo sulle metriche. Se le feature sono molte, la PCA è un modo di ridurle prima di un modello a kernel, e come per ogni altro passo l'effetto va misurato; portare una SVM in produzione, con i controlli sullo scaling che qui sono già comparsi, è l'argomento del deployment.