2.11e — CVXPY : le problème déclaré, le solveur choisi

Série 02-ML-Cours (Data Science With Agents) — Bloc B.7 de l’issue #16061 (Optimisation convexe avancée from scratch). Le Bloc A a construit les méthodes from scratch (proximal en 2.11b, ADMM en 2.11d), le début du Bloc B a pris le côté SOTA bibliothèque (sklearn en 2.11c). Ce notebook prend le troisième chemin : la modélisation déclarative. On n’écrit plus l’algorithme du tout — on écrit le problème, et on choisit le solveur par son nom.

• la même formulation Lasso écrite avec cvxpy, sur le même problème que 2.11b/2.11c (seed 42), donc comparable nombre à nombre ;
• ce que la couche de modélisation fabrique avant d’appeler le solveur : la reformulation canonique, lue sur les données réellement transmises ;
• un problème, plusieurs solveurs — CLARABEL (point intérieur), SCS (premier ordre), OSQP (QP) : même optimum, tarifs différents ;
• la dualité lue dans les variables duales : le certificat KKT des blocs A et B sort directement du solveur, sans le recoder ;
• la confrontation avec ADMM (2.11d) et avec les méthodes from scratch : les uns et les autres convergent-ils vers le même optimum global ?

Chaîne pédagogique : 2.11-Regularisation-Sparse-LASSO (la règle) -> 2.11b-Proximal-Operators-From-Scratch (bloc A.3) -> 2.11c-Lasso-SOTA-Comparison (bloc B.6) -> 2.11d-Optimisation-ADMM-From-Scratch (bloc A.2) -> 2.11e (ce notebook, bloc B.7).

1. Le même problème que les blocs A et B

La règle de méthode ne change pas d’un notebook à l’autre : pour que les nombres se comparent, le problème doit être identique. On reprend donc exactement la génération de 2.11b (cellules 9-10), telle que reprise par 2.11c (cellule 2) : \(n = 200\) mesures, \(p = 500\) dimensions, \(k = 30\) coefficients non nuls, bruit additif \(0{,}01\), seed 42, matrice sous-gaussienne normalisée \(A/\sqrt{n}\) (lignes de norme unitaire — ce détail fixe la correspondance exacte \(\alpha = \lambda/n\) avec sklearn, cf 2.11c \(\S3\)).

Deux échelles de \(\lambda\) cohabitent dans ce notebook, et il faut les garder distinctes :

Symbole	Valeur	Rôle
\(\lambda_{DJ} = \sigma\sqrt{2\log p}/\sqrt{n}\)	\(0{,}0025\)	l’a priori de Donoho-Johnstone
\(\lambda = 5\lambda_{DJ}\)	\(0{,}0125\)	le point de fonctionnement retenu dans 2.11b/2.11c
\(\lambda_{\max} = \|A^{\top}b\|_\infty\)	\(3{,}19\)	le seuil d’annulation exact (cf \(\S7\))

import numpy as np
import time
import cvxpy as cp

np.random.seed(42)

# --- Generation identique a 2.11b (cellules 9-10) et 2.11c (cellule 2) ---
n, p, k = 200, 500, 30
noise_level = 0.01
rng = np.random.default_rng(42)
A = rng.standard_normal((n, p)) / np.sqrt(n)
support = rng.choice(p, size=k, replace=False)
x_true = np.zeros(p)
x_true[support] = rng.standard_normal(k)
b = A @ x_true + noise_level * rng.standard_normal(n)

# Les trois echelles de lambda (cf tableau du S1)
sigma = noise_level
lam_theory = sigma * np.sqrt(2.0 * np.log(p)) / np.sqrt(n)
lam = 5.0 * lam_theory
lam_max = np.max(np.abs(A.T @ b))


def objectif(x, lam_=None):
    """Objectif Lasso brut : 0.5||Ax - b||^2 + lam*||x||_1."""
    l = lam if lam_ is None else lam_
    return 0.5 * np.linalg.norm(A @ x - b) ** 2 + l * np.abs(x).sum()


print(f'Probleme : n={n}, p={p}, k={k} (densite {k/p:.2%}), bruit={noise_level}, seed=42')
print(f'lambda_DJ = {lam_theory:.6f} | lambda retenu (5x) = {lam:.6f}')
print(f'lambda_max = max_j |A_j^T b| = {lam_max:.6f}')
print(f'cvxpy {cp.__version__} | solveurs installes : {sorted(cp.installed_solvers())}')

Probleme : n=200, p=500, k=30 (densite 6.00%), bruit=0.01, seed=42
lambda_DJ = 0.002493 | lambda retenu (5x) = 0.012465
lambda_max = max_j |A_j^T b| = 3.191173
cvxpy 1.9.2 | solveurs installes : ['CLARABEL', 'GLOP', 'HIGHS', 'OSQP', 'PDLP', 'SCIPY', 'SCS']

2. Le modeling déclaratif : dire le problème, pas la méthode

En 2.11b/2.11c/2.11d, le code est l’algorithme : une boucle de balayages pour la coordinate descent, une boucle de pas proximaux pour ISTA, trois mises à jour alternées pour ADMM. Ici, la ligne de code est l’énoncé mathématique :

\[\min_x \; \tfrac12 \|Ax - b\|_2^2 + \lambda \|x\|_1\]

Rien de plus n’est écrit — pas de pas de gradient, pas de variable auxiliaire, pas de critère d’arrêt. cvxpy sépare trois responsabilités qui étaient confondues dans le from scratch :

1. le modèle — l’expression 0.5*cp.sum_squares(A @ x - b) + lam*cp.norm1(x), qui est l’objectif tel qu’un humain l’écrirait ;
2. la vérification de convexité — le ruleset DCP (Disciplined Convex Programming), qui certifie à la compilation que le problème est convexe, sans jamais l’évaluer ;
3. la résolution — un solveur nommé, appelé par son nom, qui reçoit une forme canonique et rend une solution.

Le point 2 est le vrai apport pédagogique : le from scratch ne se trompe pas de convexité parce qu’on a fait l’analyse à la main ; le déclaratif l’exige. Une inversion de signe dans l’objectif n’est pas un résultat bizarre au bout de 2000 itérations — c’est un refus immédiat, avant tout calcul.

# Le probleme, ecrit tel quel. Aucune instruction d'algorithme.
x_model = cp.Variable(p)
cout = 0.5 * cp.sum_squares(A @ x_model - b) + lam * cp.norm1(x_model)
prob_lasso = cp.Problem(cp.Minimize(cout))

# Le ruleset DCP certifie la convexite avant tout calcul.
print('Le probleme est-il conforme au ruleset DCP ? ', prob_lasso.is_dcp())

# Garde-fou pedagogique : inverser le signe du terme L1 rend l'objectif concave.
prob_concave = cp.Problem(
    cp.Minimize(0.5 * cp.sum_squares(A @ x_model - b) - lam * cp.norm1(x_model)))
print('Avec -lam*norm1 (objectif concave) : DCP =', prob_concave.is_dcp(),
      '-> refuse avant toute evaluation')

# Resolution par un solveur nomme
t0 = time.time()
prob_lasso.solve(solver=cp.CLARABEL)
t_clarabel = time.time() - t0
x_cvx = x_model.value

print()
print(f'CLARABEL : status = {prob_lasso.status}')
print(f'  objectif     = {objectif(x_cvx):.8f}')
print(f'  |x|_0 (1e-6) = {int(np.sum(np.abs(x_cvx) > 1e-6))}   (verite : k={k})')
print(f'  iterations   = {prob_lasso.solver_stats.num_iters}')
print(f'  temps (mur)  = {t_clarabel:.3f}s')

