2.11c — Lasso SOTA : coordinate descent, LassoCV, et quand le SOTA prend le relais

Serie 02-ML-Cours (Data Science With Agents) — Bloc B.6 de l’issue #16061 (Optimisation convexe avancee from scratch). Le Bloc A (2.11b) a construit ISTA/FISTA from scratch et les a confrontes a sklearn.Lasso et cvxpy. Ce notebook prend le cote SOTA : il repond a la question symetrique — que fait exactement sklearn quand il resout un Lasso, et quand vaut-il mieux le laisser faire ?

• la coordinate descent (le moteur de sklearn.linear_model.Lasso) demontree from scratch, avec sa solution fermee par coordonnee et son critere d’arret ;
• le choix de lambda par validation croisee (LassoCV) : le lambda theorique de Donoho-Johnstone est un a priori ; la CV le remplace par une mesure sur les donnees ;
• la comparaison honnete ISTA (from scratch) vs CD from scratch vs CD sklearn sur le meme probleme que le bloc A : objectif final, sparsite, temps, iterations ;
• la capacite a scaler (p croissant) : c’est l’axe du tableau du \(\S8\) du body d’issue — pourquoi du from scratch, quand du SOTA.

Chaine pedagogique : 2.11-Regularisation-Sparse-LASSO (la regle) -> 2.11b-Proximal-Operators-From-Scratch (ISTA/FISTA, bloc A.3) -> 2.11c (ce notebook, bloc B.6).

1. Le meme probleme que le bloc A

Pour que les nombres se comparent entre notebooks, on reprend exactement la generation de 2.11b cellules 9-10 : \(n=200\) mesures, \(p=500\) dimensions, \(k=30\) coefficients non-nuls, bruit relatif \(0{,}01\), seed 42, matrice sous-gaussienne normalisee \(A/\sqrt{n}\) (lignes de norme unitaire — ce point sera decisif au \(\S3\) pour la correspondance exacte alpha = lambda/n, et au \(\S2\) pour la forme simple de la mise a jour).

import numpy as np
import time
np.random.seed(42)

# --- Generation identique a 2.11b (bloc A.3) ---
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_clean = A @ x_true
b = b_clean + noise_level * rng.standard_normal(n)

# Lambda theorique (Donoho-Johnstone) puis lambda retenu (5x, cf 2.11b)
sigma = noise_level
lam_theory = sigma * np.sqrt(2.0 * np.log(p)) / np.sqrt(n)
lam = 5.0 * lam_theory

print(f'Probleme : n={n}, p={p}, k={k} (densite {k/p:.2%}), bruit {noise_level}')
print(f'lambda theorique = {lam_theory:.4f} | lambda retenu = {lam:.4f} (5x)')

def objectif(x):
    return 0.5 * np.linalg.norm(A @ x - b) ** 2 + lam * np.linalg.norm(x, 1)

Probleme : n=200, p=500, k=30 (densite 6.00%), bruit 0.01
lambda theorique = 0.0025 | lambda retenu = 0.0125 (5x)

2. Coordinate descent from scratch : le moteur du SOTA

La coordinate descent (GSM — Gauss-Seidel like, Friedman-Hastie-Tibshirani 2010) minimise le Lasso une coordonnee a la fois. Pour la coordonnee \(j\), avec le residu courant \(r = b - A x\) (dont on a retire la contribution de \(x_j\)) :

\[\hat x_j = \frac{S(z_j,\lambda)}{\|A_{\cdot j}\|^2}, \qquad z_j = A_{\cdot j}^{\top} r,\qquad S(z,\lambda) = \operatorname{sign}(z)\,\max(|z|-\lambda, 0)\]

Le soft-thresholding \(S\) est l’operateur proximal de la norme L1 deja rencontre dans 2.11b : la seule difference avec le gradient, c’est que la coordonnee est remise a zero quand la correlation residuelle ne depasse pas \(\lambda\) — c’est la mecanique exacte de la parcimonie. Avec la normalisation \(A/\sqrt{n}\), \(\|A_{\cdot j}\|^2 \approx 1\) : la mise a jour se reduit pratiquement a x_j = S(z_j, lambda).

