Modèle linéaire, moindres carrés, régularisation

La page précédente a montré que le MLE, sous bruit Gaussien blanc, minimise une erreur quadratique, sans jamais expliciter sa solution quand le modèle est linéaire en $\mathit{\theta }$. Cette page ferme ce cas en quatre temps : le modèle linéaire et sa solution analytique, la généralisation aux fonctions de base, la facilité trompeuse des modèles trop riches, et la régularisation qui les dompte.

1. Le modèle linéaire

Reprenons la situation qui ouvre le chapitre : un voltmètre bruité mesure $N=10$ fois une tension continue inconnue ${\theta }_0$. La $n$-ième mesure s'écrit $x_n={\theta }_0+w_n$, où $w_n$ est le bruit de mesure. Si l'alimentation dérive lentement au fil du temps, le modèle s'enrichit d'une pente : $x_n={\theta }_0+{\theta }_{1}n+w_n$, la sortie glisse linéairement avec l'instant de mesure $n$. Ce second modèle est la régression linéaire : une droite ajustée à des points bruités. Ces deux modèles partagent une propriété commune : la sortie dépend linéairement des paramètres ${\theta }_l$.

Forme scalaire, puis forme matricielle

Sous forme scalaire, le modèle linéaire général s'écrit :

$$x_n=\sum _{l=0}^{p-1}{\theta }_l{h}_l[n]+w_n,n=1,\dots ,N,$$

où :

• $x_n$ est la $n$-ième sortie mesurée,
• $h_0[n],\dots ,h_{p-1}[n]$ sont les composantes du modèle : des séquences connues, $h_l[n]$ désignant la $n$-ième valeur de la composante $l$,
• ${\theta }_0,\dots ,{\theta }_{p-1}$ sont les paramètres inconnus à estimer,
• $w_n$ est le bruit de mesure.

Ce modèle marque le passage au multidimensionnel : non plus un $\theta$ scalaire comme dans les pages précédentes, mais $p$ paramètres ${\theta }_l$, chacun pondérant une composante connue de la sortie. La tension continue en est le cas $p=1$ avec $h_0[n]=1$ ; la dérive, le cas $p=2$ avec $h_0[n]=1$ et $h_1[n]=n$.

Empiler les $N$ équations scalaires permet d'obtenir une écriture compacte.

Modèle linéaire

$$\mathbf{x}=\mathbf{H}\mathit{\theta }+\mathbf{w},\mathbf{w}\sim \mathcal{N}(\mathbf{0},{\sigma }^2\mathbf{I}),$$

où :

• $\mathbf{x}=(x_1,\dots ,x_N{)}^T$ est le vecteur des sorties,
• $\mathbf{H}$ (dimension $N\times p$, connue, de rang plein, $N>p$) porte le modèle : sa colonne $l$ empile les $N$ valeurs de la composante $l$, $(h_l[1],\dots ,h_l[N]{)}^T$,
• $\mathit{\theta }\in {\mathbb{R}}^p$ est le vecteur des paramètres.

Exemples

Le modèle linéaire est fréquemment rencontré en ingénierie. Les exemples suivants montrent la polyvalence de ce modèle.

Estimation d'une composante continue

• Forme scalaire :

$$x_n={\theta }_0+w_n$$

• Forme matricielle :

$$\left(\begin{array}{c}{x}_1\\ x_2\\ ⋮\\ x_N\end{array}\right)=\left(\begin{array}{c}1\\ 1\\ ⋮\\ 1\end{array}\right){\theta }_0+\mathbf{w}$$

Régression linéaire

• Forme scalaire :

$$x_n={\theta }_0+{\theta }_{1}n+w_n$$

• Forme matricielle :