Le probleme est-il conforme au ruleset DCP ?  True
Avec -lam*norm1 (objectif concave) : DCP = False -> refuse avant toute evaluation

CLARABEL : status = optimal
  objectif     = 0.36064507
  |x|_0 (1e-6) = 93   (verite : k=30)
  iterations   = 12
  temps (mur)  = 0.689s

Lecture — un modèle déclaratif, et une solution au sixième chiffre

is_dcp() rend True sur la formulation correcte et False sur la variante concave : la couche de modélisation ne résout pas seulement, elle vérifie. C’est une différence de nature avec le from scratch — là-bas, une erreur de signe se paie en itérations silencieuses ; ici, elle se paie à la compilation.

Le solveur prend ce modèle et rend une solution dont l’objectif s’affiche à la huitième décimale, en 12 itérations de méthode de point intérieur. Ce chiffre est à lire avec le \(\S5\) : la coordinate descent de 2.11c converge en 82 balayages, mais chaque balayage est une boucle Python sur 500 coordonnées. Un point intérieur fait peu d’itérations très coûteuses, une méthode de premier ordre fait beaucoup d’itérations bon marché — c’est exactement ce que le \(\S4\) met en tableau.

La parcimonie mesurée est de 93 coefficients non nuls (au seuil \(10^{-6}\)) pour \(k = 30\) vrais : on est au \(\lambda\) retenu de 2.11b, celui qui sur-parcimonie volontairement le support — le même régime que le « 92 » de 2.11c. Les deux nombres ne diffèrent que par la position de la dernière coordonnée marginale autour du seuil, pas par le point de fonctionnement.

3. Ce que le solveur voit : la reformulation canonique

Le modèle déclaré compte 500 inconnues. Le solveur, lui, ne reçoit pas ce problème tel quel : cvxpy le canonicalise d’abord, c’est-à-dire qu’il le réécrit dans la forme que tous les solveurs coniques savent lire :

\[\min_x \; \tfrac12 x^{\top}Px + c^{\top}x \quad \text{s.c.} \quad Ax + s = b, \; s \in K\]

Deux reformulations mécaniques ont lieu, et get_problem_data les rend mesurables plutôt que mystérieuses :

• \(\|x\|_1 = \sum_j t_j\) sous les \(2p\) contraintes \(-t_j \le x_j \le t_j\) — l’épigraphe de la norme, qui transforme une valeur absolue non différentiable en un cône ;
• le terme quadratique \(\tfrac12\|Ax - b\|^2\) est porté par une matrice \(P\), le résidu devenant une variable à part entière reliée par une égalité.

Le coût de cette traduction est le sujet de cette section : le problème transmis n’a pas la taille du problème écrit.

data_conique, chaine_canonique, inverse = prob_lasso.get_problem_data(cp.CLARABEL)
Ac, Pc, cc = data_conique['A'], data_conique['P'], data_conique['c']

print('Ce que CLARABEL recoit reellement :')
print(f'  variables canoniques : {Ac.shape[1]}   (le modele en declare {p})')
print(f'  lignes canoniques    : {Ac.shape[0]}')
print(f'  contraintes          : {data_conique["dims"]}')
print(f'  A canonique          : {Ac.shape}, {Ac.nnz} non-nuls '
      f'({100 * Ac.nnz / (Ac.shape[0] * Ac.shape[1]):.1f}% de densite, stockage creux)')
print(f'  P (partie quadratique): {Pc.shape}, {Pc.nnz} non-nuls')
cnz = np.asarray(cc)[np.asarray(cc) != 0]
print(f'  cout lineaire c      : {cnz.size} composantes non nulles, '
      f'toutes egales a {np.unique(np.round(cnz, 8))[0]:.8f} (= lambda)')

Ce que CLARABEL recoit reellement :
  variables canoniques : 1200   (le modele en declare 500)
  lignes canoniques    : 1200
  contraintes          : 200 equalities, 1000 inequalities, 0 exponential cones, 
SOC constraints: [], PSD constraints: [],
 3d power cones [], [].
  A canonique          : (1200, 1200), 102200 non-nuls (7.1% de densite, stockage creux)
  P (partie quadratique): (1200, 1200), 200 non-nuls
  cout lineaire c      : 500 composantes non nulles, toutes egales a 0.01246456 (= lambda)

Lecture — 500 inconnues deviennent 1200 variables

Le problème transmis à CLARABEL porte 1200 variables pour 500 inconnues déclarées, et 1200 lignes de contraintes réparties en 200 égalités (les résidus \(r = Ax - b\) du terme quadratique) et 1000 inégalités (\(2 \times 500\) : l’épigraphe de la norme \(\ell_1\), une paire \(-t_j \le x_j \le t_j\) par coordonnée).

Le vecteur de coût linéaire le confirme sans ambiguïté : ses 500 composantes non nulles valent toutes exactement \(\lambda\) — ce sont les variables d’épigraphe \(t_j\), celles que la reformulation a introduites, et que le modélisateur n’a jamais écrites. La matrice \(P\) est extrêmement creuse (200 non-nuls : le terme quadratique ne touche que les résidus), et \(A\) l’est à 7 %.

C’est le contrat du déclaratif, énoncé honnêtement : on paie en structure canonique ce qu’on économise en lignes de code. Un modèle déclaratif produit un problème 2,4 fois plus gros que le problème mathématique. Le from scratch n’a pas ce coût de traduction — mais il n’a pas non plus la garantie de convexité, ni le choix du solveur. C’est cet arbitrage que le \(\S9\) chiffre à l’échelle.

4. Un problème, trois solveurs : le choix est un paramètre

Le même modèle, résolu par trois solveurs qui appartiennent à trois familles algorithmiques distinctes :

Solveur	Famille	Ce qu’il fait du problème canonique
CLARABEL	point intérieur primal-dual	peu d’itérations, algèbre linéaire dense, haute précision
SCS	premier ordre (splitting d’opérateurs)	beaucoup d’itérations bon marché, précision réglable
OSQP	ADMM (quadratique)	itérations très bon marché sur la forme QP

cvxpy accepte les trois sur le même objet Problem : changer de solveur est un argument de fonction, pas une réécriture. C’est la démonstration la plus directe de la séparation modeling / solving annoncée au \(\S2\).

parametres = {
    'CLARABEL': {},
    'SCS': {'eps': 1e-9, 'max_iters': 200000},
    'OSQP': {'eps_abs': 1e-9, 'eps_rel': 1e-9, 'max_iter': 200000},
}

print(f'{"solveur":9s} {"status":8s} {"objectif":>12s} {"|x|_0":>6s} {"iterations":>10s} '
      f'{"temps":>8s} {"|x-x_CLARABEL|_2":>16s}')
resultats = {}
for nom, kw in parametres.items():
    v = cp.Variable(p)
    pr = cp.Problem(cp.Minimize(0.5 * cp.sum_squares(A @ v - b) + lam * cp.norm1(v)))
    t0 = time.time()
    pr.solve(solver=nom, **kw)
    t = time.time() - t0
    xx = v.value
    resultats[nom] = (pr.status, objectif(xx), int(np.sum(np.abs(xx) > 1e-6)),
                      pr.solver_stats.num_iters, t, np.linalg.norm(xx - x_cvx))
    st, ob, nn, it, tt, ec = resultats[nom]
    print(f'{nom:9s} {st:8s} {ob:12.8f} {nn:6d} {it:10d} {tt:7.3f}s {ec:16.2e}')