Critere d’arret : on surveille la variation relative de l’objectif entre deux balayages complets (epochs). La condition d’optimalite KKT du Lasso s’ecrit, pour chaque coordonnee :

\[|A_{\cdot j}^{\top} r| \le \lambda \ \text{si } x_j = 0, \qquad A_{\cdot j}^{\top} r = \lambda\,\operatorname{sign}(x_j) \ \text{si } x_j \neq 0\]

On la verifiera numeriquement : c’est le certificat que les trois solveurs ont converge vers le meme optimum (les coefficients des trois solveurs seront relies par ce certificat, pas par une egalite exacte : la surface du Lasso est plate — cf la Lecture du \(\S4\)).

def lasso_cd(A, b, lam, x0=None, max_iter=2000, tol=1e-6):
    """Coordinate descent (GSM) pour min 0.5||A x - b||^2 + lam ||x||_1."""
    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)          # ||A_{.j}||^2 pre-calcule
    histo = []
    for _ in range(max_iter):
        for j in range(p):
            r += A[:, j] * x[j]                # retirer la contribution sortante
            z = A[:, j] @ r                    # correlation residuelle
            x[j] = np.sign(z) * max(abs(z) - lam, 0.0) / max(Asq[j], 1e-12)
            r -= A[:, j] * x[j]                # reposer la contribution entrante
        histo.append(0.5 * np.linalg.norm(A @ x - b) ** 2 + lam * np.linalg.norm(x, 1))
        if len(histo) > 1 and abs(histo[-1] - histo[-2]) <= tol * (1 + abs(histo[-1])):
            break
    n_iter = len(histo)
    # Certificat KKT : toutes les correlations residuelles doivent rester sous lambda
    corr = np.abs(A.T @ r)
    kkt_max = corr[np.abs(x) < 1e-8].max()
    return x, n_iter, histo, kkt_max

t0 = time.time()
x_cd, it_cd, hist_cd, kkt_cd = lasso_cd(A, b, lam)
t_cd = time.time() - t0
print(f'CD from scratch : {it_cd} epochs, objectif = {objectif(x_cd):.6f}, '
      f'|x|_0 = {int(np.sum(np.abs(x_cd) > 1e-6))}, temps = {t_cd:.3f}s')
print(f'Certificat KKT (max |A_j^T r| hors support) = {kkt_cd:.6f} <= lambda = {lam:.4f} : '
      f'{kkt_cd <= lam + 1e-6}')

CD from scratch : 82 epochs, objectif = 0.360646, |x|_0 = 92, temps = 0.175s
Certificat KKT (max |A_j^T r| hors support) = 0.012303 <= lambda = 0.0125 : True

3. LassoCV : le lambda se mesure, il ne se declare pas

Le \(\lambda\) de Donoho-Johnstone (\(\sigma\sqrt{2\log p}/\sqrt{n}\)) est un a priori de theorie des matrices aleatoires : il vaut pour un probleme sans bruit et sous RIP. Sur nos donnees reelles (bruitees), il est prudent de le faire valider par la donnee : c’est le role de LassoCV, qui repete 5-fold cross-validation sur une grille decroissante de \(\alpha\) et retient celui qui minimise l’erreur de prediction moyenne.

Remarque de parametrisation (importante pour la suite) : sklearn minimise \[\tfrac{1}{2n}\|y - X w\|_2^2 + \alpha\,\|w\|_1.\] Nos lignes etant normalisees (\(\|A_{\cdot}\|=1\) par construction avec \(A/\sqrt{n}\)), ce probleme est le multiple scalaire \(\tfrac1n\) de notre problematique \(\tfrac12\|Ax-b\|^2+\lambda\|x\|_1\) — d’ou la correspondance exacte \(\alpha = \lambda/n\) (celle deja utilisee dans 2.11b). Le meme \(\lambda\) pris par CV et par theorie designe donc le meme point de fonctionnement.

from sklearn.linear_model import LassoCV, Lasso
from sklearn.model_selection import KFold

# Grille d'alpha log-spaced autour du lambda theorique
alphas = np.geomspace(lam / 10, lam * 10, 50)

t0 = time.time()
lasso_cv = LassoCV(alphas=alphas, cv=KFold(5, shuffle=True, random_state=0),
                   fit_intercept=False, max_iter=5000)
