def _dual_obj_2d(alpha_j, c, yi, yj, vi, vj, Kii, Kjj, Kij):
"""Objectif dual W restreint au sous-probleme (i, j), en fonction de alpha_j.
alpha_i est elimine par la contrainte d'egalite (alpha_i y_i + alpha_j y_j = c),
et v_i = sum_{k != i, j} alpha_k y_k K(x_i, x_k) (idem pour j). Ne restent que
les termes de W qui couplent les deux variables ; les autres sont constants et
s'annulent des lors qu'on compare deux valeurs de alpha_j.
"""
s = yi * yj
alpha_i = yi * (c - alpha_j * yj)
return (alpha_i + alpha_j
- 0.5 * (Kii * alpha_i * alpha_i + Kjj * alpha_j * alpha_j
+ 2.0 * s * Kij * alpha_i * alpha_j)
- alpha_i * yi * vi - alpha_j * yj * vj)
def _take_step(i, j, y, alpha, b, K, C, tol, eps, E_cache):
"""Sous-probleme 2D analytique (Platt, eq. 12.7-12.11).
Tente de mettre a jour (alpha[i], alpha[j]). Retourne (changed, b_new).
Si pas de progres : retourne (0, b) sans modification.
"""
if i == j:
return 0, b
alpha_i_old = alpha[i].copy()
alpha_j_old = alpha[j].copy()
yi, yj = y[i], y[j]
# bornes L, H imposees par 0 <= alpha <= C et alpha_i y_i + alpha_j y_j = const
if yi != yj:
L = max(0.0, alpha[j] - alpha[i])
H = min(C, C + alpha[j] - alpha[i])
else:
L = max(0.0, alpha[i] + alpha[j] - C)
H = min(C, alpha[i] + alpha[j])
if L >= H:
return 0, b
Kii, Kjj, Kij = K[i, i], K[j, j], K[i, j]
s = yi * yj
# courbure de W le long de la contrainte : eta = ||phi(x_i) - phi(x_j)||^2 >= 0
eta = Kii + Kjj - 2.0 * Kij
if eta > 0:
# W est strictement concave en alpha_j : maximum interieur -> formule close
alpha_j_new = alpha_j_old + yj * (E_cache[i] - E_cache[j]) / eta
alpha_j_new = min(H, max(L, alpha_j_new))
else:
# eta <= 0 : W est affine (eta = 0) ou convexe (eta < 0) en alpha_j, donc son
# maximum sur [L, H] est atteint a une BORNE. Platt (1998, App. A) evalue
# l'objectif aux deux bornes et retient la meilleure. Diviser par eta (ou par
# un max(eta, 1e-12)) donnerait un pas arbitrairement grand, dans une
# direction qui ne caracterise plus le maximum.
c = alpha_i_old * yi + alpha_j_old * yj
vi = E_cache[i] + yi - b - alpha_i_old * yi * Kii - alpha_j_old * yj * Kij
vj = E_cache[j] + yj - b - alpha_j_old * yj * Kjj - alpha_i_old * yi * Kij
W_L = _dual_obj_2d(L, c, yi, yj, vi, vj, Kii, Kjj, Kij)
W_H = _dual_obj_2d(H, c, yi, yj, vi, vj, Kii, Kjj, Kij)
if W_L > W_H:
alpha_j_new = L
elif W_L < W_H:
alpha_j_new = H
else:
alpha_j_new = alpha_j_old # objectif plat : aucun mouvement utile
if abs(alpha_j_new - alpha_j_old) < eps * (alpha_j_old + eps):
return 0, b
alpha_i_new = alpha_i_old + s * (alpha_j_old - alpha_j_new)
# mettre a jour alpha
alpha[i] = alpha_i_new
alpha[j] = alpha_j_new
# recalculer b (Platt eq. 12.7-12.8)
b1 = b - E_cache[i] - yi * (alpha_i_new - alpha_i_old) * Kii - yj * (alpha_j_new - alpha_j_old) * Kij
b2 = b - E_cache[j] - yi * (alpha_i_new - alpha_i_old) * Kij - yj * (alpha_j_new - alpha_j_old) * Kjj
if 0 < alpha_i_new < C:
b_new = b1
elif 0 < alpha_j_new < C:
b_new = b2
else:
b_new = 0.5 * (b1 + b2)
# mettre a jour le cache d'erreurs : E_k = sum_l alpha_l y_l K_lk + b - y_k
delta_i = alpha_i_new - alpha_i_old
delta_j = alpha_j_new - alpha_j_old
E_cache += yi * delta_i * K[:, i] + yj * delta_j * K[:, j]
# le nouveau biais decale TOUTES les erreurs cachees de (b_new - b)
E_cache += (b_new - b)
return 1, b_new
def smo_fit(X, y, C=1.0, gamma=1.0, tol=1e-3, max_passes=20, eps=1e-3, verbose=False):
"""SMO de Platt pour SVM binaire soft-margin avec noyau RBF.
Parametres
----------
X : (n, d) array d'entrainement
y : (n,) array d'etiquettes dans {-1, +1}
C : parametre de soft-margin (slack penalty)
gamma : parametre du noyau RBF
tol : tolerance du test KKT par exemple (section 5)
max_passes : nombre max de passes completes sans changement
eps : tolerance relative de rejet d'un pas sur alpha_j
Retourne
--------
dict avec cles : alpha (n,), b (scalaire), support_mask (bool n,),
n_iter (int, nombre d'appels a examine_example -- n par passe effectuee)
"""
n = X.shape[0]
alpha = np.zeros(n)
b = 0.0
K = rbf_kernel(X, X, gamma)
# cache d'erreurs E[i] = sum_j alpha_j y_j K[j,i] + b - y[i]
E_cache = -y.copy()
def examine_example(i):
"""Tente un pas sur l'exemple i : test KKT, puis cascade de Platt.
Le biais est une variable d'etat de smo_fit : chaque pas qui progresse
le met a jour, et le sous-probleme suivant doit repartir de la nouvelle
valeur (sinon les erreurs cachees E sont decalees de b, et la solution
duale converge vers un modele dont le score de decision est biaise).
"""
nonlocal b
Ei = E_cache[i]
yi = y[i]
ri = yi * Ei # ri = y_i (f(x_i) - y_i) = y_i f(x_i) - 1
# KKT : alpha=0 => r>=0, 0<alpha<C => r=0, alpha=C => r<=0 (a tol pres)
violates = (ri < -tol and alpha[i] < C) or (ri > tol and alpha[i] > 0)
if not violates:
return 0
# Heuristique 1 de Platt : parmi les exemples NON BORNES (0 < alpha_j < C),
# prendre celui qui maximise |Ei - Ej| (le mouvement attendu le plus ample).
non_bound = np.where((alpha > 0) & (alpha < C))[0]
non_bound = non_bound[non_bound != i]
if len(non_bound) > 0:
j = non_bound[np.argmax(np.abs(E_cache[non_bound] - Ei))]
ok, b = _take_step(i, j, y, alpha, b, K, C, tol, eps, E_cache)
if ok:
return 1
# Heuristique 2 (balayage complet) : tous les j != i, dans l'ordre. Garantit
# qu'aucune violation KKT ne reste inexploree faute de candidat non borne.
for j in range(n):
if j == i:
continue
ok, b = _take_step(i, j, y, alpha, b, K, C, tol, eps, E_cache)
if ok:
return 1
return 0
passes = 0
n_iter = 0
while passes < max_passes:
num_changed = 0
# --- balayage complet des n exemples ---
for i in range(n):
n_iter += 1
if examine_example(i):
num_changed += 1
if num_changed == 0:
passes += 1
else:
passes = 0
if verbose:
print(f' passe={passes}, num_changed={num_changed}, |SV|={(alpha > 1e-6).sum()}')
support_mask = alpha > 1e-6
return {
'alpha': alpha,
'b': float(b),
'support_mask': support_mask,
'n_iter': n_iter,
}
# --- auto-controle A : le biais appris est-il celui de la solution duale ? ---
# Jeu jouet lineairement separable (droite x1 + x2 = 1) : 4 points negatifs,
# 2 positifs. A l'optimum, tout point tel que 0 < alpha_k < C est SUR la marge,
# donc doit verifier y_k f(x_k) = 1 -- une condition qui echoue si b reste a 0.
_Xs = np.array([[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0], [2.0, 2.0], [3.0, 2.0]])
_ys = np.array([-1.0, -1.0, -1.0, -1.0, 1.0, 1.0])
_m = smo_fit(_Xs, _ys, C=1.0, gamma=1.0, tol=1e-4, max_passes=20)
_sm = _m['support_mask']
_on_margin = (_m['alpha'] > 1e-6) & (_m['alpha'] < 1.0 - 1e-6)
_f = predict_raw(_Xs, _Xs[_sm], _ys[_sm], _m['alpha'][_sm], _m['b'], 1.0)
_ecart = np.abs(_ys[_on_margin] * _f[_on_margin] - 1.0).max() if _on_margin.any() else float('nan')
_ecart_b0 = np.abs(_ys[_on_margin] * (_f[_on_margin] - _m['b']) - 1.0).max() if _on_margin.any() else float('nan')
print(f'test jouet : b appris = {_m["b"]:+.4f}, |SV| = {_sm.sum()}/{len(_ys)}, n_iter = {_m["n_iter"]}')
print(f'test jouet : ecart max |y_k f(x_k) - 1| sur les non-bornees = {_ecart:.2e}')
print(f"test jouet : le meme ecart si l'on force b = 0 = {_ecart_b0:.2e}")
# --- auto-controle B : la branche eta <= 0 doit MAXIMISER l'objectif dual ---
# Un noyau valide verifie K_ij <= sqrt(K_ii K_jj) (Cauchy-Schwarz), donc
# eta = K_ii + K_jj - 2 K_ij = ||phi_i - phi_j||^2 >= 0 : le regime eta <= 0 ne se
# rencontre qu'avec des points CONFONDUS (eta = 0 exact) ou du bruit flottant. On
# teste donc separement (1) que le raccourci 2D de la branche egale l'objectif dual
# complet, (2) qu'une paire confondue donne bien eta = 0, et (3) que la comparaison
# aux bornes maximise reellement -- sur un etat eta < 0 fabrique pour cibler
# l'arithmetique (aucun etat a noyau valide ne le produit).
def _W_full(alpha, yy, KK):
"""Objectif dual complet, recalcule SANS _dual_obj_2d (reference independante)."""
ay = alpha * yy
return float(alpha.sum() - 0.5 * ay @ KK @ ay)
# (1) etat reel coherent : le raccourci 2D doit egaler l'objectif complet
_rng = np.random.default_rng(7)
_Xw = _rng.normal(size=(7, 3))
_yw = np.where(_rng.random(7) > 0.5, 1.0, -1.0)
_Kw = rbf_kernel(_Xw, _Xw, 1.7)
_aw = np.round(_rng.random(7) * 0.5, 3)
_bw = 0.23
_Ew = (_Kw @ (_aw * _yw) + _bw) - _yw
_i, _j = 2, 5
_c = _aw[_i] * _yw[_i] + _aw[_j] * _yw[_j]
_vi = _Ew[_i] + _yw[_i] - _bw - _aw[_i] * _yw[_i] * _Kw[_i, _i] - _aw[_j] * _yw[_j] * _Kw[_i, _j]
_vj = _Ew[_j] + _yw[_j] - _bw - _aw[_j] * _yw[_j] * _Kw[_j, _j] - _aw[_i] * _yw[_i] * _Kw[_i, _j]
_w0 = _dual_obj_2d(_aw[_j], _c, _yw[_i], _yw[_j], _vi, _vj, _Kw[_i, _i], _Kw[_j, _j], _Kw[_i, _j])
_ecart_2d = 0.0
for _t in (0.0, 0.17, 0.41, 0.66):
_a2 = _aw.copy()
_a2[_j] = _t
_a2[_i] = _yw[_i] * (_c - _t * _yw[_j]) # contrainte d'egalite respectee
_ecart_2d = max(_ecart_2d,
abs((_W_full(_a2, _yw, _Kw) - _W_full(_aw, _yw, _Kw))
- (_dual_obj_2d(_t, _c, _yw[_i], _yw[_j], _vi, _vj,
_Kw[_i, _i], _Kw[_j, _j], _Kw[_i, _j]) - _w0)))
assert _ecart_2d < 1e-10, _ecart_2d
print(f'branche eta<=0 : raccourci 2D == objectif dual complet (ecart max {_ecart_2d:.1e})')
# (2) eta = 0 est atteignable avec un noyau RBF : deux points confondus
_Xp = np.array([[0.0, 0.0], [0.0, 0.0], [2.0, 1.0]])
_Kp = rbf_kernel(_Xp, _Xp, 1.0)
_eta_dup = _Kp[0, 0] + _Kp[1, 1] - 2.0 * _Kp[0, 1]
print(f'branche eta<=0 : deux points confondus -> eta = {_eta_dup:.1e} (cas reel de la branche)')
# (3) etat eta < 0 fabrique : le maximum doit tomber sur une borne, sans pas arbitraire
_Kd = np.array([[1.0, 1.5], [1.5, 1.0]])
_yd = np.array([1.0, 1.0])
for _Ed in (np.array([0.05, 0.0]), np.array([0.0, 0.0])):
_ad = np.array([0.05, 0.95])
_take_step(0, 1, _yd, _ad, 0.0, _Kd, 1.0, 1e-3, 1e-3, _Ed.copy())
_grille = np.linspace(0.0, 1.0, 20001)
_vals = np.array([_W_full(np.array([1.0 - _t, _t]), _yd, _Kd) for _t in _grille])
_argmax = _grille[_vals.argmax()]
print(f'branche eta<0 : E={_Ed} -> alpha_j choisi = {_ad[1]:.2f}, argmax de l objectif = {_argmax:.2f}')
assert abs(_ad[1] - _argmax) < 1e-3, (_ad[1], _argmax)
print('branche eta<=0 : OK (maximum pris a une borne)')