solveur   status       objectif  |x|_0 iterations    temps |x-x_CLARABEL|_2
CLARABEL  optimal    0.36064507     93         12   0.828s         0.00e+00
SCS       optimal    0.36064506     91        450   0.325s         4.63e-05
OSQP      optimal    0.36064506     91       2775   1.117s         4.63e-05

Lecture — même optimum, trois tarifs, trois récits d’itération

Les trois solveurs rendent le même statut (optimal) et le même objectif — les écarts se lisent à la huitième décimale, du même ordre que la tolérance demandée. Les solutions coïncident à \(\sim 5\cdot10^{-5}\) en norme : comme en 2.11c \(\S4\), la surface du Lasso est plate le long des directions faiblement corrélées, et trois solveurs qui visent le même optimum s’arrêtent à des points équivalents.

Ce qui les sépare vraiment, c’est le récit d’itération, et il est structurel :

• CLARABEL fait une douzaine de pas de Newton sur le système d’optimalité — chaque pas résout un système linéaire dense, mais il n’en faut qu’une poignée ;
• SCS fait plusieurs centaines d’itérations d’opérateurs proximaux : chaque itération est une multiplication matrice-vecteur, très bon marché ;
• OSQP en fait plusieurs milliers, chacune encore moins chère, sur la forme quadratique du problème.

Aucun n’est « meilleur » dans l’absolu : le point intérieur gagne sur la précision par défaut, les méthodes de premier ordre gagnent quand le problème devient trop gros pour une factorisation dense — c’est exactement la ligne de partage que le \(\S9\) mesure sur \(p\) croissant.

Note d’honnêteté sur la parcimonie : |x|_0 varie de quelques unités d’un solveur à l’autre. Ce n’est pas une divergence de solution mais un effet de seuil : les coordonnées marginales valent \(\sim 10^{-6}\), et les compter ou non dépend du seuil de comptage, pas du solveur. Le critère qui, lui, ne triche pas est l’écart d’objectif.

5. La confrontation : cvxpy contre from scratch et sklearn, même \(\lambda\)

2.11c a confronté la coordinate descent from scratch, ISTA et sklearn sur ce problème. On ajoute ici le quatrième moteur — le solveur conique — sur le même \(\lambda\), avec les mêmes définitions d’objectif et de parcimonie. La question posée est celle de l’acceptance du bloc B.7 : les solveurs from scratch convergent-ils vers le même optimum global ?

Les trois from scratch sont repris tels quels de leurs notebooks d’origine, sans retouche : la coordinate descent de 2.11c \(\S2\) (boucle GSM, certificat KKT) et l’ISTA compacte de 2.11b (cellule 5), telle que reproduite en 2.11c \(\S4\).

def lasso_cd(A, b, lam, x0=None, max_iter=2000, tol=1e-6):
    """Coordinate descent (GSM) — reprise de 2.11c S2, meme boucle, meme critere d'arret."""
    n_, p_ = A.shape
    x = np.zeros(p_) if x0 is None else x0.copy()
    r = b - A @ x
    Asq = np.einsum('ij,ij->j', A, A)
    obj_local = lambda xx: 0.5 * np.linalg.norm(A @ xx - b) ** 2 + lam * np.abs(xx).sum()
    histo = []
    for _ in range(max_iter):
        for j in range(p_):
            r += A[:, j] * x[j]
            z = A[:, j] @ r
            x[j] = np.sign(z) * max(abs(z) - lam, 0.0) / max(Asq[j], 1e-12)
            r -= A[:, j] * x[j]
        histo.append(obj_local(x))
        if len(histo) > 1 and abs(histo[-1] - histo[-2]) <= tol * (1 + abs(histo[-1])):
            break
    return x, len(histo)


def prox_l1(x, t):
    return np.sign(x) * np.maximum(np.abs(x) - t, 0.0)


def ista_compact(A, b, lam, n_iter=2000):
    """ISTA — reprise de 2.11b cellule 5 / 2.11c S4 : gradient + soft-thresholding."""
    L = np.linalg.norm(A, 2) ** 2
    x = np.zeros(A.shape[1])
    for _ in range(n_iter):
        x = prox_l1(x - A.T @ (A @ x - b) / L, lam / L)
    return x


from sklearn.linear_model import Lasso

t0 = time.time()
x_cd, it_cd = lasso_cd(A, b, lam)
t_cd = time.time() - t0

t0 = time.time()
x_ista = ista_compact(A, b, lam)
t_ista = time.time() - t0

t0 = time.time()
skl = Lasso(alpha=lam / n, fit_intercept=False, max_iter=2000, tol=1e-8)  # alpha = lam/n (cf 2.11c S3)
skl.fit(A, b)
t_skl = time.time() - t0
x_skl = skl.coef_.copy()

moteurs = [('cvxpy/CLARABEL', x_cvx, t_clarabel, f'{prob_lasso.solver_stats.num_iters} iter.'),
           ('CD from scratch (2.11c)', x_cd, t_cd, f'{it_cd} epochs'),
           ('ISTA (2.11b)', x_ista, t_ista, '2000 iter.'),
           ('sklearn CD (2.11c)', x_skl, t_skl, f'{skl.n_iter_} iter.')]

print(f'{"moteur":26s} {"objectif":>12s} {"|x|_0":>6s} {"temps":>8s} {"|x-x_cvxpy|_2":>14s}  detail')
for nom, xx, tt, det in moteurs:
    print(f'{nom:26s} {objectif(xx):12.8f} {int(np.sum(np.abs(xx) > 1e-6)):6d} '
          f'{tt:7.3f}s {np.linalg.norm(xx - x_cvx):14.2e}  {det}')

x_cd_ref = resultats['CLARABEL'][1]
print()
print('Ecarts d objectif au solveur conique (le critere qui ne depend pas du seuil) :')
for nom, xx, tt, det in moteurs[1:]:
    print(f'  {nom:26s} : {objectif(xx) - x_cd_ref:+.3e}')

moteur                         objectif  |x|_0    temps  |x-x_cvxpy|_2  detail
cvxpy/CLARABEL               0.36064507     93   0.689s       0.00e+00  12 iter.
CD from scratch (2.11c)      0.36064623     92   0.519s       3.89e-03  82 epochs
ISTA (2.11b)                 0.36064506     91   0.146s       4.63e-05  2000 iter.
sklearn CD (2.11c)           0.36064506     91   0.010s       4.63e-05  116 iter.

Ecarts d objectif au solveur conique (le critere qui ne depend pas du seuil) :
  CD from scratch (2.11c)    : +1.167e-06
  ISTA (2.11b)               : -2.789e-09
  sklearn CD (2.11c)         : -2.789e-09

Lecture — quatre moteurs, un seul optimum global

Le verdict est net et c’est la thèse de l’acceptance : les quatre moteurs atteignent le même objectif. Les écarts entre objectifs se comptent en \(10^{-6}\) ou moins, alors que l’objectif lui-même vaut \(0{,}361\) — cinq ordres de grandeur en dessous. Pour un problème dont la solution n’est pas unique (surface plate), c’est le seul critère de comparaison qui ait un sens : les solutions elles-mêmes diffèrent de quelques \(10^{-3}\) en norme sans que cela change quoi que ce soit au coût.