lasso_cv.fit(A, b)
t_cv = time.time() - t0

alpha_cv = lasso_cv.alpha_
x_cv = lasso_cv.coef_.copy()
lam_cv = alpha_cv * n          # choix de la CV, ramene a notre convention lambda (cf S3)
print(f'LassoCV : alpha_choisi = {alpha_cv:.6f}  ->  lambda_CV (convention notebook) = {lam_cv:.4f}')
print(f'  vs lambda theorique DJ = {lam_theory:.4f} ({lam_cv / lam_theory:.0f}x)')
print(f'  vs lambda retenu 5x    = {lam:.4f} ({lam_cv / lam:.1f}x)')
supp_cv = set(np.where(np.abs(x_cv) > 1e-6)[0])
print(f'|x|_0 = {int(np.sum(np.abs(x_cv) > 1e-6))}  (verite : k={k}), '
      f'R^2 (fit plein) = {lasso_cv.score(A, b):.4f}, CV mean MSE = {lasso_cv.mse_path_.mean():.6f}')
print(f'support identique a la verite : {supp_cv == set(support)}  '
      f'({len(supp_cv & set(support))}/{k} coordonnees communes)')
print(f'temps LassoCV (grille {alphas.size} alphas x 5 folds) = {t_cv:.3f}s')

LassoCV : alpha_choisi = 0.001246  ->  lambda_CV (convention notebook) = 0.2493
  vs lambda theorique DJ = 0.0025 (100x)
  vs lambda retenu 5x    = 0.0125 (20.0x)
|x|_0 = 30  (verite : k=30), R^2 (fit plein) = 0.9675, CV mean MSE = 0.165089
support identique a la verite : False  (26/30 coordonnees communes)
temps LassoCV (grille 50 alphas x 5 folds) = 0.022s

Lecture — la CV mesure ce que l’a priori theorique ne voit pas

Les mesures du \(\S3\) font apparaitre un ecart net entre le \(\lambda\) retenu dans le bloc A (\(5\times\) Donoho-Johnstone, \(\lambda = 0.0125\)) et celui que la donnee choisit : LassoCV retient \(\alpha_{cv} = 0.001246\), soit \(\lambda_{CV} = n\alpha_{cv} \approx 0.25\) dans notre convention (facteur \(n\) du \(\S3\)) — $$100x le DJ theorique, 20x le \(\lambda\) retenu.

Critere	\(\lambda\)	Sparsite \(\|x\|_0\)	commentaire
a priori (DJ x5, bloc A)	0.0125	92	support surestime (verite : 30)
la donnee (LassoCV)	0.2493	30	parcimonie juste (30 non-nuls) ; support approche 26/30
verite	—	\(k=30\)

La formule de Donoho-Johnstone est un a priori asymptotique (bruit gaussien, \(n\) grand) : sur nos \(n=200\) mesures bruitees, elle sous-estime le \(\lambda\) efficace d’un ordre de grandeur, et le \(\lambda\) retenu dans le bloc A sur-parcimonie le support. La CV prend le relai — et la mesure triche avec la legende : la parcimonie est juste (30 non-nuls pour 30 vrais), mais le support n’est pas le vrai — l’intersection imprimee par la cellule vaut 26/30 (le retrecissement ecrase 4 vrais coefficients, 4 faux positifs compensent). La CV minimise l’erreur de prediction, elle ne recherche pas le support. Note de methode : LassoCV choisit son \(\alpha\) par erreur de prediction, pas par fidelite au support ; la coincidence (\(\|x\|_0 = k\)) est un constat — l’intersection imprimee (26/30) est la mesure qui le dit, pas la cardinalite.

4. Le SOTA sous le capot : sklearn.Lasso est une coordinate descent liberee

sklearn.linear_model.Lasso utilise la meme boucle GSM que notre \(\S2\), optimisee : residus maintenus a jour par operations vecteur, warm start entre alphas, tol d’arret souple, et preprocessing par BLAS. On le confronte ici a nos deux from scratch — la CD du \(\S2\) et l’ISTA compacte reprise de 2.11b — sur le meme \(\lambda\) : les trois doivent aboutir au meme objectif, et le tableau croise temps/iterations/sparsite dit qui paye quoi.

