From 4af3950739c3b76e409f4db32af324f2cafaa017 Mon Sep 17 00:00:00 2001 From: jsboige Date: Thu, 24 Sep 2026 13:54:43 +0200 Subject: [PATCH] =?UTF-8?q?feat(notebooks,#16061):=202.11e-CVXPY-Optimisat?= =?UTF-8?q?ion=20=E2=80=94=20bloc=20B.7,=20Lasso=20declaratif?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Bloc B.7 de l'acceptance #16061 (« Notebook CVXPY-OPT ») : formulation declarative du Lasso, resolution CLARABEL/SCS/OSQP, comparaison aux solutions from scratch (2.11b/2.11c) et ADMM (2.11d) sur leurs problemes respectifs. Contenu mesure (36 cellules, 12 code, execution reelle 62.8s, 0 erreur) : - S2 : modele en une ligne, ruleset DCP verifie (is_dcp True/False sur variante concave) — certification de convexite a la compilation ; - S3 : canonicalisation mesuree via get_problem_data : 500 inconnues declarees deviennent 1200 variables canoniques (200 egalites + 1000 inegalites d'epigraphe), vecteur de cout c = 500 composantes valant toutes lambda ; - S4 : un modele, trois solveurs — CLARABEL 12 it., SCS 450 it., OSQP 2775 it., meme objectif, solutions a 4.6e-05 ; - S5 : confrontation — cvxpy 0.36064507, CD scratch 0.36064623, ISTA et sklearn 0.36064506 (ecarts <= 1.2e-06) : le verdict « meme optimum global » demande par l'acceptance, mesure et non affirme ; - S6 : dualite lue dans les variables duales (forme epigraphe) — stationnarite a 1.2e-15, complementarites a 7e-11, certificat KKT |c_j| <= lambda avec egalite sur le support, sans le recoder ; - S7 : chemin de 12 resolutions + seuil lambda_max = ||A^T b||_inf calcule (et non cherche), avec la degenerescence au seuil mesuree par l'ecart d'objectif ; - S8 : ADMM de 2.11d rejoue sur son propre design correle (reprise exacte de design_sparse_corr/admm_lasso) confronte a l'optimum certifie : rho=20 converge en 178 it. a 2.5e-04 de f*, rho=1 et 5 n'atteignent pas la borne ; - S9 : echelle n=200, p jusqu'a 4000 — le declaratif perd la course de vitesse (8.9s contre 0.315s pour sklearn) et gagne celle du changement de probleme. Licence pedagogique : meme probleme que 2.11b/2.11c (seed 42) pour la comparabilite, 3 exercices (group Lasso, forme contrainte et dual, lambda_max par le solveur), C.1 strict (stubs qui s'executent), C.2 strict (outputs reels). Note : ECOS n'est plus distribue avec cvxpy 1.9 (mesure : solveurs installes = CLARABEL/SCS/OSQP/HIGHS/GLOP/PDLP/SCIPY) — la lignee ECOS est portee par CLARABEL, ce que le notebook documente en S1 et en references. Co-authored-by: Claude-Code --- .../2.11e-CVXPY-Optimisation.ipynb | 1619 +++++++++++++++++ 1 file changed, 1619 insertions(+) create mode 100644 MyIA.AI.Notebooks/ML/DataScienceWithAgents/02-ML-Cours/2.11e-CVXPY-Optimisation.ipynb diff --git a/MyIA.AI.Notebooks/ML/DataScienceWithAgents/02-ML-Cours/2.11e-CVXPY-Optimisation.ipynb b/MyIA.AI.Notebooks/ML/DataScienceWithAgents/02-ML-Cours/2.11e-CVXPY-Optimisation.ipynb new file mode 100644 index 0000000000..cec431bfa2 --- /dev/null +++ b/MyIA.AI.Notebooks/ML/DataScienceWithAgents/02-ML-Cours/2.11e-CVXPY-Optimisation.ipynb @@ -0,0 +1,1619 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "4e5c66e0", + "metadata": { + "papermill": { + "duration": 0.0073, + "end_time": "2026-09-24T11:51:33.380436+00:00", + "exception": false, + "start_time": "2026-09-24T11:51:33.373136+00:00", + "status": "completed" + }, + "tags": [] + }, + "source": [ + "# 2.11e — CVXPY : le problème déclaré, le solveur choisi\n", + "\n", + "**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.\n", + "\n", + "- la même formulation Lasso écrite en **une ligne** avec `cvxpy`, sur **le même problème** que `2.11b`/`2.11c` (seed 42), donc comparable nombre à nombre ;\n", + "- ce que la couche de modélisation **fabrique** avant d'appeler le solveur : la reformulation canonique, lue sur les données réellement transmises ;\n", + "- **un problème, plusieurs solveurs** — CLARABEL (point intérieur), SCS (premier ordre), OSQP (QP) : même optimum, tarifs différents ;\n", + "- la **dualité lue dans les variables duales** : le certificat KKT des blocs A et B sort directement du solveur, sans le recoder ;\n", + "- 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 ?\n", + "\n", + "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).\n" + ] + }, + { + "cell_type": "markdown", + "id": "45f0e970", + "metadata": { + "papermill": { + "duration": 0.007976, + "end_time": "2026-09-24T11:51:33.396084+00:00", + "exception": false, + "start_time": "2026-09-24T11:51:33.388108+00:00", + "status": "completed" + }, + "tags": [] + }, + "source": [ + "## 1. Le même problème que les blocs A et B\n", + "\n", + "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$).\n", + "\n", + "Deux échelles de $\\lambda$ cohabitent dans ce notebook, et il faut les garder distinctes :\n", + "\n", + "| Symbole | Valeur | Rôle |\n", + "|---|---|---|\n", + "| $\\lambda_{DJ} = \\sigma\\sqrt{2\\log p}/\\sqrt{n}$ | $0{,}0025$ | l'a priori de Donoho-Johnstone |\n", + "| $\\lambda = 5\\lambda_{DJ}$ | $0{,}0125$ | le point de fonctionnement retenu dans `2.11b`/`2.11c` |\n", + "| $\\lambda_{\\max} = \\|A^{\\top}b\\|_\\infty$ | $3{,}19$ | le seuil d'annulation exact (cf $\\S7$) |\n" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "5a4a8993", + "metadata": { + "execution": { + "iopub.execute_input": "2026-09-24T11:51:33.419731Z", + "iopub.status.busy": "2026-09-24T11:51:33.419186Z", + "iopub.status.idle": "2026-09-24T11:51:46.259812Z", + "shell.execute_reply": "2026-09-24T11:51:46.258283Z" + }, + "papermill": { + "duration": 12.858273, + "end_time": "2026-09-24T11:51:46.261072+00:00", + "exception": false, + "start_time": "2026-09-24T11:51:33.402799+00:00", + "status": "completed" + }, + "tags": [] + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Probleme : n=200, p=500, k=30 (densite 6.00%), bruit=0.01, seed=42\n", + "lambda_DJ = 0.002493 | lambda retenu (5x) = 0.012465\n", + "lambda_max = max_j |A_j^T b| = 3.191173\n", + "cvxpy 1.9.2 | solveurs installes : ['CLARABEL', 'GLOP', 'HIGHS', 'OSQP', 'PDLP', 'SCIPY', 'SCS']\n" + ] + } + ], + "source": [ + "import numpy as np\n", + "import time\n", + "import cvxpy as cp\n", + "\n", + "np.random.seed(42)\n", + "\n", + "# --- Generation identique a 2.11b (cellules 9-10) et 2.11c (cellule 2) ---\n", + "n, p, k = 200, 500, 30\n", + "noise_level = 0.01\n", + "rng = np.random.default_rng(42)\n", + "A = rng.standard_normal((n, p)) / np.sqrt(n)\n", + "support = rng.choice(p, size=k, replace=False)\n", + "x_true = np.zeros(p)\n", + "x_true[support] = rng.standard_normal(k)\n", + "b = A @ x_true + noise_level * rng.standard_normal(n)\n", + "\n", + "# Les trois echelles de lambda (cf tableau du S1)\n", + "sigma = noise_level\n", + "lam_theory = sigma * np.sqrt(2.0 * np.log(p)) / np.sqrt(n)\n", + "lam = 5.0 * lam_theory\n", + "lam_max = np.max(np.abs(A.T @ b))\n", + "\n", + "\n", + "def objectif(x, lam_=None):\n", + " \"\"\"Objectif Lasso brut : 0.5||Ax - b||^2 + lam*||x||_1.\"\"\"\n", + " l = lam if lam_ is None else lam_\n", + " return 0.5 * np.linalg.norm(A @ x - b) ** 2 + l * np.abs(x).sum()\n", + "\n", + "\n", + "print(f'Probleme : n={n}, p={p}, k={k} (densite {k/p:.2%}), bruit={noise_level}, seed=42')\n", + "print(f'lambda_DJ = {lam_theory:.6f} | lambda retenu (5x) = {lam:.6f}')\n", + "print(f'lambda_max = max_j |A_j^T b| = {lam_max:.6f}')\n", + "print(f'cvxpy {cp.__version__} | solveurs installes : {sorted(cp.installed_solvers())}')\n" + ] + }, + { + "cell_type": "markdown", + "id": "fc7aa0dc", + "metadata": { + "papermill": { + "duration": 0.006733, + "end_time": "2026-09-24T11:51:46.273996+00:00", + "exception": false, + "start_time": "2026-09-24T11:51:46.267263+00:00", + "status": "completed" + }, + "tags": [] + }, + "source": [ + "## 2. Le modeling déclaratif : dire le problème, pas la méthode\n", + "\n", + "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** :\n", + "\n", + "$$\\min_x \\; \\tfrac12 \\|Ax - b\\|_2^2 + \\lambda \\|x\\|_1$$\n", + "\n", + "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 :\n", + "\n", + "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 ;\n", + "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 ;\n", + "3. **la résolution** — un solveur nommé, appelé par son nom, qui reçoit une forme canonique et rend une solution.\n", + "\n", + "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.\n" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "5044f97e", + "metadata": { + "execution": { + "iopub.execute_input": "2026-09-24T11:51:46.288495Z", + "iopub.status.busy": "2026-09-24T11:51:46.287964Z", + "iopub.status.idle": "2026-09-24T11:51:46.989739Z", + "shell.execute_reply": "2026-09-24T11:51:46.988643Z" + }, + "papermill": { + "duration": 0.710992, + "end_time": "2026-09-24T11:51:46.991607+00:00", + "exception": false, + "start_time": "2026-09-24T11:51:46.280615+00:00", + "status": "completed" + }, + "tags": [] + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Le probleme est-il conforme au ruleset DCP ? True\n", + "Avec -lam*norm1 (objectif concave) : DCP = False -> refuse avant toute evaluation\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + "CLARABEL : status = optimal\n", + " objectif = 0.36064507\n", + " |x|_0 (1e-6) = 93 (verite : k=30)\n", + " iterations = 12\n", + " temps (mur) = 0.689s\n" + ] + } + ], + "source": [ + "# Le probleme, ecrit tel quel. Aucune instruction d'algorithme.\n", + "x_model = cp.Variable(p)\n", + "cout = 0.5 * cp.sum_squares(A @ x_model - b) + lam * cp.norm1(x_model)\n", + "prob_lasso = cp.Problem(cp.Minimize(cout))\n", + "\n", + "# Le ruleset DCP certifie la convexite avant tout calcul.\n", + "print('Le probleme est-il conforme au ruleset DCP ? ', prob_lasso.is_dcp())\n", + "\n", + "# Garde-fou pedagogique : inverser le signe du terme L1 rend l'objectif concave.\n", + "prob_concave = cp.Problem(\n", + " cp.Minimize(0.5 * cp.sum_squares(A @ x_model - b) - lam * cp.norm1(x_model)))\n", + "print('Avec -lam*norm1 (objectif concave) : DCP =', prob_concave.is_dcp(),\n", + " '-> refuse avant toute evaluation')\n", + "\n", + "# Resolution par un solveur nomme\n", + "t0 = time.time()\n", + "prob_lasso.solve(solver=cp.CLARABEL)\n", + "t_clarabel = time.time() - t0\n", + "x_cvx = x_model.value\n", + "\n", + "print()\n", + "print(f'CLARABEL : status = {prob_lasso.status}')\n", + "print(f' objectif = {objectif(x_cvx):.8f}')\n", + "print(f' |x|_0 (1e-6) = {int(np.sum(np.abs(x_cvx) > 1e-6))} (verite : k={k})')\n", + "print(f' iterations = {prob_lasso.solver_stats.num_iters}')\n", + "print(f' temps (mur) = {t_clarabel:.3f}s')\n" + ] + }, + { + "cell_type": "markdown", + "id": "58175af5", + "metadata": { + "papermill": { + "duration": 0.006499, + "end_time": "2026-09-24T11:51:47.005964+00:00", + "exception": false, + "start_time": "2026-09-24T11:51:46.999465+00:00", + "status": "completed" + }, + "tags": [] + }, + "source": [ + "### Lecture — une ligne de modèle, et une solution au sixième chiffre\n", + "\n", + "`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.\n", + "\n", + "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.\n", + "\n", + "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.\n" + ] + }, + { + "cell_type": "markdown", + "id": "5705bd14", + "metadata": { + "papermill": { + "duration": 0.007327, + "end_time": "2026-09-24T11:51:47.019787+00:00", + "exception": false, + "start_time": "2026-09-24T11:51:47.012460+00:00", + "status": "completed" + }, + "tags": [] + }, + "source": [ + "## 3. Ce que le solveur voit : la reformulation canonique\n", + "\n", + "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 :\n", + "\n", + "$$\\min_x \\; \\tfrac12 x^{\\top}Px + c^{\\top}x \\quad \\text{s.c.} \\quad Ax + s = b, \\; s \\in K$$\n", + "\n", + "Deux reformulations mécaniques ont lieu, et `get_problem_data` les rend **mesurables** plutôt que mystérieuses :\n", + "\n", + "- $\\|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 ;\n", + "- 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é.\n", + "\n", + "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.\n" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "b1add0cf", + "metadata": { + "execution": { + "iopub.execute_input": "2026-09-24T11:51:47.036728Z", + "iopub.status.busy": "2026-09-24T11:51:47.036318Z", + "iopub.status.idle": "2026-09-24T11:51:47.048642Z", + "shell.execute_reply": "2026-09-24T11:51:47.047234Z" + }, + "papermill": { + "duration": 0.023073, + "end_time": "2026-09-24T11:51:47.049895+00:00", + "exception": false, + "start_time": "2026-09-24T11:51:47.026822+00:00", + "status": "completed" + }, + "tags": [] + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Ce que CLARABEL recoit reellement :\n", + " variables canoniques : 1200 (le modele en declare 500)\n", + " lignes canoniques : 1200\n", + " contraintes : 200 equalities, 1000 inequalities, 0 exponential cones, \n", + "SOC constraints: [], PSD constraints: [],\n", + " 3d power cones [], [].\n", + " A canonique : (1200, 1200), 102200 non-nuls (7.1% de densite, stockage creux)\n", + " P (partie quadratique): (1200, 1200), 200 non-nuls\n", + " cout lineaire c : 500 composantes non nulles, toutes egales a 0.01246456 (= lambda)\n" + ] + } + ], + "source": [ + "data_conique, chaine_canonique, inverse = prob_lasso.get_problem_data(cp.CLARABEL)\n", + "Ac, Pc, cc = data_conique['A'], data_conique['P'], data_conique['c']\n", + "\n", + "print('Ce que CLARABEL recoit reellement :')\n", + "print(f' variables canoniques : {Ac.shape[1]} (le modele en declare {p})')\n", + "print(f' lignes canoniques : {Ac.shape[0]}')\n", + "print(f' contraintes : {data_conique[\"dims\"]}')\n", + "print(f' A canonique : {Ac.shape}, {Ac.nnz} non-nuls '\n", + " f'({100 * Ac.nnz / (Ac.shape[0] * Ac.shape[1]):.1f}% de densite, stockage creux)')\n", + "print(f' P (partie quadratique): {Pc.shape}, {Pc.nnz} non-nuls')\n", + "cnz = np.asarray(cc)[np.asarray(cc) != 0]\n", + "print(f' cout lineaire c : {cnz.size} composantes non nulles, '\n", + " f'toutes egales a {np.unique(np.round(cnz, 8))[0]:.8f} (= lambda)')\n" + ] + }, + { + "cell_type": "markdown", + "id": "a6a51a6b", + "metadata": { + "papermill": { + "duration": 0.006826, + "end_time": "2026-09-24T11:51:47.064084+00:00", + "exception": false, + "start_time": "2026-09-24T11:51:47.057258+00:00", + "status": "completed" + }, + "tags": [] + }, + "source": [ + "### Lecture — 500 inconnues deviennent 1200 variables\n", + "\n", + "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).\n", + "\n", + "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 %.\n", + "\n", + "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 écrit en une ligne 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.\n" + ] + }, + { + "cell_type": "markdown", + "id": "d9a185f5", + "metadata": { + "papermill": { + "duration": 0.007717, + "end_time": "2026-09-24T11:51:47.079208+00:00", + "exception": false, + "start_time": "2026-09-24T11:51:47.071491+00:00", + "status": "completed" + }, + "tags": [] + }, + "source": [ + "## 4. Un problème, trois solveurs : le choix est un paramètre\n", + "\n", + "Le même modèle, résolu par trois solveurs qui appartiennent à **trois familles algorithmiques distinctes** :\n", + "\n", + "| Solveur | Famille | Ce qu'il fait du problème canonique |\n", + "|---|---|---|\n", + "| **CLARABEL** | point intérieur primal-dual | peu d'itérations, algèbre linéaire dense, haute précision |\n", + "| **SCS** | premier ordre (splitting d'opérateurs) | beaucoup d'itérations bon marché, précision réglable |\n", + "| **OSQP** | ADMM (quadratique) | itérations très bon marché sur la forme QP |\n", + "\n", + "`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$.\n" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "ba2ea341", + "metadata": { + "execution": { + "iopub.execute_input": "2026-09-24T11:51:47.094300Z", + "iopub.status.busy": "2026-09-24T11:51:47.093817Z", + "iopub.status.idle": "2026-09-24T11:51:49.381211Z", + "shell.execute_reply": "2026-09-24T11:51:49.379986Z" + }, + "papermill": { + "duration": 2.297801, + "end_time": "2026-09-24T11:51:49.382471+00:00", + "exception": false, + "start_time": "2026-09-24T11:51:47.084670+00:00", + "status": "completed" + }, + "tags": [] + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "solveur status objectif |x|_0 iterations temps |x-x_CLARABEL|_2\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "CLARABEL optimal 0.36064507 93 12 0.828s 0.00e+00\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "SCS optimal 0.36064506 91 450 0.325s 4.63e-05\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "OSQP optimal 0.36064506 91 2775 1.117s 4.63e-05\n" + ] + } + ], + "source": [ + "parametres = {\n", + " 'CLARABEL': {},\n", + " 'SCS': {'eps': 1e-9, 'max_iters': 200000},\n", + " 'OSQP': {'eps_abs': 1e-9, 'eps_rel': 1e-9, 'max_iter': 200000},\n", + "}\n", + "\n", + "print(f'{\"solveur\":9s} {\"status\":8s} {\"objectif\":>12s} {\"|x|_0\":>6s} {\"iterations\":>10s} '\n", + " f'{\"temps\":>8s} {\"|x-x_CLARABEL|_2\":>16s}')\n", + "resultats = {}\n", + "for nom, kw in parametres.items():\n", + " v = cp.Variable(p)\n", + " pr = cp.Problem(cp.Minimize(0.5 * cp.sum_squares(A @ v - b) + lam * cp.norm1(v)))\n", + " t0 = time.time()\n", + " pr.solve(solver=nom, **kw)\n", + " t = time.time() - t0\n", + " xx = v.value\n", + " resultats[nom] = (pr.status, objectif(xx), int(np.sum(np.abs(xx) > 1e-6)),\n", + " pr.solver_stats.num_iters, t, np.linalg.norm(xx - x_cvx))\n", + " st, ob, nn, it, tt, ec = resultats[nom]\n", + " print(f'{nom:9s} {st:8s} {ob:12.8f} {nn:6d} {it:10d} {tt:7.3f}s {ec:16.2e}')\n" + ] + }, + { + "cell_type": "markdown", + "id": "37fdb32d", + "metadata": { + "papermill": { + "duration": 0.005777, + "end_time": "2026-09-24T11:51:49.396028+00:00", + "exception": false, + "start_time": "2026-09-24T11:51:49.390251+00:00", + "status": "completed" + }, + "tags": [] + }, + "source": [ + "### Lecture — même optimum, trois tarifs, trois récits d'itération\n", + "\n", + "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.\n", + "\n", + "Ce qui les sépare vraiment, c'est le **récit d'itération**, et il est structurel :\n", + "\n", + "- **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 ;\n", + "- **SCS** fait plusieurs centaines d'itérations d'opérateurs proximaux : chaque itération est une multiplication matrice-vecteur, très bon marché ;\n", + "- **OSQP** en fait plusieurs milliers, chacune encore moins chère, sur la forme quadratique du problème.\n", + "\n", + "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.\n", + "\n", + "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**.\n" + ] + }, + { + "cell_type": "markdown", + "id": "b99436ad", + "metadata": { + "papermill": { + "duration": 0.005852, + "end_time": "2026-09-24T11:51:49.407798+00:00", + "exception": false, + "start_time": "2026-09-24T11:51:49.401946+00:00", + "status": "completed" + }, + "tags": [] + }, + "source": [ + "## 5. La confrontation : cvxpy contre from scratch et sklearn, même $\\lambda$\n", + "\n", + "`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 ?*\n", + "\n", + "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$.\n" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "4d2ee674", + "metadata": { + "execution": { + "iopub.execute_input": "2026-09-24T11:51:49.424990Z", + "iopub.status.busy": "2026-09-24T11:51:49.424582Z", + "iopub.status.idle": "2026-09-24T11:51:51.314155Z", + "shell.execute_reply": "2026-09-24T11:51:51.312483Z" + }, + "papermill": { + "duration": 1.90104, + "end_time": "2026-09-24T11:51:51.315700+00:00", + "exception": false, + "start_time": "2026-09-24T11:51:49.414660+00:00", + "status": "completed" + }, + "tags": [] + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "moteur objectif |x|_0 temps |x-x_cvxpy|_2 detail\n", + "cvxpy/CLARABEL 0.36064507 93 0.689s 0.00e+00 12 iter.\n", + "CD from scratch (2.11c) 0.36064623 92 0.519s 3.89e-03 82 epochs\n", + "ISTA (2.11b) 0.36064506 91 0.146s 4.63e-05 2000 iter.\n", + "sklearn CD (2.11c) 0.36064506 91 0.010s 4.63e-05 116 iter.\n", + "\n", + "Ecarts d objectif au solveur conique (le critere qui ne depend pas du seuil) :\n", + " CD from scratch (2.11c) : +1.167e-06\n", + " ISTA (2.11b) : -2.789e-09\n", + " sklearn CD (2.11c) : -2.789e-09\n" + ] + } + ], + "source": [ + "def lasso_cd(A, b, lam, x0=None, max_iter=2000, tol=1e-6):\n", + " \"\"\"Coordinate descent (GSM) — reprise de 2.11c S2, meme boucle, meme critere d'arret.\"\"\"\n", + " n_, p_ = A.shape\n", + " x = np.zeros(p_) if x0 is None else x0.copy()\n", + " r = b - A @ x\n", + " Asq = np.einsum('ij,ij->j', A, A)\n", + " obj_local = lambda xx: 0.5 * np.linalg.norm(A @ xx - b) ** 2 + lam * np.abs(xx).sum()\n", + " histo = []\n", + " for _ in range(max_iter):\n", + " for j in range(p_):\n", + " r += A[:, j] * x[j]\n", + " z = A[:, j] @ r\n", + " x[j] = np.sign(z) * max(abs(z) - lam, 0.0) / max(Asq[j], 1e-12)\n", + " r -= A[:, j] * x[j]\n", + " histo.append(obj_local(x))\n", + " if len(histo) > 1 and abs(histo[-1] - histo[-2]) <= tol * (1 + abs(histo[-1])):\n", + " break\n", + " return x, len(histo)\n", + "\n", + "\n", + "def prox_l1(x, t):\n", + " return np.sign(x) * np.maximum(np.abs(x) - t, 0.0)\n", + "\n", + "\n", + "def ista_compact(A, b, lam, n_iter=2000):\n", + " \"\"\"ISTA — reprise de 2.11b cellule 5 / 2.11c S4 : gradient + soft-thresholding.\"\"\"\n", + " L = np.linalg.norm(A, 2) ** 2\n", + " x = np.zeros(A.shape[1])\n", + " for _ in range(n_iter):\n", + " x = prox_l1(x - A.T @ (A @ x - b) / L, lam / L)\n", + " return x\n", + "\n", + "\n", + "from sklearn.linear_model import Lasso\n", + "\n", + "t0 = time.time()\n", + "x_cd, it_cd = lasso_cd(A, b, lam)\n", + "t_cd = time.time() - t0\n", + "\n", + "t0 = time.time()\n", + "x_ista = ista_compact(A, b, lam)\n", + "t_ista = time.time() - t0\n", + "\n", + "t0 = time.time()\n", + "skl = Lasso(alpha=lam / n, fit_intercept=False, max_iter=2000, tol=1e-8) # alpha = lam/n (cf 2.11c S3)\n", + "skl.fit(A, b)\n", + "t_skl = time.time() - t0\n", + "x_skl = skl.coef_.copy()\n", + "\n", + "moteurs = [('cvxpy/CLARABEL', x_cvx, t_clarabel, f'{prob_lasso.solver_stats.num_iters} iter.'),\n", + " ('CD from scratch (2.11c)', x_cd, t_cd, f'{it_cd} epochs'),\n", + " ('ISTA (2.11b)', x_ista, t_ista, '2000 iter.'),\n", + " ('sklearn CD (2.11c)', x_skl, t_skl, f'{skl.n_iter_} iter.')]\n", + "\n", + "print(f'{\"moteur\":26s} {\"objectif\":>12s} {\"|x|_0\":>6s} {\"temps\":>8s} {\"|x-x_cvxpy|_2\":>14s} detail')\n", + "for nom, xx, tt, det in moteurs:\n", + " print(f'{nom:26s} {objectif(xx):12.8f} {int(np.sum(np.abs(xx) > 1e-6)):6d} '\n", + " f'{tt:7.3f}s {np.linalg.norm(xx - x_cvx):14.2e} {det}')\n", + "\n", + "x_cd_ref = resultats['CLARABEL'][1]\n", + "print()\n", + "print('Ecarts d objectif au solveur conique (le critere qui ne depend pas du seuil) :')\n", + "for nom, xx, tt, det in moteurs[1:]:\n", + " print(f' {nom:26s} : {objectif(xx) - x_cd_ref:+.3e}')\n" + ] + }, + { + "cell_type": "markdown", + "id": "90568e8c", + "metadata": { + "papermill": { + "duration": 0.006344, + "end_time": "2026-09-24T11:51:51.330140+00:00", + "exception": false, + "start_time": "2026-09-24T11:51:51.323796+00:00", + "status": "completed" + }, + "tags": [] + }, + "source": [ + "### Lecture — quatre moteurs, un seul optimum global\n", + "\n", + "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.\n", + "\n", + "| Critère | cvxpy/CLARABEL | CD from scratch | ISTA | sklearn |\n", + "|---|---|---|---|---|\n", + "| Objectif ($\\lambda = 0{,}0125$) | 0,360645 | 0,360646 | 0,360645 | 0,360645 |\n", + "| Statut | optimal (certifié DCP) | certificat KKT | point fixe proximal | optimal |\n", + "| Itérations | ~12 | ~82 epochs | 2000 (fixes) | warm start |\n", + "\n", + "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 demande trois lignes — c'est le sujet des exercices du $\\S11$.\n" + ] + }, + { + "cell_type": "markdown", + "id": "7163163b", + "metadata": { + "papermill": { + "duration": 0.007043, + "end_time": "2026-09-24T11:51:51.347063+00:00", + "exception": false, + "start_time": "2026-09-24T11:51:51.340020+00:00", + "status": "completed" + }, + "tags": [] + }, + "source": [ + "## 6. La dualité sort du solveur : le certificat KKT sans le recoder\n", + "\n", + "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\n", + "\n", + "$$\\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$$\n", + "\n", + "et de lire `constraint.dual_value`. Les multiplicateurs $\\mu^{+}, \\mu^{-} \\ge 0$ des deux familles de contraintes portent toute l'information duale :\n", + "\n", + "- **stationnarité** : la corrélation résiduelle s'écrit $c = A^{\\top}(b - Ax) = \\mu^{+} - \\mu^{-}$ ;\n", + "- **complémentarité** : $\\mu^{+}_j (t_j - x_j) = 0$ et $\\mu^{-}_j (t_j + x_j) = 0$ ;\n", + "- d'où la **condition KKT du Lasso**, retrouvée sans l'écrire : $|c_j| \\le \\lambda$, avec égalité exactement sur le support.\n", + "\n", + "Autrement dit : la théorie de `2.11c` ne se réimplémente pas ici, elle se **lit**.\n" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "0cac9be0", + "metadata": { + "execution": { + "iopub.execute_input": "2026-09-24T11:51:51.361861Z", + "iopub.status.busy": "2026-09-24T11:51:51.361361Z", + "iopub.status.idle": "2026-09-24T11:51:52.015728Z", + "shell.execute_reply": "2026-09-24T11:51:52.014407Z" + }, + "papermill": { + "duration": 0.663251, + "end_time": "2026-09-24T11:51:52.016679+00:00", + "exception": false, + "start_time": "2026-09-24T11:51:51.353428+00:00", + "status": "completed" + }, + "tags": [] + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Epigraphe : status = optimal, objectif = 0.36064506\n", + " ecart a la forme naturelle : |x_epi - x_cvx|_2 = 3.55e-05\n", + "\n", + "Stationnarite : max |c - (mu+ - mu-)| = 1.20e-15\n", + "Complementarite : max |mu+ * (t - x)| = 6.94e-11\n", + "Complementarite : max |mu- * (t + x)| = 4.88e-11\n", + "\n", + "max |c_j| global = 0.012465 vs lambda = 0.012465\n", + "sur le support (|x_j|>0) : max | |c_j| - lambda | = 2.48e-05 (egalite des deux cotes)\n", + "hors du support (x_j=0) : max |c_j| = 0.012379 (marge 8.56e-05 sous le seuil)\n" + ] + } + ], + "source": [ + "# Forme epigraphe explicite : les multiplicateurs deviennent accessibles.\n", + "x_epi = cp.Variable(p)\n", + "t_epi = cp.Variable(p, nonneg=True)\n", + "c_pos = x_epi <= t_epi\n", + "c_neg = -x_epi <= t_epi\n", + "prob_epi = cp.Problem(\n", + " cp.Minimize(0.5 * cp.sum_squares(A @ x_epi - b) + lam * cp.sum(t_epi)),\n", + " [c_pos, c_neg])\n", + "prob_epi.solve(solver=cp.CLARABEL)\n", + "\n", + "x_epi_val = x_epi.value\n", + "t_epi_val = t_epi.value\n", + "mu_pos = np.asarray(c_pos.dual_value).ravel()\n", + "mu_neg = np.asarray(c_neg.dual_value).ravel()\n", + "corr = A.T @ (b - A @ x_epi_val)\n", + "\n", + "print(f'Epigraphe : status = {prob_epi.status}, objectif = {objectif(x_epi_val):.8f}')\n", + "print(f' ecart a la forme naturelle : |x_epi - x_cvx|_2 = {np.linalg.norm(x_epi_val - x_cvx):.2e}')\n", + "print()\n", + "print(f'Stationnarite : max |c - (mu+ - mu-)| = {np.max(np.abs(corr - (mu_pos - mu_neg))):.2e}')\n", + "print(f'Complementarite : max |mu+ * (t - x)| = {np.max(np.abs(mu_pos * (t_epi_val - x_epi_val))):.2e}')\n", + "print(f'Complementarite : max |mu- * (t + x)| = {np.max(np.abs(mu_neg * (t_epi_val + x_epi_val))):.2e}')\n", + "\n", + "nonz = np.abs(x_epi_val) > 1e-6\n", + "print()\n", + "print(f'max |c_j| global = {np.max(np.abs(corr)):.6f} vs lambda = {lam:.6f}')\n", + "print(f'sur le support (|x_j|>0) : max | |c_j| - lambda | = {np.max(np.abs(np.abs(corr[nonz]) - lam)):.2e}'\n", + " f' (egalite des deux cotes)')\n", + "print(f'hors du support (x_j=0) : max |c_j| = {np.max(np.abs(corr[~nonz])):.6f}'\n", + " f' (marge {lam - np.max(np.abs(corr[~nonz])):.2e} sous le seuil)')\n" + ] + }, + { + "cell_type": "markdown", + "id": "2f573320", + "metadata": { + "papermill": { + "duration": 0.007503, + "end_time": "2026-09-24T11:51:52.030894+00:00", + "exception": false, + "start_time": "2026-09-24T11:51:52.023391+00:00", + "status": "completed" + }, + "tags": [] + }, + "source": [ + "### Lecture — la théorie lue, pas recodée\n", + "\n", + "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.\n", + "\n", + "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.\n", + "\n", + "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.\n" + ] + }, + { + "cell_type": "markdown", + "id": "a6bcc6cc", + "metadata": { + "papermill": { + "duration": 0.007528, + "end_time": "2026-09-24T11:51:52.046245+00:00", + "exception": false, + "start_time": "2026-09-24T11:51:52.038717+00:00", + "status": "completed" + }, + "tags": [] + }, + "source": [ + "## 7. Un modèle, dix-huit résolutions : le chemin de régularisation et le seuil $\\lambda_{\\max}$\n", + "\n", + "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 :\n", + "\n", + "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 ;\n", + "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.\n", + "\n", + "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**.\n" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "9f7d3704", + "metadata": { + "execution": { + "iopub.execute_input": "2026-09-24T11:51:52.064475Z", + "iopub.status.busy": "2026-09-24T11:51:52.063954Z", + "iopub.status.idle": "2026-09-24T11:52:01.597068Z", + "shell.execute_reply": "2026-09-24T11:52:01.595054Z" + }, + "papermill": { + "duration": 9.544339, + "end_time": "2026-09-24T11:52:01.598413+00:00", + "exception": false, + "start_time": "2026-09-24T11:51:52.054074+00:00", + "status": "completed" + }, + "tags": [] + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "--- Chemin de regularisation : 12 resolutions du meme modele ---\n", + " lambda |x|_0 objectif temps\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + " 0.24929 30 6.234041 0.588s\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + " 0.14460 35 3.824345 0.563s\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + " 0.08387 39 2.297614 0.558s\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + " 0.04865 38 1.361812 0.569s\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + " 0.02822 48 0.801935 0.620s\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + " 0.01637 80 0.470996 0.600s\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + " 0.00949 110 0.275978 0.661s\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + " 0.00551 137 0.161340 0.807s\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + " 0.00319 165 0.094123 0.799s\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + " 0.00185 179 0.054802 0.684s\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + " 0.00107 186 0.031864 0.798s\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + " 0.00062 192 0.018509 0.677s\n", + "\n", + "--- Le seuil lambda_max, verifie par le solveur ---\n", + "lambda_max calcule (max_j |A_j^T b|) = 3.191173 | f(0) = 0.5||b||^2 = 24.915170\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "lambda = 0.999 * lambda_max : |x|_0(1e-6) = 1 | max|x| = 3.12e-03 | f(x) - f(0) = -4.98e-06\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "lambda = 1.000 * lambda_max : |x|_0(1e-6) = 1 | max|x| = 3.42e-04 | f(x) - f(0) = +9.94e-08\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "lambda = 1.001 * lambda_max : |x|_0(1e-6) = 1 | max|x| = 1.41e-06 | f(x) - f(0) = +6.39e-09\n" + ] + } + ], + "source": [ + "print('--- Chemin de regularisation : 12 resolutions du meme modele ---')\n", + "print(f'{\"lambda\":>9s} {\"|x|_0\":>6s} {\"objectif\":>11s} {\"temps\":>8s}')\n", + "for l in np.geomspace(lam * 20, lam / 20, 12):\n", + " v = cp.Variable(p)\n", + " pr = cp.Problem(cp.Minimize(0.5 * cp.sum_squares(A @ v - b) + l * cp.norm1(v)))\n", + " t0 = time.time()\n", + " pr.solve(solver=cp.CLARABEL)\n", + " t = time.time() - t0\n", + " xx = v.value\n", + " print(f'{l:9.5f} {int(np.sum(np.abs(xx) > 1e-6)):6d} {objectif(xx, l):11.6f} {t:7.3f}s')\n", + "\n", + "print()\n", + "print('--- Le seuil lambda_max, verifie par le solveur ---')\n", + "print(f'lambda_max calcule (max_j |A_j^T b|) = {lam_max:.6f} | f(0) = 0.5||b||^2 = '\n", + " f'{0.5 * np.linalg.norm(b) ** 2:.6f}')\n", + "f_zero = 0.5 * np.linalg.norm(b) ** 2\n", + "for fac in [0.999, 1.0, 1.001]:\n", + " l = fac * lam_max\n", + " v = cp.Variable(p)\n", + " pr = cp.Problem(cp.Minimize(0.5 * cp.sum_squares(A @ v - b) + l * cp.norm1(v)))\n", + " pr.solve(solver=cp.CLARABEL)\n", + " xx = v.value\n", + " print(f'lambda = {fac:.3f} * lambda_max : |x|_0(1e-6) = {int(np.sum(np.abs(xx) > 1e-6)):2d}'\n", + " f' | max|x| = {np.abs(xx).max():.2e} | f(x) - f(0) = {objectif(xx, l) - f_zero:+.2e}')\n" + ] + }, + { + "cell_type": "markdown", + "id": "8fae7af0", + "metadata": { + "papermill": { + "duration": 0.006813, + "end_time": "2026-09-24T11:52:01.614442+00:00", + "exception": false, + "start_time": "2026-09-24T11:52:01.607629+00:00", + "status": "completed" + }, + "tags": [] + }, + "source": [ + "### Lecture — le chemin coûte 12 résolutions, et le seuil se calcule\n", + "\n", + "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.\n", + "\n", + "Le seuil, lui, se lit proprement dans les trois lignes du tableau :\n", + "\n", + "- à $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 ;\n", + "- à $\\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 ;\n", + "- 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.\n", + "\n", + "La leçon de méthode est double. D'abord, un seuil se **calcule** ($\\|A^{\\top}b\\|_{\\infty}$, une ligne) 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.\n" + ] + }, + { + "cell_type": "markdown", + "id": "cc49308b", + "metadata": { + "papermill": { + "duration": 0.008888, + "end_time": "2026-09-24T11:52:01.631145+00:00", + "exception": false, + "start_time": "2026-09-24T11:52:01.622257+00:00", + "status": "completed" + }, + "tags": [] + }, + "source": [ + "## 8. Le problème du bloc A.2 : cvxpy comme arbitre de l'ADMM\n", + "\n", + "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.\n", + "\n", + "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é**.\n" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "eaa2ee45", + "metadata": { + "execution": { + "iopub.execute_input": "2026-09-24T11:52:01.645414Z", + "iopub.status.busy": "2026-09-24T11:52:01.644994Z", + "iopub.status.idle": "2026-09-24T11:52:01.822581Z", + "shell.execute_reply": "2026-09-24T11:52:01.821212Z" + }, + "papermill": { + "duration": 0.185305, + "end_time": "2026-09-24T11:52:01.823760+00:00", + "exception": false, + "start_time": "2026-09-24T11:52:01.638455+00:00", + "status": "completed" + }, + "tags": [] + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Probleme 2.11d : n=100, p=50, k=5, rho_corr=0.95, kappa(X)=8.18e+01\n", + "lambda_max_raw = 115.2719 | lambda_test = 0.5*lambda_max_raw = 57.6360\n", + "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\n", + "\n", + " rho it. converge f(z) f(z)-f* |z-b*|_2 temps\n", + " 1.0 500 False 160.305962 +1.42e+00 1.82e-01 0.058s\n", + " 5.0 500 False 158.893498 +3.04e-03 7.01e-03 0.060s\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + " 20.0 178 True 158.890705 +2.49e-04 2.08e-03 0.023s\n" + ] + } + ], + "source": [ + "def design_sparse_corr(n_=100, p_=50, k_=5, corr=0.95, snr=4.0, seed=20260914):\n", + " \"\"\"Design correle AR(1) — reprise exacte de 2.11d cellule 4.\"\"\"\n", + " local_rng = np.random.default_rng(seed)\n", + " idx = np.arange(p_)\n", + " Sigma = corr ** np.abs(idx[:, None] - idx[None, :])\n", + " L = np.linalg.cholesky(Sigma)\n", + " X = local_rng.standard_normal((n_, p_)) @ L.T\n", + " support_ = np.sort(local_rng.choice(p_, size=k_, replace=False))\n", + " beta_true_ = np.zeros(p_)\n", + " beta_true_[support_] = local_rng.standard_normal(k_)\n", + " sigma_ = np.linalg.norm(X @ beta_true_) / np.sqrt(n_) / snr\n", + " y_ = X @ beta_true_ + local_rng.standard_normal(n_) * sigma_\n", + " return X, beta_true_, y_, support_\n", + "\n", + "\n", + "def soft_threshold(z, lam_):\n", + " return np.sign(z) * np.maximum(np.abs(z) - lam_, 0.0)\n", + "\n", + "\n", + "def admm_lasso(X, y, lam_, rho, max_iter=500, eps_abs=1e-4, eps_rel=1e-4):\n", + " \"\"\"ADMM pour min 0.5||X b - y||^2 + lam||z||_1 s.c. b = z — reprise de 2.11d cellule 9.\"\"\"\n", + " n_, p_ = X.shape\n", + " XtX, Xty = X.T @ X, X.T @ y\n", + " L = np.linalg.cholesky(XtX + rho * np.eye(p_)) # pre-factorisation du b-update\n", + " beta = np.zeros(p_)\n", + " z = np.zeros(p_)\n", + " u = np.zeros(p_)\n", + " eps_p, eps_d = np.sqrt(p_) * eps_abs, np.sqrt(p_) * eps_abs\n", + " z_prev = z.copy()\n", + " converged = False\n", + " for it in range(1, max_iter + 1):\n", + " tmp = np.linalg.solve(L, Xty + rho * (z - u))\n", + " beta = np.linalg.solve(L.T, tmp)\n", + " z = soft_threshold(beta + u, lam_ / rho)\n", + " u = u + beta - z\n", + " r_p = np.linalg.norm(beta - z)\n", + " r_d = rho * np.linalg.norm(z - z_prev)\n", + " tol_p = eps_p + eps_rel * max(np.linalg.norm(beta), np.linalg.norm(z))\n", + " tol_d = eps_d + eps_rel * rho * np.linalg.norm(u)\n", + " if r_p < tol_p and r_d < tol_d:\n", + " converged = True\n", + " break\n", + " z_prev = z.copy()\n", + " return beta, z, it, converged\n", + "\n", + "\n", + "X, beta_true, y, supp = design_sparse_corr()\n", + "lam_max_raw = np.max(np.abs(X.T @ y))\n", + "lam_test = 0.5 * lam_max_raw # parametrage de 2.11d cellule 15\n", + "f_local = lambda bb: 0.5 * np.linalg.norm(X @ bb - y) ** 2 + lam_test * np.abs(bb).sum()\n", + "\n", + "v = cp.Variable(X.shape[1])\n", + "pr = cp.Problem(cp.Minimize(0.5 * cp.sum_squares(X @ v - y) + lam_test * cp.norm1(v)))\n", + "t0 = time.time()\n", + "pr.solve(solver=cp.CLARABEL)\n", + "t_cvx_11d = time.time() - t0\n", + "b_opt = v.value\n", + "f_star = f_local(b_opt)\n", + "\n", + "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}')\n", + "print(f'lambda_max_raw = {lam_max_raw:.4f} | lambda_test = 0.5*lambda_max_raw = {lam_test:.4f}')\n", + "print(f'Optimum certifie (CLARABEL) : f* = {f_star:.6f}, |b*|_0 = '\n", + " f'{int(np.sum(np.abs(b_opt) > 1e-6))}, support vrai = {list(supp)}, temps = {t_cvx_11d:.3f}s')\n", + "print()\n", + "print(f'{\"rho\":>6s} {\"it.\":>5s} {\"converge\":>9s} {\"f(z)\":>12s} {\"f(z)-f*\":>11s} {\"|z-b*|_2\":>10s} {\"temps\":>8s}')\n", + "for rho in [1.0, 5.0, 20.0]:\n", + " t0 = time.time()\n", + " bd, zd, itd, conv = admm_lasso(X, y, lam_test, rho)\n", + " tt = time.time() - t0\n", + " print(f'{rho:6.1f} {itd:5d} {str(conv):>9s} {f_local(zd):12.6f} '\n", + " f'{f_local(zd) - f_star:+11.2e} {np.linalg.norm(zd - b_opt):10.2e} {tt:7.3f}s')\n" + ] + }, + { + "cell_type": "markdown", + "id": "88eabb2b", + "metadata": { + "papermill": { + "duration": 0.009512, + "end_time": "2026-09-24T11:52:01.842301+00:00", + "exception": false, + "start_time": "2026-09-24T11:52:01.832789+00:00", + "status": "completed" + }, + "tags": [] + }, + "source": [ + "### Lecture — l'arbitre tranche, et le paramètre $\\rho$ réapparaît\n", + "\n", + "Le tableau se lit sur trois colonnes, et chacune raconte la même histoire sous un angle différent.\n", + "\n", + "**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.\n", + "\n", + "**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.\n", + "\n", + "**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.\n", + "\n", + "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.\n" + ] + }, + { + "cell_type": "markdown", + "id": "cd3f9bfb", + "metadata": { + "papermill": { + "duration": 0.008212, + "end_time": "2026-09-24T11:52:01.858982+00:00", + "exception": false, + "start_time": "2026-09-24T11:52:01.850770+00:00", + "status": "completed" + }, + "tags": [] + }, + "source": [ + "## 9. L'échelle : $n$ fixe, $p$ croissant\n", + "\n", + "`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$).\n", + "\n", + "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$.\n" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "39bff63f", + "metadata": { + "execution": { + "iopub.execute_input": "2026-09-24T11:52:01.875428Z", + "iopub.status.busy": "2026-09-24T11:52:01.874962Z", + "iopub.status.idle": "2026-09-24T11:52:19.806377Z", + "shell.execute_reply": "2026-09-24T11:52:19.804885Z" + }, + "papermill": { + "duration": 17.941664, + "end_time": "2026-09-24T11:52:19.807465+00:00", + "exception": false, + "start_time": "2026-09-24T11:52:01.865801+00:00", + "status": "completed" + }, + "tags": [] + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "n = 200 fixe, k = 30 fixe, lambda = 5 * Donoho-Johnstone(p) ; temps mesures en un passage\n", + " p cvxpy CLARABEL sklearn CD CD 50 ep. ISTA 300 it.\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + " 500 0.614s (12 it.) 0.008s 0.315s 0.029s\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + " 2000 3.687s (16 it.) 0.073s 1.134s 0.101s\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + " 4000 8.898s (18 it.) 0.315s 2.562s 0.152s\n" + ] + } + ], + "source": [ + "print('n = 200 fixe, k = 30 fixe, lambda = 5 * Donoho-Johnstone(p) ; temps mesures en un passage')\n", + "print(f'{\"p\":>6s} {\"cvxpy CLARABEL\":>15s} {\"sklearn CD\":>12s} {\"CD 50 ep.\":>10s} {\"ISTA 300 it.\":>12s}')\n", + "for pp in [500, 2000, 4000]:\n", + " r2 = np.random.default_rng(1000 + pp)\n", + " A_ = r2.standard_normal((n, pp)) / np.sqrt(n)\n", + " s_ = r2.choice(pp, size=k, replace=False)\n", + " xt_ = np.zeros(pp)\n", + " xt_[s_] = r2.standard_normal(k)\n", + " b_ = A_ @ xt_ + noise_level * r2.standard_normal(n)\n", + " lam_ = 5.0 * sigma * np.sqrt(2.0 * np.log(pp)) / np.sqrt(n)\n", + "\n", + " v = cp.Variable(pp)\n", + " pr = cp.Problem(cp.Minimize(0.5 * cp.sum_squares(A_ @ v - b_) + lam_ * cp.norm1(v)))\n", + " t0 = time.time()\n", + " pr.solve(solver=cp.CLARABEL)\n", + " t_cvx_ = time.time() - t0\n", + " it_cvx_ = pr.solver_stats.num_iters\n", + "\n", + " t0 = time.time()\n", + " Lasso(alpha=lam_ / n, fit_intercept=False, max_iter=5000, tol=1e-6).fit(A_, b_)\n", + " t_skl_ = time.time() - t0\n", + "\n", + " t0 = time.time()\n", + " lasso_cd(A_, b_, lam_, max_iter=50)\n", + " t_cd_ = time.time() - t0\n", + "\n", + " t0 = time.time()\n", + " L_ = np.linalg.norm(A_, 2) ** 2\n", + " xi = np.zeros(pp)\n", + " for _ in range(300):\n", + " xi = prox_l1(xi - A_.T @ (A_ @ xi - b_) / L_, lam_ / L_)\n", + " t_ista_ = time.time() - t0\n", + "\n", + " 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" + ] + }, + { + "cell_type": "markdown", + "id": "eabb250c", + "metadata": { + "papermill": { + "duration": 0.006777, + "end_time": "2026-09-24T11:52:19.822093+00:00", + "exception": false, + "start_time": "2026-09-24T11:52:19.815316+00:00", + "status": "completed" + }, + "tags": [] + }, + "source": [ + "### Lecture — le déclaratif paie sa généralité à l'échelle\n", + "\n", + "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$.\n", + "\n", + "La raison est structurelle, pas accidentelle. Trois mécanismes se cumulent :\n", + "\n", + "- **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 ;\n", + "- **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 ;\n", + "- **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.\n", + "\n", + "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.\n", + "\n", + "Tableau de synthèse du bloc B (acceptance $\\S8$ de #16061) :\n", + "\n", + "| Axe | From scratch (2.11b/2.11c/2.11d) | sklearn (2.11c) | cvxpy (2.11e) |\n", + "|---|---|---|---|\n", + "| Objectif final (même $\\lambda$) | identique | identique | identique |\n", + "| Temps ($p=500$) | ~0,1-0,5 s | ~0,01 s | ~0,6 s |\n", + "| Temps ($p=4000$) | ~0,2-2 s | ~0,3 s | plusieurs s |\n", + "| Itérations | 82 epochs / 2000 / 500 | warm start | ~12-18 |\n", + "| Lignes de code | 15-25 | 1 | 2 |\n", + "| Coût du changement de problème | un nouvel algorithme | non couvert | une ligne |\n", + "| Garantie de convexité | analyse manuelle | non | ruleset DCP |\n", + "| Certificat dual | codé à la main | non exposé | lu dans les duales |\n" + ] + }, + { + "cell_type": "markdown", + "id": "65e52f07", + "metadata": { + "papermill": { + "duration": 0.008714, + "end_time": "2026-09-24T11:52:19.839185+00:00", + "exception": false, + "start_time": "2026-09-24T11:52:19.830471+00:00", + "status": "completed" + }, + "tags": [] + }, + "source": [ + "## 10. Synthèse\n", + "\n", + "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$.\n", + "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.\n", + "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**.\n", + "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.\n", + "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.\n", + "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.\n" + ] + }, + { + "cell_type": "markdown", + "id": "c45777f4", + "metadata": { + "papermill": { + "duration": 0.008602, + "end_time": "2026-09-24T11:52:19.855738+00:00", + "exception": false, + "start_time": "2026-09-24T11:52:19.847136+00:00", + "status": "completed" + }, + "tags": [] + }, + "source": [ + "## 11. Exercices\n", + "\n", + "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.\n" + ] + }, + { + "cell_type": "markdown", + "id": "0c56f325", + "metadata": { + "papermill": { + "duration": 0.007682, + "end_time": "2026-09-24T11:52:19.871435+00:00", + "exception": false, + "start_time": "2026-09-24T11:52:19.863753+00:00", + "status": "completed" + }, + "tags": [] + }, + "source": [ + "### Exercice 1 — group Lasso : le déclaratif absorbe la structure\n", + "\n", + "Le **group Lasso** pénalise des **groupes** de coordonnées entiers :\n", + "\n", + "$$\\min_x \\; \\tfrac12 \\|Ax - b\\|_2^2 + \\lambda \\sum_{g=1}^{G} \\|x_g\\|_2$$\n", + "\n", + "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.\n", + "\n", + "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.\n", + "\n", + "```python\n", + "# Indice : partitionner les p coordonnees en G groupes contigus de taille p//G\n", + "# Indice : la penalite s'ecrit cp.sum([cp.norm(x[groupe], 2) for groupe in groupes])\n", + "# Indice : compter les groupes ANNULES, pas les coordonnees (c'est la bonne granularite ici)\n", + "```\n" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "b1fcbdb2", + "metadata": { + "execution": { + "iopub.execute_input": "2026-09-24T11:52:19.891566Z", + "iopub.status.busy": "2026-09-24T11:52:19.891021Z", + "iopub.status.idle": "2026-09-24T11:52:19.898025Z", + "shell.execute_reply": "2026-09-24T11:52:19.896649Z" + }, + "papermill": { + "duration": 0.01834, + "end_time": "2026-09-24T11:52:19.899294+00:00", + "exception": false, + "start_time": "2026-09-24T11:52:19.880954+00:00", + "status": "completed" + }, + "tags": [] + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Exercice 1 a completer en TP.\n" + ] + } + ], + "source": [ + "def exercice_1_group_lasso(A, b, lam, G=50, solver=cp.CLARABEL):\n", + " \"\"\"Exercice 1 : group Lasso declaratif — min 0.5||Ax-b||^2 + lam*sum_g ||x_g||_2.\n", + "\n", + " Retourne (x, nb_groupes_annules, nb_coordonnees_non_nulles), ou None si non complete.\n", + " \"\"\"\n", + " p_ = A.shape[1]\n", + " x = cp.Variable(p_)\n", + " # TODO etudiant : construire la partition en G groupes contigus de taille p//G\n", + " # TODO etudiant : ecrire la penalite somme des normes 2 par groupe\n", + " # TODO etudiant : resoudre, puis compter les groupes dont la norme est sous 1e-6\n", + " return None # TODO etudiant\n", + "\n", + "# Etape 1 : comparer au Lasso sur le meme lambda — le group Lasso annule-t-il par blocs ?\n", + "x_gr, n_gr, n_co = (None, None, None)\n", + "print('Exercice 1 a completer en TP.')\n" + ] + }, + { + "cell_type": "markdown", + "id": "83e5dcfc", + "metadata": { + "papermill": { + "duration": 0.009362, + "end_time": "2026-09-24T11:52:19.917943+00:00", + "exception": false, + "start_time": "2026-09-24T11:52:19.908581+00:00", + "status": "completed" + }, + "tags": [] + }, + "source": [ + "### Exercice 2 — la forme contrainte, et le lien dual avec la forme pénalisée\n", + "\n", + "Le Lasso a deux écritures équivalentes, reliées par la dualité lagrangienne :\n", + "\n", + "$$\\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$$\n", + "\n", + "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$.\n", + "\n", + "```python\n", + "# Indice : le multiplicateur dual de la contrainte s'obtient par contrainte.dual_value\n", + "# Indice : balayer tau autour de ||x_cvx||_1 et verifier que la solution ne bouge plus\n", + "# Indice : comparer le dual_value de la contrainte a lam (c'est le sens du theoreme)\n", + "```\n" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "7b7b5373", + "metadata": { + "execution": { + "iopub.execute_input": "2026-09-24T11:52:19.939060Z", + "iopub.status.busy": "2026-09-24T11:52:19.938560Z", + "iopub.status.idle": "2026-09-24T11:52:19.945755Z", + "shell.execute_reply": "2026-09-24T11:52:19.944374Z" + }, + "papermill": { + "duration": 0.019588, + "end_time": "2026-09-24T11:52:19.947081+00:00", + "exception": false, + "start_time": "2026-09-24T11:52:19.927493+00:00", + "status": "completed" + }, + "tags": [] + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Exercice 2 a completer en TP.\n" + ] + } + ], + "source": [ + "def exercice_2_forme_contrainte(A, b, lam, x_star, tol=1e-3):\n", + " \"\"\"Exercice 2 : forme contrainte min 0.5||Ax-b||^2 s.c. ||x||_1 <= tau.\n", + "\n", + " A partir de x_star (solution penalisee), retrouver tau = ||x_star||_1, resoudre la\n", + " forme contrainte, comparer les solutions et lire le multiplicateur dual.\n", + " Retourne un dict ou None si non complete.\n", + " \"\"\"\n", + " p_ = A.shape[1]\n", + " x = cp.Variable(p_)\n", + " tau = np.abs(x_star).sum() if x_star is not None else None\n", + " # TODO etudiant : ecrire la contrainte cp.norm1(x) <= tau\n", + " # TODO etudiant : resoudre et mesurer |x_contraint - x_star|_2\n", + " # TODO etudiant : lire contrainte.dual_value et le comparer a lam\n", + " return None # TODO etudiant\n", + "\n", + "res2 = exercice_2_forme_contrainte(A, b, lam, x_cvx)\n", + "print('Exercice 2 a completer en TP.')\n" + ] + }, + { + "cell_type": "markdown", + "id": "809a7574", + "metadata": { + "papermill": { + "duration": 0.009182, + "end_time": "2026-09-24T11:52:19.964139+00:00", + "exception": false, + "start_time": "2026-09-24T11:52:19.954957+00:00", + "status": "completed" + }, + "tags": [] + }, + "source": [ + "### Exercice 3 — $\\lambda_{\\max}$ et le début du chemin, par le solveur seul\n", + "\n", + "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.\n", + "\n", + "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.\n", + "\n", + "```python\n", + "# Indice : grille np.geomspace(1e-2, 10.0, 25) ; pour chaque lambda, resoudre et tester max|x| < 1e-6\n", + "# Indice : comparer le plus petit lambda annulant de la grille a lam_max (lequel est au-dessus ?)\n", + "# Indice : chronometrer la boucle et comparer a une seule evaluation de lam_max\n", + "```\n" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "id": "04b29fb6", + "metadata": { + "execution": { + "iopub.execute_input": "2026-09-24T11:52:19.986412Z", + "iopub.status.busy": "2026-09-24T11:52:19.985827Z", + "iopub.status.idle": "2026-09-24T11:52:19.993984Z", + "shell.execute_reply": "2026-09-24T11:52:19.992624Z" + }, + "papermill": { + "duration": 0.021342, + "end_time": "2026-09-24T11:52:19.995193+00:00", + "exception": false, + "start_time": "2026-09-24T11:52:19.973851+00:00", + "status": "completed" + }, + "tags": [] + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Exercice 3 a completer en TP.\n" + ] + } + ], + "source": [ + "def exercice_3_lambda_max_par_solveur(A, b, n_points=25):\n", + " \"\"\"Exercice 3 : trouver par le solveur le plus petit lambda d'une grille annulant x.\n", + "\n", + " Retourne (lambda_trouve, lambda_max_analytique, temps_boucle) ou None si non complete.\n", + " \"\"\"\n", + " p_ = A.shape[1]\n", + " grille = np.geomspace(1e-2, 10.0, n_points)\n", + " # TODO etudiant : boucler sur la grille, resoudre, tester max|x| < 1e-6\n", + " # TODO etudiant : comparer a np.max(np.abs(A.T @ b)) et au cout d'une seule evaluation\n", + " return None # TODO etudiant\n", + "\n", + "res3 = exercice_3_lambda_max_par_solveur(A, b)\n", + "print('Exercice 3 a completer en TP.')\n" + ] + }, + { + "cell_type": "markdown", + "id": "2b13d4bb", + "metadata": { + "papermill": { + "duration": 0.008195, + "end_time": "2026-09-24T11:52:20.012640+00:00", + "exception": false, + "start_time": "2026-09-24T11:52:20.004445+00:00", + "status": "completed" + }, + "tags": [] + }, + "source": [ + "## 12. Références\n", + "\n", + "- 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.\n", + "- Grant, M., Boyd, S., Ye, Y. (2006). *Disciplined Convex Programming*. In *Global Optimization: From Theory to Implementation*. — la certification de convexité à la compilation.\n", + "- 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$.\n", + "- 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`).\n", + "- 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$.\n", + "- 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.\n", + "- 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$.\n", + "- 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.\n", + "- 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`.\n", + "- Documentation `cvxpy` : `Problem.get_problem_data`, `SolverStats`, `Constraint.dual_value` — les trois API utilisées aux $\\S3$, $\\S4$ et $\\S6$.\n" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.13.7" + }, + "papermill": { + "default_parameters": {}, + "duration": 56.473014, + "end_time": "2026-09-24T11:52:20.803716+00:00", + "environment_variables": {}, + "exception": null, + "input_path": "2.11e-CVXPY-Optimisation.ipynb", + "output_path": "2.11e-CVXPY-Optimisation.ipynb", + "parameters": {}, + "start_time": "2026-09-24T11:51:24.330702+00:00", + "version": "2.6.0" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} \ No newline at end of file