Critère	cvxpy/CLARABEL	CD from scratch	ISTA	sklearn
Objectif (\(\lambda = 0{,}0125\))	0,360645	0,360646	0,360645	0,360645
Statut	optimal (certifié DCP)	certificat KKT	point fixe proximal	optimal
Itérations	~12	~82 epochs	2000 (fixes)	warm start

Ce que chaque colonne apporte, et qu’aucune autre ne remplace : le from scratch explique la mécanique (soft-threshold par coordonnée, point fixe proximal) ; sklearn optimise la boucle ; cvxpy change le problème qu’on peut écrire, pas seulement la vitesse à laquelle on le résout. Un group Lasso, une contrainte \(\|x\|_1 \le \tau\), un objectif à trois termes : le from scratch demanderait trois algorithmes, le déclaratif se contente de l’énoncé — c’est le sujet des exercices du \(\S11\).

6. La dualité sort du solveur : le certificat KKT sans le recoder

En 2.11c \(\S2\), le certificat d’optimalité était codé à la main : après convergence, on vérifiait que \(|A_{\cdot j}^{\top} r| \le \lambda\) hors du support. Avec cvxpy, ce certificat est un produit du solveur : il suffit d’écrire la forme épigraphe

\[\min_{x,t} \; \tfrac12\|Ax - b\|^2 + \lambda \sum_j t_j \quad \text{s.c.} \quad -t_j \le x_j \le t_j\]

et de lire constraint.dual_value. Les multiplicateurs \(\mu^{+}, \mu^{-} \ge 0\) des deux familles de contraintes portent toute l’information duale :

• stationnarité : la corrélation résiduelle s’écrit \(c = A^{\top}(b - Ax) = \mu^{+} - \mu^{-}\) ;
• complémentarité : \(\mu^{+}_j (t_j - x_j) = 0\) et \(\mu^{-}_j (t_j + x_j) = 0\) ;
• d’où la condition KKT du Lasso, retrouvée sans l’écrire : \(|c_j| \le \lambda\), avec égalité exactement sur le support.

Autrement dit : la théorie de 2.11c ne se réimplémente pas ici, elle se lit.

# Forme epigraphe explicite : les multiplicateurs deviennent accessibles.
x_epi = cp.Variable(p)
t_epi = cp.Variable(p, nonneg=True)
c_pos = x_epi <= t_epi
c_neg = -x_epi <= t_epi
prob_epi = cp.Problem(
    cp.Minimize(0.5 * cp.sum_squares(A @ x_epi - b) + lam * cp.sum(t_epi)),
    [c_pos, c_neg])
prob_epi.solve(solver=cp.CLARABEL)

x_epi_val = x_epi.value
t_epi_val = t_epi.value
mu_pos = np.asarray(c_pos.dual_value).ravel()
mu_neg = np.asarray(c_neg.dual_value).ravel()
corr = A.T @ (b - A @ x_epi_val)

print(f'Epigraphe  : status = {prob_epi.status}, objectif = {objectif(x_epi_val):.8f}')
print(f'  ecart a la forme naturelle : |x_epi - x_cvx|_2 = {np.linalg.norm(x_epi_val - x_cvx):.2e}')
print()
print(f'Stationnarite   : max |c - (mu+ - mu-)|        = {np.max(np.abs(corr - (mu_pos - mu_neg))):.2e}')
print(f'Complementarite : max |mu+ * (t - x)|          = {np.max(np.abs(mu_pos * (t_epi_val - x_epi_val))):.2e}')
print(f'Complementarite : max |mu- * (t + x)|          = {np.max(np.abs(mu_neg * (t_epi_val + x_epi_val))):.2e}')

nonz = np.abs(x_epi_val) > 1e-6
print()
print(f'max |c_j| global          = {np.max(np.abs(corr)):.6f}   vs lambda = {lam:.6f}')
print(f'sur le support (|x_j|>0)  : max | |c_j| - lambda | = {np.max(np.abs(np.abs(corr[nonz]) - lam)):.2e}'
      f'   (egalite des deux cotes)')
print(f'hors du support (x_j=0)   : max |c_j| = {np.max(np.abs(corr[~nonz])):.6f}'
      f'   (marge {lam - np.max(np.abs(corr[~nonz])):.2e} sous le seuil)')

Epigraphe  : status = optimal, objectif = 0.36064506
  ecart a la forme naturelle : |x_epi - x_cvx|_2 = 3.55e-05

Stationnarite   : max |c - (mu+ - mu-)|        = 1.20e-15
Complementarite : max |mu+ * (t - x)|          = 6.94e-11
Complementarite : max |mu- * (t + x)|          = 4.88e-11

max |c_j| global          = 0.012465   vs lambda = 0.012465
sur le support (|x_j|>0)  : max | |c_j| - lambda | = 2.48e-05   (egalite des deux cotes)
hors du support (x_j=0)   : max |c_j| = 0.012379   (marge 8.56e-05 sous le seuil)

Lecture — la théorie lue, pas recodée

Les trois lectures se tiennent aux tolérances du solveur : la stationnarité est satisfaite à \(10^{-15}\) (c’est une égalité linéaire, un point intérieur la tient à la précision machine), et les deux complémentarités à \(10^{-10}\). Le certificat KKT de 2.11c se retrouve intact, mais il n’a pas été écrit : il est sorti des multiplicateurs.

La nuance la plus instructive est sur le support. Sur les coordonnées non nulles, \(|c_j| = \lambda\) exactement (à \(2{,}5\cdot10^{-5}\), la précision réelle du point intérieur) ; hors du support, \(|c_j|\) reste sous \(\lambda\) avec une marge de \(\sim 10^{-4}\). Autrement dit la signature du Lasso — corrélation exactement au seuil sur le support, strictement en dessous ailleurs — est visible dans la sortie du solveur, alors qu’elle demandait en 2.11c une boucle de vérification après coup.

Réserve d’honnêteté, utile pour les exercices : \(\mu^{+} + \mu^{-} = \lambda\) n’est pas vérifié point par point quand \(x_j \approx 0\). À ces coordonnées, \(t_j = 0\) et les deux contraintes sont actives simultanément : le couple \((\mu^{+}_j, \mu^{-}_j)\) n’est alors pas unique (toute répartition de somme \(\lambda\) convient), et le solveur en choisit une. C’est une dégénérescence du dual, pas une erreur — le certificat qui compte, \(c = \mu^{+} - \mu^{-}\), reste exact.

7. Un modèle, dix-huit résolutions : le chemin de régularisation et le seuil \(\lambda_{\max}\)

La force du déclaratif se voit mieux quand le problème change sans que l’algorithme change. Deux exercices de style dans cette section, sur le même modèle :

1. le chemin : résoudre pour une grille de \(\lambda\) décroissants — c’est ce que fait LassoCV en interne (2.11c \(\S3\)), et c’est aussi la première marche des solveurs LARS ;
2. le seuil : \(\lambda_{\max} = \|A^{\top}b\|_{\infty}\) annule la solution. Le from scratch le subit (il converge vers zéro) ; le déclaratif le calcule à partir de la condition KKT \(|c_j| \le \lambda\), et le solveur n’est là que pour confirmer.

Le seuil mérite qu’on s’y arrête, parce que c’est un point de dégénérescence : à \(\lambda = \lambda_{\max}\) exactement, la solution n’est pas unique — \(x = 0\) est optimal, et quelques coordonnées peuvent l’être aussi. La bonne mesure n’est donc pas le nombre de non-nuls, mais l’écart d’objectif à zéro.