$$\left(\begin{array}{c}{x}_1\\ x_2\\ x_3\\ ⋮\\ x_N\end{array}\right)=\left(\begin{array}{cc}1& 1\\ 1& 2\\ 1& 3\\ ⋮& ⋮\\ 1& N\end{array}\right)\left(\begin{array}{c}{\theta }_0\\ {\theta }_1\end{array}\right)+\mathbf{w}$$

Régression polynomiale d'ordre $L$

• Forme scalaire :

$$x_n=\sum _{l=0}^L{\theta }_l{n}^l+w_n$$

• Forme matricielle, illustrée pour $L=4$ ($p=5$ colonnes, une par puissance de l'instant, de $n^0$ à $n^4$) ; le passage à un $L$ quelconque ne fait qu'ajouter des colonnes $n^l$ de la même manière :

$$\left(\begin{array}{c}{x}_1\\ x_2\\ x_3\\ ⋮\\ x_N\end{array}\right)=\left(\begin{array}{ccccc}1& 1& 1& 1& 1\\ 1& 2& 2^2& 2^3& 2^4\\ 1& 3& 3^2& 3^3& 3^4\\ ⋮& ⋮& ⋮& ⋮& ⋮\\ 1& N& N^2& N^3& N^4\end{array}\right)\left(\begin{array}{c}{\theta }_0\\ {\theta }_1\\ {\theta }_2\\ {\theta }_3\\ {\theta }_4\end{array}\right)+\mathbf{w}$$

Cette matrice, dont chaque colonne élève les instants à une puissance croissante, se nomme matrice de Vandermonde (np.vander en NumPy). La régression polynomiale est l'illustration principale de cette page : elle reviendra à chaque section.

Estimation des coefficients d'un filtre FIR

• Forme scalaire, avec $s$ un signal d'entrée connu et la convention $s_k=0$ pour $k\le 0$ :

$$x_n=\sum _{l=0}^L{\theta }_l{s}_{n-l}+w_n$$

• Forme matricielle :

$$\left(\begin{array}{c}{x}_1\\ x_2\\ ⋮\\ x_N\end{array}\right)=\left(\begin{array}{cccc}{s}_1& 0& \cdots & 0\\ s_2& s_1& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ s_N& s_{N-1}& \cdots & s_{N-L}\end{array}\right)\left(\begin{array}{c}{\theta }_0\\ {\theta }_1\\ ⋮\\ {\theta }_L\end{array}\right)+\mathbf{w}$$

Cette matrice $\mathbf{H}$ est construite à partir du signal connu $s$ (l'entrée du filtre).

2. Estimation des paramètres

Hypothèses

Pour estimer $\mathit{\theta }$, nous allons considérer que le signal peut être décrit par le modèle suivant :

$$\mathbf{x}=\mathbf{H}\mathit{\theta }+\mathbf{w}$$

Nous retenons l'estimateur des moindres carrés qui, lorsque le bruit est Gaussien i.i.d., est équivalent à l'estimateur MLE.

Expression

Proposition, Solution des moindres carrés

$$\hat{\mathit{\theta }}={\mathbf{H}}^{+}\mathbf{x}=({\mathbf{H}}^T\mathbf{H}{)}^{-1}{\mathbf{H}}^T\mathbf{x},$$

où :

• ${\mathbf{H}}^{+}=({\mathbf{H}}^T\mathbf{H}{)}^{-1}{\mathbf{H}}^T$ est la pseudo-inverse de $\mathbf{H}$.

Démonstration

Rappel de dérivation matricielle. Pour $\mathbf{a}\in {\mathbb{R}}^p$ constant et $\mathbf{A}$ symétrique constante,

$$\begin{array}{rl}{\mathrm{\nabla }}_{\mathit{\theta }}({\mathbf{a}}^T\mathit{\theta })& =\mathbf{a},\\ {\mathrm{\nabla }}_{\mathit{\theta }}({\mathit{\theta }}^T\mathbf{A}\mathit{\theta })& =2\mathbf{A}\mathit{\theta }.\end{array}$$

La solution. En développant $J(\mathit{\theta })=\Vert \mathbf{x}-\mathbf{H}\mathit{\theta }{\Vert }^2$ :

$$J(\mathit{\theta })={\mathbf{x}}^T\mathbf{x}-2{\mathit{\theta }}^T{\mathbf{H}}^T\mathbf{x}+{\mathit{\theta }}^T{\mathbf{H}}^T\mathbf{H}\mathit{\theta }.$$

Annuler le gradient donne

$$\mathrm{\nabla }J(\mathit{\theta })=-2{\mathbf{H}}^T\mathbf{x}+2{\mathbf{H}}^T\mathbf{H}\mathit{\theta }=\mathbf{0},$$

soit les équations normales ${\mathbf{H}}^T\mathbf{H}\mathit{\theta }={\mathbf{H}}^T\mathbf{x}$, d'où la solution $\hat{\mathit{\theta }}=({\mathbf{H}}^T\mathbf{H}{)}^{-1}{\mathbf{H}}^T\mathbf{x}$.

${\mathbf{H}}^T\mathbf{H}$ est inversible car $\mathbf{H}$ est de rang plein et $N>p$ : le système comporte davantage d'équations que d'inconnues.

Notons qu'en pratique l'inverse $({\mathbf{H}}^T\mathbf{H}{)}^{-1}$ ne doit pas être calculée explicitement : l'inversion est coûteuse et potentiellement instable numériquement. Il est préférable d'utiliser les fonctions optimisées np.linalg.lstsq(H, x) ou np.linalg.solve(H.T H, H.T x).

Propriétés

Propriétés

$\hat{\mathit{\theta }}$ est sans biais, de covariance ${\sigma }^2({\mathbf{H}}^T\mathbf{H}{)}^{-1}$. Cette covariance atteint la borne CRLB : $\hat{\mathit{\theta }}$ est efficace.

Démonstration

1. Sans biais, covariance ${\sigma }^2({\mathbf{H}}^T\mathbf{H}{)}^{-1}$. En substituant $\mathbf{x}=\mathbf{H}\mathit{\theta }+\mathbf{w}$ :

$$\hat{\mathit{\theta }}=({\mathbf{H}}^T\mathbf{H}{)}^{-1}{\mathbf{H}}^T(\mathbf{H}\mathit{\theta }+\mathbf{w})=\mathit{\theta }+({\mathbf{H}}^T\mathbf{H}{)}^{-1}{\mathbf{H}}^T\mathbf{w}.$$

$E[\hat{\mathit{\theta }}]=\mathit{\theta }$ (le second terme est d'espérance nulle), pour tout $\mathit{\theta }$ : l'estimateur est sans biais. Sa covariance vaut

$$\mathrm{cov}(\hat{\mathit{\theta }})=({\mathbf{H}}^T\mathbf{H}{)}^{-1}{\mathbf{H}}^T\mathrm{cov}(\mathbf{w})\mathbf{H}({\mathbf{H}}^T\mathbf{H}{)}^{-1}={\sigma }^2({\mathbf{H}}^T\mathbf{H}{)}^{-1}.$$

2. Atteint la CRLB. La log-vraisemblance est $\ell (\mathit{\theta })=-\frac{N}{2}\ln (2\pi {\sigma }^2)-\frac{1}{2{\sigma }^2}\Vert \mathbf{x}-\mathbf{H}\mathit{\theta }{\Vert }^2$, de hessienne ${\mathrm{\nabla }}_{\mathit{\theta }}^2\ell =-\frac{1}{{\sigma }^2}{\mathbf{H}}^T\mathbf{H}$. La matrice d'information de Fisher (page Fondations) est donc $\mathbf{I}(\mathit{\theta })={\mathbf{H}}^T\mathbf{H}/{\sigma }^2$, constante : la CRLB correspondante, $\mathbf{I}(\mathit{\theta }{)}^{-1}={\sigma }^2({\mathbf{H}}^T\mathbf{H}{)}^{-1}$, est exactement la covariance établie au point 1. $\hat{\mathit{\theta }}$ atteint la borne : c'est un estimateur efficace, de variance minimale parmi les estimateurs sans biais.

3. C'est la solution des moindres carrés. Un ingénieur cherchant, sans formuler aucune hypothèse Gaussienne, le $\mathit{\theta }$ qui minimise $\Vert \mathbf{x}-\mathbf{H}\mathit{\theta }{\Vert }^2$ obtiendrait par calcul direct exactement la même expression $({\mathbf{H}}^T\mathbf{H}{)}^{-1}{\mathbf{H}}^T\mathbf{x}$.

Interprétation géométrique

Le produit $\mathbf{H}\hat{\mathit{\theta }}$ admet une interprétation géométrique : il correspond à la projection orthogonale de $\mathbf{x}$ sur le sous-espace engendré par les colonnes de $\mathbf{H}$ ($\mathrm{span}(\mathbf{H})$).

$$\mathbf{H}\hat{\mathit{\theta }}=\mathbf{P}\mathbf{x},$$

où :

• $\mathbf{P}=\mathbf{H}({\mathbf{H}}^T\mathbf{H}{)}^{-1}{\mathbf{H}}^T$ est le projecteur orthogonal sur $\mathrm{span}(\mathbf{H})$ (${\mathbf{P}}^2=\mathbf{P}$, ${\mathbf{P}}^T=\mathbf{P}$).

Géométriquement, le résidu $\mathbf{x}-\mathbf{H}\hat{\mathit{\theta }}=(\mathbf{I}-\mathbf{P})\mathbf{x}$ est orthogonal à $\mathrm{span}(\mathbf{H})$ : ${\mathbf{H}}^T(\mathbf{x}-\mathbf{H}\hat{\mathit{\theta }})={\mathbf{H}}^T\mathbf{x}-{\mathbf{H}}^T\mathbf{H}\hat{\mathit{\theta }}=\mathbf{0}$. Si cette lecture est juste, le vecteur des sorties, sa projection et le résidu doivent former un triangle rectangle, l'angle droit posé sur le sous-espace.

Figure 2.6, Géométrie des moindres carrés : $\mathbf{x}$, sa projection $\mathbf{H}\hat{\mathit{\theta }}$ sur $\mathrm{span}(\mathbf{H})$, et le résidu $\mathbf{x}-\mathbf{H}\hat{\mathit{\theta }}$ (script figures/ch2/03-projection-geometrie.py).

Exemple : régression polynomiale d'ordre 4

Reprenons la régression polynomiale du §1, avec $L=4$. Une sortie $x_n$ suit un polynôme d'ordre 4 de l'instant $n$, noyé dans un bruit additif Gaussien ($\sigma =1$) ; un polynôme du même ordre est ajusté par moindres carrés sur $N=25$ mesures.

Figure 2.7, Régression polynomiale d'ordre 4 ajustée par moindres carrés sur $N=25$ mesures, bruit Gaussien $\sigma =1$ (script figures/ch2/03-regression-polynomiale-exemple.py).

Malgré un bruit d'écart-type 1, la régression retrouve fidèlement l'allure du signal : les coefficients estimés, $\hat{\mathit{\theta }}\approx (3,05,0,33,-\mathrm{0,054},\mathrm{0,001},\mathrm{0,000}2)$, restent proches des coefficients vrais ${\mathit{\theta }}^{\star }=(3,0,4,-0,05,0,\mathrm{0,000}2)$.

from scratch

import numpy as np
rng = np.random.default_rng(2026)

N, L = 25, 4
n = np.arange(N) - (N - 1) / 2   # instants centrés : conditionnement de H bien meilleur
H = np.column_stack([n**l for l in range(L + 1)])

theta_star = np.array([3.0, 0.4, -0.05, 0.0, 0.0002])
x = H  theta_star + rng.normal(0.0, 1.0, N)

theta_hat, *_ = np.linalg.lstsq(H, x, rcond=None)
theta_hat
# >>> array([ 3.0479,  0.3298, -0.0543,  0.0008,  0.0002])

numpy

H = np.vander(n, L + 1, increasing=True)   # équivalent à la boucle ci-dessus, plus concis
theta_hat, *_ = np.linalg.lstsq(H, x, rcond=None)

3. Le modèle sur fonctions de base

En dehors de quelques cas simples, une sortie réelle s'écrit rarement directement sous forme linéaire. Une solution consiste à projeter d'abord l'entrée par une fonction non linéaire connue, puis à rester linéaire dans cet espace projeté. Cette technique augmente la capacité du modèle sans quitter le cadre des moindres carrés.

Notons $t_n$ la grandeur en laquelle les composantes sont évaluées et écrivons la composante $h_l[n]$ sous la forme $h_l[n]={\phi }_l(t_n)$ où ${\phi }_l$ est une fonction connue et $t_n$ l'entrée associée à la $n$-ième mesure. Pour les signaux du §1, l'entrée est l'instant, $t_n=n$ ; mais rien n'impose ni que l'entrée soit le temps, ni que les fonctions soient des puissances. L'entrée peut même être multidimensionnelle, ${\mathbf{t}}_n\in {\mathbb{R}}^d$, notée en gras lorsqu'elle est vectorielle. Les données forment alors des couples $(t_n,x_n)$, une entrée et la sortie mesurée correspondante, exactement la structure des jeux de données du machine learning.

Hypothèses

Les échantillons reçus peuvent se décomposer sous la forme :

$$x_n=\sum _{j=0}^{M-1}{\theta }_j{\phi }_j(t_n)+w_n.$$

Sous forme matricielle, le modèle reste linéaire :

$$\mathbf{x}=\mathbf{H}\mathit{\theta }+\mathbf{w},$$

où :

• ${\phi }_0,\dots ,{\phi }_{M-1}$ sont $M$ fonctions connues de l'entrée, les fonctions de base (linear basis functions),
• $\mathbf{H}=[{\mathit{\phi }}_0\cdots {\mathit{\phi }}_{M-1}]$, avec ${\mathit{\phi }}_j=({\phi }_j(t_1),\dots ,{\phi }_j(t_N){)}^T$ la colonne $j$.

Chaque ${\phi }_j$ projette l'entrée sur un axe d'un nouvel espace de représentation, et c'est dans cet espace projeté, non dans l'espace d'origine, que la régression redevient linéaire. En rassemblant les projections dans un vecteur, $\mathit{\phi }(\mathbf{t})=({\phi }_0(\mathbf{t}),\dots ,{\phi }_{M-1}(\mathbf{t}){)}^T$, le modèle s'écrit d'un trait, $x_n=\mathit{\phi }({\mathbf{t}}_n{)}^T\mathit{\theta }+w_n$, et la chaîne complète se résume en un schéma.

Figure 2.8, Chaîne du modèle sur fonctions de base : l'entrée ${\mathbf{t}}_n$ est projetée en dimension $M$ par $\mathit{\phi }$, le produit scalaire avec $\mathit{\theta }$ pondère et somme ces projections, et l'ajout du bruit $w_n$ produit la sortie $x_n$ (source figures/ch2/03-schema-fonctions-base.tex).

Exemples : trois familles de bases

À titre d'exemple, les modèles suivants montrent différents choix de fonctions de base.

Base polynomiale

$$\mathit{\phi }(t_n)=\left(\begin{array}{c}1\\ t_n\\ ⋮\\ t_n^{M-1}\end{array}\right)$$

• la composante continue, la droite et la régression polynomiale du §1 en sont des cas particuliers (avec $t_n=n$).

Base de Fourier

Les fonctions trigonométriques forment la deuxième famille classique. Le cas le plus simple est l'estimation des amplitudes d'une somme de sinusoïdes de pulsations ${\omega }_k$ connues :

$$\mathit{\phi }(t_n)=\left(\begin{array}{c}\cos ({\omega }_1{t}_n)\\ ⋮\\ \cos ({\omega }_{M/2}{t}_n)\\ \sin ({\omega }_1{t}_n)\\ ⋮\\ \sin ({\omega }_{M/2}{t}_n)\end{array}\right)$$

Base gaussienne

$${\phi }_j(t)=\exp (-\frac{(t-c_j{)}^2}{2{\ell }^2}),$$

où :

• $c_j$ est le centre de la $j$-ième fonction de base,
• $\ell$ est leur largeur, commune.

Par exemple, $M=9$ centres régulièrement espacés sur $[-1,1]$ ($c_j=-1+0,25j$) et une largeur $\ell =0,25$ pavent l'intervalle : chaque ${\phi }_j$ ne « voit » que les entrées proches de son centre. Contrairement aux bases polynomiale et de Fourier, dont chaque fonction agit sur tout l'axe, la base gaussienne est locale : un coefficient ${\theta }_j$ ne façonne la courbe qu'au voisinage de $c_j$. Trois représentants de chaque famille, tracés côte à côte, rendent cette différence visible.

Figure 2.9, Trois représentants de chaque famille de fonctions de base sur $[-1,1]$ : polynomiale ($x$, $x^2$, $x^3$), Fourier ($\cos \pi x$, $\sin \pi x$, $\cos 2\pi x$), gaussienne (centres $-0,75$, $0$, $0,75$, largeur $\ell =0,25$) (script figures/ch2/03-fonctions-de-base.py).

Expression

Proposition, Solution des moindres carrés

$$\hat{\mathit{\theta }}={\mathbf{H}}^{+}\mathbf{x}={\mathbf{G}}^{-1}{\mathbf{H}}^T\mathbf{x},$$

où :

• $\mathbf{G}={\mathbf{H}}^T\mathbf{H}$ est la matrice de Gram des fonctions de base, de dimension $M\times M$. Son terme $(j,j^{\prime })$ est le produit scalaire de deux fonctions de base échantillonnées sur les entrées :

$$G_{j,j^{\prime }}={\mathit{\phi }}_j^T{\mathit{\phi }}_{j^{\prime }}=\sum _{n=1}^N{\phi }_j(t_n){\phi }_{j^{\prime }}(t_n).$$

La matrice de Gram est symétrique et positive. Elle n'est inversible que si les fonctions de base sont linéairement indépendantes une fois échantillonnées sur les $N$ entrées, c'est-à-dire si aucune colonne de $\mathbf{H}$ n'est combinaison des autres. Deux fonctions de base presque identiques sur l'échantillon rendent deux colonnes de $\mathbf{H}$ presque colinéaires, $\mathbf{G}$ presque singulière, et son inversion instable. C'est exactement ce qui s'observe ci-dessous, lorsque le nombre de fonctions de base devient trop grand.

Surapprentissage

Les fonctions de base rendent l'enrichissement du modèle trivial : une fonction de base de plus est une colonne de plus dans $\mathbf{H}$, une ligne de code avec np.vander. Prenons $N=18$ mesures d'un cosinus bruité, $x_n=\cos (2\pi \cdot 0,6t_n)+w_n$ avec $\sigma =0,3$ et des entrées $t_n$ régulièrement espacées dans $[-1,1]$, et ajustons des polynômes d'ordre $L$ croissant. Plus de fonctions de base, est-ce toujours mieux ?

Figure 2.10, Régression polynomiale d'ordre $L\in \{1,4,8,12\}$ sur $N=18$ mesures d'un cosinus bruité ($\sigma =0,3$), avec l'écart RMSE au signal vrai de chaque ajustement (script figures/ch2/03-surapprentissage-ordre.py).

Lorsque nous augmentons l'ordre, la covariance ${\sigma }^2({\mathbf{H}}^T\mathbf{H}{)}^{-1}$ du §2 augmente à mesure que les colonnes de $\mathbf{H}$ se multiplient et se ressemblent : le conditionnement de $\mathbf{H}$ passe de $1,7$ à $L=1$ à $4\cdot {10}^4$ à $L=12$. Les coefficients estimés explosent en conséquence : $\underset{l}{max}|{\hat{\theta }}_l|\approx 3000$ à $L=12$, contre $6$ à $L=4$. L'estimateur reste sans biais mais sa variance est devenue telle que chaque tirage de bruit produit un polynôme entièrement différent.

Confusion classique

« Plus de paramètres, donc un meilleur modèle. » Non : plus de paramètres améliore toujours l'ajustement aux données observées, jamais nécessairement l'écart au signal vrai. Le premier se mesure sur les points déjà acquis, le second sur ce que le modèle prétend prédire ; la figure 2.10 les montre évoluer en sens contraires dès $L>4$. Juger un modèle sur son erreur d'ajustement seule revient à noter un étudiant sur des annales qu'il a apprises par cœur.

4. La régularisation

Principe

Le polynôme de degré 12 du §3 pose un problème : $p=13$ paramètres pour $N=18$ points, des coefficients à $3000$ et une courbe qui oscille entre les mesures. Comment retenir les coefficients sans revenir à la main à un ordre plus faible ? La solution est l'estimateur du maximum a posteriori : munir $\mathit{\theta }$ d'un prior $p(\mathit{\theta })$ qui pénalise les coefficients invraisemblables, et minimiser $J_{\mathrm{MAP}}$.

Expression

Pour le modèle linéaire à bruit Gaussien blanc, le terme d'attache aux données de $J_{\mathrm{MAP}}$ est la loss des moindres carrés du §2, et l'estimateur prend une forme spécifique.

Régularisation

L'estimateur MAP minimise la loss des données augmentée d'un terme de régularisation :

$${\hat{\mathit{\theta }}}_{\mathrm{MAP}}=\arg \underset{\mathit{\theta }}{min}{J}_{\mathrm{MAP}}(\mathit{\theta }),$$

où

$$J_{\mathrm{MAP}}(\mathit{\theta })=\Vert \mathbf{x}-\mathbf{H}\mathit{\theta }{\Vert }^2+R(\mathit{\theta }),$$

et où $R(\mathit{\theta })=-\ln p(\mathit{\theta })$ est le terme de régularisation, imposé par le prior.

Exemple de prior

Le choix du prior est un choix de régularisation, exactement comme le choix du bruit était un choix de loss (§3 de la page précédente). Trois priors classiques, i.i.d. sur les composantes de $\mathit{\theta }$, donnent trois régularisations très utilisées :

Prior sur ${\theta }_i$	$J_{\mathrm{MAP}}(\mathit{\theta })$	Nom	Effet
Gaussien $\mathcal{N}(0,{\sigma }_0^2)$	$\Vert \mathbf{x}-\mathbf{H}\mathit{\theta }{\Vert }^2+\lambda \Vert \mathit{\theta }{\Vert }_2^2$	Ridge (L2)	rétrécit tous les coefficients, solution analytique
Laplace$(0,b)$	$\Vert \mathbf{x}-\mathbf{H}\mathit{\theta }{\Vert }^2+\lambda \Vert \mathit{\theta }{\Vert }_1$	LASSO (L1)	annule certains coefficients, parcimonie, pas de solution analytique
aucune loi propre (heuristique)	$\Vert \mathbf{x}-\mathbf{H}\mathit{\theta }{\Vert }^2+\lambda \Vert \mathit{\theta }{\Vert }_0$	Sélection de sous-ensemble (L0)	parcimonie exacte, mais combinatoire ($2^p$ sous-ensembles à tester)

Le cas Gaussien admet une solution fermée : c'est la régression ridge.

Proposition, Solution de la régression ridge

$${\hat{\mathit{\theta }}}_{\mathrm{ridge}}=({\mathbf{H}}^T\mathbf{H}+\lambda \mathbf{I}{)}^{-1}{\mathbf{H}}^T\mathbf{x},$$

où :

• $\lambda >0$ est le poids de régularisation, hérité du prior ($\lambda ={\sigma }^2/{\sigma }_0^2$).

Démonstration

En annulant le gradient de $J_{\mathrm{MAP}}(\mathit{\theta })=\Vert \mathbf{x}-\mathbf{H}\mathit{\theta }{\Vert }^2+\lambda \Vert \mathit{\theta }{\Vert }_2^2$ avec les règles de dérivation matricielle du §2 :

$$\mathrm{\nabla }{J}_{\mathrm{MAP}}(\mathit{\theta })=-2{\mathbf{H}}^T\mathbf{x}+2{\mathbf{H}}^T\mathbf{H}\mathit{\theta }+2\lambda \mathit{\theta }=\mathbf{0},$$

d'où la solution ${\hat{\mathit{\theta }}}_{\mathrm{ridge}}=({\mathbf{H}}^T\mathbf{H}+\lambda \mathbf{I}{)}^{-1}{\mathbf{H}}^T\mathbf{x}$.

La matrice ${\mathbf{H}}^T\mathbf{H}+\lambda \mathbf{I}$ est toujours inversible pour $\lambda >0$, même quand ${\mathbf{H}}^T\mathbf{H}$ est singulière ou mal conditionnée : la régularisation stabilise aussi le calcul numérique. ∎

Le prior Laplace produit le même rétrécissement, mais annule exactement certains coefficients plutôt que de les réduire : c'est la parcimonie, exploitée par le LASSO pour sélectionner des variables. Le prix en est la perte de la solution analytique : la minimisation devient numérique, comme au §6 de la page précédente. Le critère idéal, pénaliser $\Vert \mathit{\theta }{\Vert }_0$, le nombre de coefficients non nuls, est combinatoire : il faudrait essayer les $2^p$ sous-ensembles de coefficients possibles.

Exemple : régression polynomiale

Reprenons exactement le polynôme de degré 12 du §3 ($p=13$ paramètres, $N=18$ points du cosinus bruité) et faisons varier $\lambda$, en partant de $\lambda =0$, l'absence de régularisation. Si le prior Gaussien joue son rôle, un $\lambda$ bien choisi doit retenir les coefficients sans écraser le signal.

Figure 2.11, Régression ridge de degré 12 sur les $N=18$ mesures de la figure 2.10, pour $\lambda \in \{0,{10}^{-4},{10}^{-2},1\}$ (script figures/ch2/03-ridge-lambda.py).

À $\lambda =0$, les coefficients à $3000$ du §3 sont présents et la courbe oscille entre les mesures : c'est le panneau $L=12$ de la figure 2.10. Un $\lambda$ petit (${10}^{-4}$) suffit à ramener ces coefficients à une échelle raisonnable ($\approx 9$) et retrouve le cosinus sous-jacent. Un $\lambda$ trop grand ($1$) écrase le signal vers une courbe quasi plate : le biais domine, c'est le sous-apprentissage.

Le compromis se lit numériquement, pas seulement visuellement : l'écart au signal vrai (RMSE) vaut $0,43$ à $\lambda =0$, descend à $0,09$ à $\lambda ={10}^{-4}$, avant de remonter à $0,12$ ($\lambda ={10}^{-2}$) puis $0,41$ ($\lambda =1$), la même courbe en U que celle du compromis biais-variance de la MSE (page Fondations), ici visible directement sur l'allure de la courbe ajustée.