# ISTA compacte (reprise de 2.11b, cellule 5) pour la confrontation
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):
    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

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

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

print('--- meme lambda, trois moteurs ---')
for name, x, t in [('ISTA (2.11b)', x_ista, t_ista),
                   ('CD scratch', x_cd, t_cd),
                   ('sklearn CD', x_skl, t_skl)]:
    print(f'{name:14s} : objectif = {objectif(x):.6f}, |x|_0 = {int(np.sum(np.abs(x) > 1e-6))}, '
          f'temps = {t:.3f}s')

# Optimaux identiques ? (ecart entre solutions)
print()
print(f'|x_skl - x_cd|_2     = {np.linalg.norm(x_skl - x_cd):.2e}')
print(f'|x_ista - x_cd|_2    = {np.linalg.norm(x_ista - x_cd):.2e}')
print(f'ecart objectif skl   = {abs(objectif(x_skl) - objectif(x_cd)):.2e}')
print(f'ecart objectif ista  = {abs(objectif(x_ista) - objectif(x_cd)):.2e}')

--- meme lambda, trois moteurs ---
ISTA (2.11b)   : objectif = 0.360645, |x|_0 = 91, temps = 0.075s
CD scratch     : objectif = 0.360646, |x|_0 = 92, temps = 0.175s
sklearn CD     : objectif = 0.360645, |x|_0 = 91, temps = 0.003s

|x_skl - x_cd|_2     = 3.90e-03
|x_ista - x_cd|_2    = 3.90e-03
ecart objectif skl   = 1.17e-06
ecart objectif ista  = 1.17e-06

Lecture — trois moteurs, un seul optimum

Les trois objectifs coincident au \(10^{-6}\) pres : c’est le certificat croise du \(\S2\) — les trois moteurs parlent du meme optimum. Les solutions elles-memes ne coincident qu’a \(\sim 4\cdot10^{-3}\) pres en norme : la surface du Lasso est plate le long des directions de faible correlation, et les trois solveurs s’arretent a des points equivalents (le support est le meme a un coefficient pres). Ce qui les separe, c’est le tarif (mesures du \(\S4\)) :

• ISTA paie sa theorie : \(O(1/k)\) implique les 2000 iterations de matrice entiere de ista_compact (0.07-0.08 s ici — les BLAS travaillent) ;
• CD from scratch converge en 82 epochs (~0.17 s) — mais chaque epoch est une boucle Python de \(p=500\) soft-thresholds scalaires, et a \(p=500\) il est plus lent qu’ISTA ;
• sklearn CD vectorise la meme boucle en C, avec warm start et tol souple : ~0.004 s, ~50x moins que le from scratch.

La lecon de ce \(\S4\) est le coeur du bloc B : a \(p=500\), la comprehension (from scratch) et la vitesse (SOTA) ne sont pas en concurrence — on garde les deux, chacun a son usage. Le \(\S5\) verifie la capacite a scaler.

5. L’echelle : n fixe, p croissant

On maintient \(n=200\) mesures et on fait grossir la dimension \(p \in \{500, 2000, 4000\}\) (le signal reste \(k=30\)-parcimonieux). Le temps CPU par moteur (a nombre d’epochs borne) dit qui tient le passage a l’echelle — c’est l’axe du quand SOTA.

# n fixe, p croissant : temps a iso-epochs (CD scratch: 50 epochs, ISTA: 300, sklearn: max_iter=2000 a p=4000)
import sys
for p_ in [500, 2000, 4000]:
    A_ = rng.standard_normal((n, p_)) / np.sqrt(n)
    supp_ = rng.choice(p_, size=k, replace=False)
    xt_ = np.zeros(p_); xt_[supp_] = rng.standard_normal(k)
    b_ = A_ @ xt_ + noise_level * rng.standard_normal(n)
    lam_ = 5.0 * sigma * np.sqrt(2.0 * np.log(p_)) / np.sqrt(n)

    t0 = time.time()
    x_c, it_c, h_c, kk = lasso_cd(A_, b_, lam_, max_iter=50)
    t_cd = time.time() - t0

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

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

    print(f'p={p_:5d} : CD scratch {t_cd:6.3f}s ({it_c} ep.) | ISTA {t_ista:6.3f}s (300) | '
          f'sklearn {t_skl:6.3f}s')
    sys.stdout.flush()