print('--- Chemin de regularisation : 12 resolutions du meme modele ---')
print(f'{"lambda":>9s} {"|x|_0":>6s} {"objectif":>11s} {"temps":>8s}')
for l in np.geomspace(lam * 20, lam / 20, 12):
    v = cp.Variable(p)
    pr = cp.Problem(cp.Minimize(0.5 * cp.sum_squares(A @ v - b) + l * cp.norm1(v)))
    t0 = time.time()
    pr.solve(solver=cp.CLARABEL)
    t = time.time() - t0
    xx = v.value
    print(f'{l:9.5f} {int(np.sum(np.abs(xx) > 1e-6)):6d} {objectif(xx, l):11.6f} {t:7.3f}s')

print()
print('--- Le seuil lambda_max, verifie par le solveur ---')
print(f'lambda_max calcule (max_j |A_j^T b|) = {lam_max:.6f}   |   f(0) = 0.5||b||^2 = '
      f'{0.5 * np.linalg.norm(b) ** 2:.6f}')
f_zero = 0.5 * np.linalg.norm(b) ** 2
for fac in [0.999, 1.0, 1.001]:
    l = fac * lam_max
    v = cp.Variable(p)
    pr = cp.Problem(cp.Minimize(0.5 * cp.sum_squares(A @ v - b) + l * cp.norm1(v)))
    pr.solve(solver=cp.CLARABEL)
    xx = v.value
    print(f'lambda = {fac:.3f} * lambda_max : |x|_0(1e-6) = {int(np.sum(np.abs(xx) > 1e-6)):2d}'
          f' | max|x| = {np.abs(xx).max():.2e} | f(x) - f(0) = {objectif(xx, l) - f_zero:+.2e}')

--- Chemin de regularisation : 12 resolutions du meme modele ---
   lambda  |x|_0    objectif    temps
  0.24929     30    6.234041   0.588s
  0.14460     35    3.824345   0.563s
  0.08387     39    2.297614   0.558s
  0.04865     38    1.361812   0.569s
  0.02822     48    0.801935   0.620s
  0.01637     80    0.470996   0.600s
  0.00949    110    0.275978   0.661s
  0.00551    137    0.161340   0.807s
  0.00319    165    0.094123   0.799s
  0.00185    179    0.054802   0.684s
  0.00107    186    0.031864   0.798s
  0.00062    192    0.018509   0.677s

--- Le seuil lambda_max, verifie par le solveur ---
lambda_max calcule (max_j |A_j^T b|) = 3.191173   |   f(0) = 0.5||b||^2 = 24.915170
lambda = 0.999 * lambda_max : |x|_0(1e-6) =  1 | max|x| = 3.12e-03 | f(x) - f(0) = -4.98e-06
lambda = 1.000 * lambda_max : |x|_0(1e-6) =  1 | max|x| = 3.42e-04 | f(x) - f(0) = +9.94e-08
lambda = 1.001 * lambda_max : |x|_0(1e-6) =  1 | max|x| = 1.41e-06 | f(x) - f(0) = +6.39e-09

Lecture — le chemin coûte 12 résolutions, et le seuil se calcule

Le chemin montre la mécanique de la parcimonie à l’œuvre : \(\lambda\) divisé par 20 fait passer le support de quelques dizaines de coordonnées à près de deux cents, et chaque ligne est une résolution complète — c’est le prix du déclaratif, \(\sim 0{,}6\) s par point, contre la fraction de milliseconde du warm start de sklearn (lequel ne résout pas chaque point indépendamment : il chaîne les solutions, 2.11c exercice 3). Un point utile pour la culture du domaine : LassoCV de 2.11c a choisi son \(\alpha\) sur une grille ; le point de la grille de ce \(\S7\) le plus proche, \(\lambda = 0{,}2493\) (\(20\lambda\)), rend 30 non-nuls — la valeur exacte de \(k\), retrouvée par un autre chemin que la validation croisée.

Le seuil, lui, se lit proprement dans les trois lignes du tableau :

• à \(0{,}999\,\lambda_{\max}\), une coordonnée survit avec une valeur \(\sim 10^{-3}\) et améliore l’objectif d’environ \(5\cdot10^{-6}\) — la récompense existe, elle est minuscule ;
• à \(\lambda_{\max}\) exactement, l’écart à \(f(0)\) est de l’ordre du bruit du solveur (\(+10^{-7}\)) et une coordonnée traîne encore à \(3\cdot10^{-4}\) : c’est le point de dégénérescence annoncé, où plusieurs points sont optimaux ;
• au-delà (\(1{,}001\,\lambda_{\max}\)), la coordonnée rendue vaut \(1{,}4\cdot10^{-6}\) — indistinguable de zéro, l’écart d’objectif étant retombé à \(6\cdot10^{-9}\) : c’est le seuil de comptage (\(10^{-6}\)) qui la fait encore apparaître, pas la solution.

La leçon de méthode est double. D’abord, un seuil se calcule (\(\|A^{\top}b\|_{\infty}\)) plutôt que se chercher : la condition KKT le donne, le solveur ne fait que le confirmer. Ensuite, à la frontière d’un régime, le comptage de non-nuls est un mauvais instrument — c’est l’écart d’objectif qui départage une coordonnée vraie d’une coordonnée fantôme.

8. Le problème du bloc A.2 : cvxpy comme arbitre de l’ADMM

L’acceptance du bloc B.7 demande explicitement la comparaison avec ADMM et avec les méthodes proximales. Or 2.11d ne travaille pas sur le problème de 2.11c : il a son propre design, corrélé (\(\rho_{corr} = 0{,}95\), structure AR(1)), avec \(n = 100\), \(p = 50\), \(k = 5\) — un problème volontairement mal conditionné (\(\kappa(X) \approx 82\)), choisi là-bas parce qu’il allonge les trajectoires des méthodes alternées.

On rejoue donc ce problème, avec les deux définitions de 2.11d reprises telles quelles (design_sparse_corr et admm_lasso, mêmes hyperparamètres, même graine), et l’on ajoute le seul acteur qui manquait : l’arbitre conique. La question — ADMM converge-t-il vers l’optimum global ? — ne se tranche pas contre une autre méthode itérative, elle se tranche contre un optimum certifié.

def design_sparse_corr(n_=100, p_=50, k_=5, corr=0.95, snr=4.0, seed=20260914):
    """Design correle AR(1) — reprise exacte de 2.11d cellule 4."""
    local_rng = np.random.default_rng(seed)
    idx = np.arange(p_)
    Sigma = corr ** np.abs(idx[:, None] - idx[None, :])
    L = np.linalg.cholesky(Sigma)
    X = local_rng.standard_normal((n_, p_)) @ L.T
    support_ = np.sort(local_rng.choice(p_, size=k_, replace=False))
    beta_true_ = np.zeros(p_)
    beta_true_[support_] = local_rng.standard_normal(k_)
    sigma_ = np.linalg.norm(X @ beta_true_) / np.sqrt(n_) / snr
    y_ = X @ beta_true_ + local_rng.standard_normal(n_) * sigma_
    return X, beta_true_, y_, support_


def soft_threshold(z, lam_):
    return np.sign(z) * np.maximum(np.abs(z) - lam_, 0.0)


def admm_lasso(X, y, lam_, rho, max_iter=500, eps_abs=1e-4, eps_rel=1e-4):
    """ADMM pour min 0.5||X b - y||^2 + lam||z||_1 s.c. b = z — reprise de 2.11d cellule 9."""
    n_, p_ = X.shape
    XtX, Xty = X.T @ X, X.T @ y
    L = np.linalg.cholesky(XtX + rho * np.eye(p_))   # pre-factorisation du b-update
    beta = np.zeros(p_)
    z = np.zeros(p_)
    u = np.zeros(p_)
    eps_p, eps_d = np.sqrt(p_) * eps_abs, np.sqrt(p_) * eps_abs
    z_prev = z.copy()
    converged = False
    for it in range(1, max_iter + 1):
        tmp = np.linalg.solve(L, Xty + rho * (z - u))
        beta = np.linalg.solve(L.T, tmp)
        z = soft_threshold(beta + u, lam_ / rho)
        u = u + beta - z
        r_p = np.linalg.norm(beta - z)
        r_d = rho * np.linalg.norm(z - z_prev)
        tol_p = eps_p + eps_rel * max(np.linalg.norm(beta), np.linalg.norm(z))
        tol_d = eps_d + eps_rel * rho * np.linalg.norm(u)
        if r_p < tol_p and r_d < tol_d:
            converged = True
            break
        z_prev = z.copy()
    return beta, z, it, converged


X, beta_true, y, supp = design_sparse_corr()
lam_max_raw = np.max(np.abs(X.T @ y))
lam_test = 0.5 * lam_max_raw          # parametrage de 2.11d cellule 15
f_local = lambda bb: 0.5 * np.linalg.norm(X @ bb - y) ** 2 + lam_test * np.abs(bb).sum()

v = cp.Variable(X.shape[1])
pr = cp.Problem(cp.Minimize(0.5 * cp.sum_squares(X @ v - y) + lam_test * cp.norm1(v)))
t0 = time.time()
pr.solve(solver=cp.CLARABEL)
t_cvx_11d = time.time() - t0
b_opt = v.value
f_star = f_local(b_opt)

print(f'Probleme 2.11d : n={X.shape[0]}, p={X.shape[1]}, k=5, rho_corr=0.95, kappa(X)={np.linalg.cond(X):.2e}')
print(f'lambda_max_raw = {lam_max_raw:.4f} | lambda_test = 0.5*lambda_max_raw = {lam_test:.4f}')
print(f'Optimum certifie (CLARABEL) : f* = {f_star:.6f}, |b*|_0 = '
      f'{int(np.sum(np.abs(b_opt) > 1e-6))}, support vrai = {list(supp)}, temps = {t_cvx_11d:.3f}s')
print()
print(f'{"rho":>6s} {"it.":>5s} {"converge":>9s} {"f(z)":>12s} {"f(z)-f*":>11s} {"|z-b*|_2":>10s} {"temps":>8s}')
for rho in [1.0, 5.0, 20.0]:
    t0 = time.time()
    bd, zd, itd, conv = admm_lasso(X, y, lam_test, rho)
    tt = time.time() - t0
    print(f'{rho:6.1f} {itd:5d} {str(conv):>9s} {f_local(zd):12.6f} '
          f'{f_local(zd) - f_star:+11.2e} {np.linalg.norm(zd - b_opt):10.2e} {tt:7.3f}s')

Probleme 2.11d : n=100, p=50, k=5, rho_corr=0.95, kappa(X)=8.18e+01
lambda_max_raw = 115.2719 | lambda_test = 0.5*lambda_max_raw = 57.6360
Optimum certifie (CLARABEL) : f* = 158.890455, |b*|_0 = 2, support vrai = [np.int64(2), np.int64(13), np.int64(24), np.int64(44), np.int64(47)], temps = 0.017s

   rho   it.  converge         f(z)     f(z)-f*   |z-b*|_2    temps
   1.0   500     False   160.305962   +1.42e+00   1.82e-01   0.058s
   5.0   500     False   158.893498   +3.04e-03   7.01e-03   0.060s
  20.0   178      True   158.890705   +2.49e-04   2.08e-03   0.023s

Lecture — l’arbitre tranche, et le paramètre \(\rho\) réapparaît

Le tableau se lit sur trois colonnes, et chacune raconte la même histoire sous un angle différent.

La convergence : ADMM à \(\rho = 20\) franchit son critère d’arrêt en 178 itérations et tombe à \(2\cdot10^{-4}\) de l’optimum certifié. Aux petits \(\rho\), il atteint la borne de 500 itérations sans converger — c’est exactement le phénomène de 2.11d : sur un design corrélé (\(\kappa \approx 82\)), le réglage de la pénalité augmentée décide de tout. Le solveur conique, lui, n’a aucun hyperparamètre de convergence à régler : c’est son coût et son bénéfice.

L’écart à \(f^{*}\) est la mesure qui compte, et elle est doublement instructive. À \(\rho = 20\), l’écart résiduel est du même ordre que la tolérance demandée à ADMM (\(10^{-4}\) relatif) : ADMM ne s’arrête pas « près » de l’optimum par hasard, il s’arrête au niveau de précision qu’on lui a demandé. Descendre plus bas serait possible — en baissant eps_abs — au prix d’itérations supplémentaires. L’optimum conique n’est pas « meilleur », il est certifié : c’est la différence entre une convergence observée et une borne.

Le support : les deux méthodes rendent la même poignée de coordonnées, mais avec une réserve à garder en tête — à \(\lambda = 0{,}5\,\lambda_{\max}\), le problème est très parcimonieux et les coefficients vivants sont peu nombreux ; c’est le régime où un support se compare, contrairement au \(\lambda\) de 2.11c où 90 coefficients non nuls sur 500 sont majoritairement du bruit de seuil.

La conclusion de ce \(\S8\) est celle qui manquait au bloc B.7 : ADMM ne se contente pas d’être « un algorithme qui marche » — confronté à un optimum certifié par un solveur d’un autre type, il s’y ancre à la précision de son critère d’arrêt. C’est la validation croisée des deux blocs.

9. L’échelle : \(n\) fixe, \(p\) croissant

2.11c \(\S5\) a mesuré la croissance du coût pour la coordinate descent, ISTA et sklearn. On ajoute le solveur conique au tableau, sur les mêmes instances : \(n = 200\) fixe, \(k = 30\) fixe, \(p \in \{500, 2000, 4000\}\), avec le \(\lambda\) recalculé à chaque fois par la formule de Donoho-Johnstone (elle dépend de \(p\)).

Le solveur conique a une propriété que les autres n’ont pas : il ne change pas quand le problème change. Le même modèle, le même appel, seule la taille de \(A\) varie. C’est ce qui rend la mesure intéressante — le coût du déclaratif se paie en structure canonique (deux fois et demie la taille du problème, cf \(\S3\)), et cette facture croît avec \(p\).