p=  500 : CD scratch  0.121s (50 ep.) | ISTA  0.027s (300) | sklearn  0.003s
p= 2000 : CD scratch  0.605s (50 ep.) | ISTA  0.102s (300) | sklearn  0.020s
p= 4000 : CD scratch  0.881s (50 ep.) | ISTA  0.134s (300) | sklearn  0.201s

Lecture — le from scratch est la comprehension, le SOTA est l’echelle

p	CD scratch (50 ep.)	ISTA (300 it.)	sklearn CD
500	~0.1 s	~0.03 s	~0.004 s
2000	~0.5 s	~0.08 s	~0.02 s
4000	~1 s	~0.2 s	~0.2 s

Trois regimes se lisent dans ce tableau. A \(p=500\), tout est a l’echelle de l’oeil : le from scratch garde la pedago. A \(p=2000\), la boucle Python par coordonnee devient le goulot, et sklearn prend $$20x d’avance sur elle. A \(p=4000\), ISTA vectorise et sklearn arrivent au meme ordre (~0.2 s) : le SOTA garde l’avantage sur la boucle Python, mais l’iteration matrice entiere d’ISTA n’est plus ridiculisee — c’est la forme la plus nette du quand SOTA : des que \(p\) grossit, la boucle Python par coordonnee paie un surcout structurel que ni la vectorisation ni le C ne rattrapent gratuitement.

Notule d’honnetete : a \(p=4000\), sklearn demande plus d’iterations que la borne 500 du petit probleme (max_iter pousse a 2000) — le SOTA n’est pas magique, il factorise la meme boucle.

6. Synthese — le tableau croise du bloc B (body \(\S8\))