print('n = 200 fixe, k = 30 fixe, lambda = 5 * Donoho-Johnstone(p) ; temps mesures en un passage')
print(f'{"p":>6s} {"cvxpy CLARABEL":>15s} {"sklearn CD":>12s} {"CD 50 ep.":>10s} {"ISTA 300 it.":>12s}')
for pp in [500, 2000, 4000]:
    r2 = np.random.default_rng(1000 + pp)
    A_ = r2.standard_normal((n, pp)) / np.sqrt(n)
    s_ = r2.choice(pp, size=k, replace=False)
    xt_ = np.zeros(pp)
    xt_[s_] = r2.standard_normal(k)
    b_ = A_ @ xt_ + noise_level * r2.standard_normal(n)
    lam_ = 5.0 * sigma * np.sqrt(2.0 * np.log(pp)) / np.sqrt(n)

    v = cp.Variable(pp)
    pr = cp.Problem(cp.Minimize(0.5 * cp.sum_squares(A_ @ v - b_) + lam_ * cp.norm1(v)))
    t0 = time.time()
    pr.solve(solver=cp.CLARABEL)
    t_cvx_ = time.time() - t0
    it_cvx_ = pr.solver_stats.num_iters

    t0 = time.time()
    Lasso(alpha=lam_ / n, fit_intercept=False, max_iter=5000, tol=1e-6).fit(A_, b_)
    t_skl_ = time.time() - t0

    t0 = time.time()
    lasso_cd(A_, b_, lam_, max_iter=50)
    t_cd_ = time.time() - t0

    t0 = time.time()
    L_ = np.linalg.norm(A_, 2) ** 2
    xi = np.zeros(pp)
    for _ in range(300):
        xi = prox_l1(xi - A_.T @ (A_ @ xi - b_) / L_, lam_ / L_)
    t_ista_ = time.time() - t0

    print(f'{pp:6d} {t_cvx_:12.3f}s ({it_cvx_:2d} it.) {t_skl_:11.3f}s {t_cd_:9.3f}s {t_ista_:11.3f}s')

n = 200 fixe, k = 30 fixe, lambda = 5 * Donoho-Johnstone(p) ; temps mesures en un passage
     p  cvxpy CLARABEL   sklearn CD  CD 50 ep. ISTA 300 it.
   500        0.614s (12 it.)       0.008s     0.315s       0.029s
  2000        3.687s (16 it.)       0.073s     1.134s       0.101s
  4000        8.898s (18 it.)       0.315s     2.562s       0.152s

Lecture — le déclaratif paie sa généralité à l’échelle

Le classement est sans ambiguïté, et il faut le dire tel qu’il est : sur ces instances creuses (\(n \ll p\), \(k\) très petit devant \(p\)), sklearn et ISTA écrasent le solveur conique. À \(p = 4000\), la coordinate descent de sklearn fait le travail en trois dixièmes de seconde quand CLARABEL en demande près de neuf — un facteur proche de trente, et il se creuse avec \(p\).

La raison est structurelle, pas accidentelle. Trois mécanismes se cumulent :

• la traduction : le problème transmis est 2,4 fois plus gros que le problème écrit (\(\S3\)) — l’épigraphe et les résidus se paient en variables ;
• l’ordre de l’algorithme : un point intérieur résout un système dense à chaque itération ; les méthodes de premier ordre et la coordinate descent ne font que des produits matrice-vecteur, à \(O(np)\) l’itération ;
• l’adaptation au problème : sklearn est un algorithme spécialisé Lasso, avec warm start et règles d’activation. Le solveur conique est générique — il résout le même Lasso que le group Lasso que le problème à contraintes, sans rien savoir de la parcimonie.

C’est le troisième point qui réconcilie ce tableau avec le \(\S5\), et c’est la réponse honnête au pourquoi du déclaratif : ce n’est pas la vitesse qui le justifie. Sur ce problème, le from scratch et sklearn gagnent. Le déclaratif gagne dès que le problème change — et le coût de ce changement, dans les autres colonnes, est une session de développement.

Tableau de synthèse du bloc B (acceptance \(\S8\) de #16061) :

Axe	From scratch (2.11b/2.11c/2.11d)	sklearn (2.11c)	cvxpy (2.11e)
Objectif final (même \(\lambda\))	identique	identique	identique
Temps (\(p=500\))	~0,1-0,5 s	~0,01 s	~0,6 s
Temps (\(p=4000\))	~0,2-2 s	~0,3 s	plusieurs s
Itérations	82 epochs / 2000 / 500	warm start	~12-18
Coût du changement de problème	un nouvel algorithme	non couvert	réécrire l’énoncé
Garantie de convexité	analyse manuelle	non	ruleset DCP
Certificat dual	codé à la main	non exposé	lu dans les duales

10. Synthèse

1. Le déclaratif change la nature du code. On écrit le problème, pas l’algorithme ; le solveur devient un paramètre. La contrepartie immédiate est le coût de traduction en forme canonique : 500 inconnues déclarées, 1200 variables transmises, 1000 inégalités d’épigraphe pour une seule norme \(\ell_1\).
2. Le même problème, quatre moteurs, un optimum. cvxpy, la coordinate descent from scratch, ISTA et sklearn atteignent le même objectif à \(10^{-6}\) près. C’est le verdict que l’acceptance du bloc B.7 demandait ; pour un problème à solution non unique, l’objectif est le bon critère, pas la distance entre solutions.
3. La dualité se lit. Le certificat KKT de 2.11c — \(|c_j| \le \lambda\), égalité sur le support — sort directement des multiplicateurs de la forme épigraphe. La théorie n’est pas recodée, elle est exploitée.
4. Le seuil \(\lambda_{\max} = \|A^{\top}b\|_{\infty}\) se calcule. Il se déduit de la condition KKT, et la dégénérescence au seuil (\(x = 0\) optimal sans être unique) interdit de compter la parcimonie sans regarder l’écart d’objectif.
5. ADMM s’ancre à l’optimum certifié sur son propre problème (2.11d, design corrélé), à la précision de son critère d’arrêt — la validation croisée des deux blocs.
6. Le déclaratif ne gagne pas la course de vitesse, il gagne celle du changement. Sur ce Lasso creux, sklearn reste un ordre de grandeur plus rapide. La valeur de cvxpy est ailleurs : écrire un group Lasso, une contrainte \(\|x\|_1 \le \tau\) ou un objectif à trois termes sans écrire un nouvel algorithme — et sans se tromper de convexité, puisque le ruleset DCP refuse la formulation avant de la résoudre.

11. Exercices

Trois exercices de difficulté croissante, conformes à la convention C.1 (aucune erreur volontaire : les stubs s’exécutent et renvoient None). Tous portent sur le même problème que le notebook — le contexte, les données et les objets (A, b, p, lam, lam_max) sont déjà en mémoire.

Exercice 1 — group Lasso : le déclaratif absorbe la structure

Le group Lasso pénalise des groupes de coordonnées entiers :

\[\min_x \; \tfrac12 \|Ax - b\|_2^2 + \lambda \sum_{g=1}^{G} \|x_g\|_2\]

où les \(G\) groupes forment une partition de \(\{1, \dots, p\}\). La différence avec le Lasso est structurelle : le prox de \(\|\cdot\|_2\) est une contraction par bloc (pas un soft-threshold coordonnée à coordonnée), et un groupe entier est annulé ou préservé — aucune coordonnée isolée ne survit seule dans un groupe nul.

Ce que l’exercice doit mettre en évidence : en déclaratif, la pénalité est une autre expression (cp.norm(x_g, 2)), pas un autre algorithme. En from scratch, le prox change, la boucle change, le certificat change.

# Indice : partitionner les p coordonnees en G groupes contigus de taille p//G
# Indice : la penalite s'ecrit cp.sum([cp.norm(x[groupe], 2) for groupe in groupes])
# Indice : compter les groupes ANNULES, pas les coordonnees (c'est la bonne granularite ici)

def exercice_1_group_lasso(A, b, lam, G=50, solver=cp.CLARABEL):
    """Exercice 1 : group Lasso declaratif — min 0.5||Ax-b||^2 + lam*sum_g ||x_g||_2.

    Retourne (x, nb_groupes_annules, nb_coordonnees_non_nulles), ou None si non complete.
    """
    p_ = A.shape[1]
    x = cp.Variable(p_)
    # TODO etudiant : construire la partition en G groupes contigus de taille p//G
    # TODO etudiant : ecrire la penalite somme des normes 2 par groupe
    # TODO etudiant : resoudre, puis compter les groupes dont la norme est sous 1e-6
    return None  # TODO etudiant

# Etape 1 : comparer au Lasso sur le meme lambda — le group Lasso annule-t-il par blocs ?
x_gr, n_gr, n_co = (None, None, None)
print('Exercice 1 a completer en TP.')

Exercice 1 a completer en TP.

Exercice 2 — la forme contrainte, et le lien dual avec la forme pénalisée

Le Lasso a deux écritures équivalentes, reliées par la dualité lagrangienne :

\[\min_x \tfrac12\|Ax - b\|_2^2 + \lambda\|x\|_1 \quad \Longleftrightarrow \quad \min_x \tfrac12\|Ax - b\|_2^2 \;\text{ s.c. } \|x\|_1 \le \tau\]

Pour tout \(\lambda > 0\) il existe un \(\tau(\lambda)\) tel que les deux problèmes aient la même solution — et \(\tau(\lambda) = \|x^{\star}(\lambda)\|_1\), c’est-à-dire la norme \(\ell_1\) de la solution pénalisée elle-même. C’est le contenu du théorème de Lagrange appliqué à la contrainte \(\|x\|_1 \le \tau\) : le multiplicateur dual de cette contrainte est \(\lambda\).

# Indice : le multiplicateur dual de la contrainte s'obtient par contrainte.dual_value
# Indice : balayer tau autour de ||x_cvx||_1 et verifier que la solution ne bouge plus
# Indice : comparer le dual_value de la contrainte a lam (c'est le sens du theoreme)

def exercice_2_forme_contrainte(A, b, lam, x_star, tol=1e-3):
    """Exercice 2 : forme contrainte min 0.5||Ax-b||^2 s.c. ||x||_1 <= tau.

    A partir de x_star (solution penalisee), retrouver tau = ||x_star||_1, resoudre la
    forme contrainte, comparer les solutions et lire le multiplicateur dual.
    Retourne un dict ou None si non complete.
    """
    p_ = A.shape[1]
    x = cp.Variable(p_)
    tau = np.abs(x_star).sum() if x_star is not None else None
    # TODO etudiant : ecrire la contrainte cp.norm1(x) <= tau
    # TODO etudiant : resoudre et mesurer |x_contraint - x_star|_2
    # TODO etudiant : lire contrainte.dual_value et le comparer a lam
    return None  # TODO etudiant

res2 = exercice_2_forme_contrainte(A, b, lam, x_cvx)
print('Exercice 2 a completer en TP.')

Exercice 2 a completer en TP.

Exercice 3 — \(\lambda_{\max}\) et le début du chemin, par le solveur seul

Le \(\S7\) a donné \(\lambda_{\max}\) par le calcul (\(\|A^{\top}b\|_{\infty}\)) et vérifié le seuil à trois points. L’exercice reprend la question sans la formule : trouver, par le solveur seul, le plus petit \(\lambda\) d’une grille géométrique pour lequel la solution est nulle — puis comparer à \(\|A^{\top}b\|_{\infty}\), et mesurer de combien la recherche par résolutions successives est plus chère que le calcul direct.

C’est la morale du \(\S7\) retournée : le solveur confirme un seuil, il ne le découvre pas. Savoir laquelle des deux opérations on est en train de faire est ce qui distingue l’usage éclairé du déclaratif de l’usage magique.

# Indice : grille np.geomspace(1e-2, 10.0, 25) ; pour chaque lambda, resoudre et tester max|x| < 1e-6
# Indice : comparer le plus petit lambda annulant de la grille a lam_max (lequel est au-dessus ?)
# Indice : chronometrer la boucle et comparer a une seule evaluation de lam_max

def exercice_3_lambda_max_par_solveur(A, b, n_points=25):
    """Exercice 3 : trouver par le solveur le plus petit lambda d'une grille annulant x.

    Retourne (lambda_trouve, lambda_max_analytique, temps_boucle) ou None si non complete.
    """
    p_ = A.shape[1]
    grille = np.geomspace(1e-2, 10.0, n_points)
    # TODO etudiant : boucler sur la grille, resoudre, tester max|x| < 1e-6
    # TODO etudiant : comparer a np.max(np.abs(A.T @ b)) et au cout d'une seule evaluation
    return None  # TODO etudiant

res3 = exercice_3_lambda_max_par_solveur(A, b)
print('Exercice 3 a completer en TP.')

Exercice 3 a completer en TP.

Transition : la comparaison des solveurs

La chaîne d’optimisation convexe se referme en confrontant les implémentations from scratch de la série (2.11, 2.11b, 2.11d) aux solveurs SOTA appelés via CVXPY, sur les mêmes problèmes et avec les mêmes témoins.

Suivant : 2.11f-Comparaison-Optimisation-Convexe

12. Références

• Diamond, S., Boyd, S. (2016). CVXPY: A Python-Embedded Modeling Language for Convex Optimization. JMLR 17(83). — la séparation modeling / solving, le ruleset DCP, la canonicalisation.
• Grant, M., Boyd, S., Ye, Y. (2006). Disciplined Convex Programming. In Global Optimization: From Theory to Implementation. — la certification de convexité à la compilation.
• Boyd, S., Parikh, N., Chu, E., Peleato, B., Eckstein, J. (2011). Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers. Foundations and Trends in ML 3(1). — §6.4 (Lasso) : le problème du \(\S8\).
• Goulart, P., Chen, Y. (2024). Clarabel: An interior-point solver for conic programs with quadratic objectives. — CLARABEL, le solveur d’intérieur-point utilisé par défaut ici (successeur moderne d’ECOS, qui n’est plus distribué avec les versions récentes de cvxpy).
• O’Donoghue, B., Chu, E., Parikh, N., Boyd, S. (2016). Conic Optimization via Operator Splitting and Homogeneous Self-Dual Embedding. JOTA 169(3). — SCS, le solveur de premier ordre du \(\S4\).
• Stellato, B., Banjac, G., Goulart, P., Bemporad, A., Boyd, S. (2020). OSQP: An Operator Splitting Solver for Quadratic Programs. Mathematical Programming Computation 12. — OSQP, le solveur QP.
• Friedman, J., Hastie, T., Tibshirani, R. (2010). Regularization Paths for Generalized Linear Models via Coordinate Descent. JSS 33(1). — la coordinate descent de sklearn, comparée au \(\S5\).
• Yuan, M., Lin, Y. (2006). Model Selection and Estimation in Regression with Grouped Variables. JRSS-B 68(1). — le group Lasso de l’exercice 1.
• Notebooks frères : 2.11b-Proximal-Operators-From-Scratch.ipynb (bloc A.3) · 2.11c-Lasso-SOTA-Comparison.ipynb (bloc B.6) · 2.11d-Optimisation-ADMM-From-Scratch.ipynb (bloc A.2) — même problème, mêmes seeds pour 2.11b/2.11c.
• Documentation cvxpy : Problem.get_problem_data, SolverStats, Constraint.dual_value — les trois API utilisées aux \(\S3\), \(\S4\) et \(\S6\).