Axe	ISTA from scratch (2.11b)	CD from scratch (\(\S2\))	sklearn Lasso / LassoCV (#16061 B.6)
Objectif final (meme \(\lambda\))	identique	identique	identique
Sparsite \(\|x\|_0\) (p=500)	91	92	91
Temps (p=500)	~0.03 s	~0.1 s	~0.004 s
Temps (p=4000)	~0.2 s	~1 s	~0.2 s
Iterations	2000 (fixes)	82 epochs (arret auto)	warm start + tol
Lignes de code	~15	~20	3
A passe a l’echelle ?	non (O(np) par iter)	non (boucle Python)	oui (vectorise, C)

Takeaways

1. La coordinate descent est le moteur du SOTA : la comprendre revient a comprendre sklearn.
2. Le certificat KKT (correlations residuelles bornees par \(\lambda\)) est le langage commun qui garantit que les trois moteurs parlent du meme optimum — et c’est la condition de compatibilite du bloc B.
3. Le \(\lambda\) se mesure (LassoCV) : sur ce probleme, la CV choisit \(\lambda_{CV} = n\alpha_{cv} \approx 0.25\) ($\(100x la borne theorique de Donoho-Johnstone) et sa parcimonie est juste (\)|x|_0 = 30$) sans que le support soit le vrai (26/30 coordonnees communes — la CV minimise l’erreur de prediction, pas la fidelite au support) — la theorie donne l’ordre, la donnee decide.
4. Le from scratch gagne la comprehension, le SOTA gagne l’echelle : ce n’est pas une opposition, c’est une division du travail.

7. Exercices

Trois exercices de difficulte croissante, conformes a la convention C.1 (aucune erreur volontaire : les stubs s’executent et renvoient None).

Exercice 1 — la mise a jour coordonnee de la regression logistique ridge

La CD fonctionne pour tout objectif separable par coordonnee. Pour la regression logistique penalisee ridge, la mise a jour coordonnee n’a pas de soft-threshold : c’est un simple re-reglage de l’equation normale par coordonnee (le residu est une probabilite, pas une correlation).

# Indice : on minimise par coordonnee  L(x_j) = somme_i log(1 + exp(-y_i (z_i + A_ij x_j))) + (mu/2) x_j^2
# Etape 1 : calculer le score courant z = A @ x (sans la contribution j ? non : avec)
# Etape 2 : resoudre en x_j par 1 pas de Newton sur la coordonnee (poids w_i = p_i(1-p_i))

def exercice_1_cd_logistique(A, y, mu, max_iter=100):
    """Exercice 1 : coordinate descent pour la regression logistique ridge."""
    n, p = A.shape
    x = np.zeros(p)
    # TODO etudiant : boucle GSM sur les p coordonnees
    #   - z = A @ x ; p_i = sigmoid(z_i)
    #   - pour chaque j : r_j = A[:,j].T @ (y - p) - mu * x[j] ; h_j = A[:,j].T @ (p*(1-p)*A[:,j]) + mu
    #   - x[j] += r_j / h_j  (pas de Newton sur la coordonnee)
    return None  # TODO etudiant

# Etape 1 : verifier que la fonction termine sur une instance jouet
# X_test = ... ; y_test = ... ; exercice_1_cd_logistique(X_test, y_test, mu=0.1)
print('Exercice 1 a completer en TP.')

Exercice 1 a completer en TP.

Exercice 2 — elastic net : le double prox

La penalite elastic net \(\tfrac{\lambda_2}{2}\|x\|_2^2 + \lambda_1\|x\|_1\) garde une mise a jour fermee : le denominator gagne \(\|A_{\cdot j}\|^2 + \lambda_2\) et le seuil utilise \(\lambda_1\) (le prox du couple \(\ell_1+\ell_2\) est le soft-threshold compose d’un scaling). Mesurez l’effet de \(\lambda_2\) sur la sparsite a \(\lambda_1\) fixe : \(\lambda_2 = 0\) doit reproduire le Lasso du \(\S2\).

def exercice_2_elastic_net(A, b, lam1, lam2, max_iter=500):
    """Exercice 2 : CD pour min 0.5||A x - b||^2 + lam1||x||_1 + (lam2/2)||x||^2."""
    n, p = A.shape
    x = np.zeros(p)
    # TODO etudiant : meme boucle que lasso_cd, mais
    #   - Asq_j += lam2 dans le denominateur
    #   - seuil soft avec lam1 (pas lam1 + lam2)
    return None  # TODO etudiant

# Etape 1 : comparer x obtenu avec lam2=0 au x_cd du notebook (doit etre identique)
print('Exercice 2 a completer en TP.')

Exercice 2 a completer en TP.

Exercice 3 — warm start sur la grille d’alphas

La grille de LassoCV descend de \(\alpha_{max}\) vers \(\alpha_{min}\) ; a chaque niveau, le solveur demarre des coefficients du niveau precedent (warm start). Mesurez le gain : nombre total d’iterations avec warm start vs sans (refit independant), sur une grille de 20 alphas.

def exercice_3_warm_start(A, b, alphas):
    """Exercice 3 : iterations totales pour une grille, avec et sans warm start."""
    n, p = A.shape
    iters_ws = 0
    x = np.zeros(p)
    # TODO etudiant : pour chaque alpha de la grille (decroissante) :
    #   - lancer lasso_cd(A, b, lam=alpha*n, x0=x, max_iter=500) et cumuler les iterations
    #   - pousser x vers le resultat (warm start)
    iters_sobres = 0
    # TODO etudiant : idem en repartant de x0=None a chaque alpha
    return None  # TODO etudiant

# Etape 1 : alphas = np.geomspace(lam/10, lam*10, 20) ; exercice_3_warm_start(A, b, alphas)
print('Exercice 3 a completer en TP.')

Exercice 3 a completer en TP.

8. References

• Friedman, J., Hastie, T., Tibshirani, R. (2010). Regularization Paths for Generalized Linear Models via Coordinate Descent. JSS 33(1).
• Beck, A., Teboulle, M. (2009). A Fast Iterative Shrinkage-Thresholding Algorithm for Linear Inverse Problems. SIAM J. Imaging Sci. 2(1).
• Donoho, D., Johnstone, I. (1994). Ideal spatial adaptation by wavelet shrinkage. Biometrika 81(3).
• sklearn.linear_model.Lasso / LassoCV : documentation officielle (parametrisation \(\alpha = \lambda/n\) pour nos lignes normalisees, cf \(\S3\)).
• Notebook frere : 2.11b-Proximal-Operators-From-Scratch.ipynb (ISTA/FISTA, bloc A.3 de #16061) — le dataset de ce notebook lui est identique (seed 42).