Splines et approximation au sens des moindres carrés

Objectifs du chapitre

À la fin de ce chapitre, vous serez capable de:

construire l'interpolant affine par morceaux d'un tableau de données, majorer son erreur par $\frac{h ^{2}}{8} max ∣ f^{''} ∣$ et dire pourquoi sa dérivée discontinue est parfois rédhibitoire;
compter les inconnues et les conditions d'une spline cubique — $4 n$ contre $4 n - 2$ — et en déduire qu'il reste exactement deux conditions à choisir, puis écrire les deux choix usuels: spline et spline ;

Ce que le chapitre 6 laisse ouvert

Le chapitre 6 a établi deux résultats qui, mis côte à côte, forment une impasse. D'un côté, par $n + 1$ points d'abscisses distinctes passe un et un seul polynôme de degré au plus $n$ : l'interpolation polynomiale est toujours possible, et elle est unique. De l'autre, la fonction de Runge $f (x) = \frac{1}{1 + 25 x ^{2}}$ sur $[- 1, 1]$ , interpolée sur des nœuds équidistants, donne une erreur maximale de pour , de pour et de pour : . Les nœuds de Tchebychev réparent le cas particulier, mais ils supposent qu'on choisit les abscisses — ce qui n'est jamais le cas devant un tableau de mesures.

Il faut donc renoncer à quelque chose. Le chapitre 6 renonçait au mauvais candidat: il gardait un seul polynôme et déplaçait les nœuds. Ce chapitre renonce au polynôme unique. L'idée est d'une simplicité désarmante — découper l'intervalle et changer de polynôme à chaque morceau — et elle ouvre sur la totalité de ce que les bibliothèques graphiques, les machines-outils et les logiciels de dessin utilisent aujourd'hui.

La seconde moitié du chapitre renonce à autre chose encore: à passer par les points. Quand les ordonnées sont des mesures entachées d'erreur, exiger que la courbe les traverse toutes revient à exiger qu'elle reproduise le bruit. On lui demandera alors de passer près, au sens le plus économique qui soit, et cette seconde histoire aboutira aux équations normales.

L'interpolation affine par morceaux

Commençons par le cas le plus pauvre et le plus utilisé. Donnons-nous une subdivision

a = x_{0} < x_{1} < \dots < x_{n} = b, h_{i} = x_{i + 1} - x_{i}, h = i max h_{i}, (7.1)

et les valeurs $y_{i} = f (x_{i})$ .

Son erreur se déduit en trois lignes de la formule d'erreur du chapitre 6. Sur un seul intervalle $[x_{i}, x_{i + 1}]$ , l'interpolation en deux nœuds donne, pour $f$ de classe $C^{2}$ ,

f (x) - p (x) = \frac{f ^{''} ( ξ _{x} )}{2} (x - x_{i}) (x - x_{i + 1}), ξ_{x} \in] x_{i}, x_{i + 1} [,

et le produit $(x - x_{i}) (x - x_{i + 1})$ atteint son maximum en module au milieu de l'intervalle, où il vaut $h_{i}^{2} /4$ . D'où

Démonstration. Fixons $x \in [x_{i}, x_{i + 1}]$ . Le polynôme $p$ restreint à cet intervalle est l'interpolant de $f$ aux deux nœuds $x_{i}$ et ; la formule d'erreur de l'interpolation à deux nœuds (chapitre 6) fournit un strictement intérieur tel que

f (x) - p (x) = \frac{f ^{''} ( ξ _{x} )}{2 !} ω (x), ω (x) = (x - x_{i}) (x - x_{i + 1}) .

La fonction $ω$ est un trinôme négatif sur l'intervalle, de sommet en $x = (x_{i} + x_{i + 1}) /2$ , où

∣ ω ∣ = (\frac{h _{i}}{2}) (\frac{h _{i}}{2}) = \frac{h _{i}^{2}}{4} .

Donc $∣ f (x) - p (x)∣ \leq \frac{1}{2} \cdot \frac{h _{i}^{2}}{4} max ∣ f^{''} ∣ = \frac{h _{i}^{2}}{8} max ∣ f^{''} ∣$ . En majorant par et en prenant le maximum sur , on obtient (7.3).

Voilà qui règle la divergence: sur la fonction de Runge, $max ∣ f^{''} ∣ = 50$ et la borne (7.3) donne $6, 25 h^{2}$ , qui tend vers zéro. L'interpolation par morceaux converge toujours, là où l'interpolation globale peut diverger. Le prix à payer est l'ordre: $h^{2}$ au lieu de la convergence spectaculaire d'un polynôme sur des nœuds bien choisis.

Le vrai défaut de l'interpolant affine: sa dérivée

Le défaut de $p$ n'est pas son erreur, c'est sa régularité. En chaque nœud intérieur, la pente saute de $\frac{y _{i} - y _{i - 1}}{h _{i - 1}}$ à : est continue mais pas dérivable, et sa dérivée seconde est nulle à l'intérieur des morceaux et indéfinie aux nœuds.

Cela peut être sans conséquence — pour lire une table de logarithmes, personne ne s'en soucie — et cela peut être fatal. Si la courbe décrit la trajectoire d'un outil de fraisage, la discontinuité de la pente est un à-coup; si elle décrit le profil d'une came, la discontinuité de la courbure est un choc d'accélération, donc une contrainte et un bruit; si elle sert de tracé routier ou ferroviaire, elle impose au véhicule une variation instantanée d'accélération latérale. La demande naturelle est donc: une courbe de classe $C^{2}$ , formée de morceaux polynomiaux du plus bas degré possible.

Posons la question en comptant. Sur $n$ intervalles, avec des morceaux de degré $m$ , il y a $n (m + 1)$ coefficients. Les conditions d'interpolation en imposent $2 n$ (chaque morceau doit passer par ses deux extrémités), et la continuité des dérivées d'ordre $1$ à $k$ aux $n - 1$ nœuds intérieurs en impose $k (n - 1)$ . Le décompte donne, pour les trois premiers degrés:

degré $m$	inconnues $n (m + 1)$	conditions pour $C^{m - 1}$	différence
1	$2 n$	$2 n$ (interpolation seule, $C^{0}$ )	$0$

Le degré 1 est rigide: aucune liberté, et pas de dérivée. Le degré 2 offre la classe $C^{1}$ mais une seule condition libre, qui doit être posée à une extrémité — l'information s'y propage ensuite de proche en proche, ce qui rend la courbe quadratique par morceaux asymétrique et sujette à des oscillations parasites. Le degré 3 offre la classe $C^{2}$ et deux conditions libres, une par extrémité: c'est le premier degré où le problème est symétrique. C'est la raison profonde pour laquelle la spline cubique est partout, et non un accident historique.

Question 7.1

Pourquoi le degré 3 est-il le bon compromis pour une courbe interpolante de classe $C^{2}$ ?

La spline cubique interpolante

Les quatre coefficients $(a_{i}, b_{i}, c_{i}, d_{i})$ de chaque morceau sont les inconnues. Le théorème qui suit est le décompte, et c'est la plus jolie arithmétique du chapitre: elle explique d'où vient la liberté résiduelle, donc pourquoi il existe plusieurs sortes de splines.

Démonstration. Les inconnues. Chaque morceau $s_{i}$ est un polynôme de degré au plus 3, donc porte $4$ coefficients dans l'écriture (7.4). Il y a $n$ morceaux, d'où $4 n$ inconnues. Le point 1 de la définition ne fournit aucune équation: il définit les inconnues.

Les conditions d'interpolation et la continuité de $s$ . Écrivons, pour chaque $i = 0, \dots, n - 1$ ,

s_{i} (x_{i}) = y_{i} et s_{i} (x_{i + 1}) = y_{i + 1} . (7.6)

Cela fait $2$ équations par morceau, soit $2 n$ équations. Elles contiennent à la fois le point 2 et la continuité de $s$ . En effet, $s$ est polynomiale donc continue à l'intérieur de chaque morceau; en un nœud intérieur $x_{i}$ ( $1 \leq i \leq n - 1$ ), (7.6) donne

s (x_{i}^{-}) = s_{i - 1} (x_{i}) = y_{i} = s_{i} (x_{i}) = s (x_{i}^{+}),

donc $s$ est continue en $x_{i}$ et y vaut $y_{i}$ . Aux deux bords, (7.6) donne $s (a) = y_{0}$ et $s (b) = y_{n}$ . Réciproquement, si est continue et interpole, alors chaque vaut en et en : (7.6) est exactement équivalente au point 2 joint à la continuité de . : chacune contraint un morceau différent ou une extrémité différente d'un même morceau.

La continuité de $s^{'}$ . À l'intérieur d'un morceau, $s^{'}$ existe et est continue. Le seul risque est aux $n - 1$ nœuds intérieurs $x_{1}, \dots, x_{n - 1}$ , où l'on impose

s_{i - 1}^{'} (x_{i}) = s_{i}^{'} (x_{i}), i = 1, \dots, n - 1, (7.7)

soit $n - 1$ équations. Rien n'est imposé en $x_{0}$ ni en $x_{n}$ , qui n'ont pas de voisin de l'autre côté.

La continuité de $s^{''}$ . Même raisonnement, aux mêmes nœuds:

s_{i - 1}^{''} (x_{i}) = s_{i}^{''} (x_{i}), i = 1, \dots, n - 1, (7.8)

soit encore $n - 1$ équations. Et c'est tout: $s^{'''}$ n'a pas à être continue — elle ne l'est d'ailleurs pas en général, car $s_{i}^{'''} = 6 d_{i}$ est constante par morceau et saute d'un morceau à l'autre. Une spline cubique est de classe et pas davantage, sauf coïncidence.

Le total. En additionnant,

interpolation 2 n + s^{'} continue (n - 1) + s^{''} continue (n - 1) = 4 n - 2, (7.9)

pour $4 n$ inconnues. L'ensemble des solutions, s'il est non vide, est donc un sous-espace affine de dimension au moins $4 n - (4 n - 2) = 2$ . Le théorème 7.3 montrera qu'en ajoutant deux équations bien choisies on obtient une solution et une seule: la dimension vaut donc exactement 2, et deux conditions supplémentaires ne sont pas seulement nécessaires, elles suffisent. $□$

Le système tridiagonal des moments

Résoudre directement les $4 n - 2$ équations en les $4 n$ coefficients serait maladroit. Le bon changement d'inconnues est classique et tient à une remarque: la dérivée seconde d'une spline cubique est affine par morceaux et continue, donc entièrement déterminée par ses valeurs aux nœuds.

Sur $[x_{i}, x_{i + 1}]$ , $s^{''}$ est affine et vaut $M_{i}$ en $x_{i}$ , en ; donc

s^{''} (x) = M_{i} \frac{x _{i + 1} - x}{h _{i}} + M_{i + 1} \frac{x - x _{i}}{h _{i}} . (7.11)

Deux intégrations, puis l'ajustement des deux constantes par $s (x_{i}) = y_{i}$ et $s (x_{i + 1}) = y_{i + 1}$ , donnent la forme qui sert dans tout le chapitre:

s (x) = \frac{M _{i} ( x _{i + 1} - x ) ^{3} + M _{i + 1} ( x - x _{i} ) ^{3}}{6 h _{i}} + (\frac{y _{i}}{h _{i}} - \frac{M _{i} h _{i}}{6}) (x_{i + 1} - x) + (\frac{y _{i + 1}}{h _{i}} - \frac{M _{i + 1} h _{i}}{6}) (x - x_{i}) . (7.12)

Par construction, (7.12) satisfait déjà l'interpolation (7.6) et la continuité de $s^{''}$ (7.8), puisque $s^{''}$ y est l'interpolant affine par morceaux des $M_{i}$ , qui est continu. Il ne reste donc à imposer que la continuité de $s^{'}$ .

En dérivant (7.12) et en écrivant $s_{i - 1}^{'} (x_{i}) = s_{i}^{'} (x_{i})$ , on obtient, après simplification, la relation à trois termes

h_{i - 1} M_{i - 1} + 2 (h_{i - 1} + h_{i}) M_{i} + h_{i} M_{i + 1} = 6 (\frac{y _{i + 1} - y _{i}}{h _{i}} - \frac{y _{i} - y _{i - 1}}{h _{i - 1}}), (7.13)

pour $i = 1, \dots, n - 1$ . Le membre de droite est six fois une différence divisée du second ordre: c'est la courbure discrète des données, exactement la quantité que la spline doit répartir.

En posant $α_{i} = 2 (h_{i - 1} + h_{i})$ et $δ_{i} = \frac{y _{i + 1} - y _{i}}{h _{i}} - \frac{y _{i} - y _{i - 1}}{h _{i - 1}}$ , le système s'écrit

α_{1} h_{1} h_{1} α_{2} ⋱ ⋱ ⋱ h_{n - 2} h_{n - 2} α_{n - 1} M_{1} M_{2} ⋮ M_{n - 1} = 6 δ_{1} δ_{2} ⋮ δ_{n - 1}, (7.14)

où l'on a supposé $M_{0}$ et $M_{n}$ connus — c'est précisément le rôle des deux conditions libres.

Spline naturelle: $M_{0} = M_{n} = 0$ , et (7.14) est le système complet, de taille $n - 1$ .
Spline encastrée: $s^{'} (a) = α$ et ajoutent les deux lignes et le système devient tridiagonal de taille , de même structure.

Démonstration. Dominance diagonale. Sur la ligne $i$ de (7.14), le terme diagonal vaut $α_{i} = 2 (h_{i - 1} + h_{i})$ et la somme des modules des termes hors diagonale vaut au plus $h_{i - 1} + h_{i}$ (elle vaut exactement cela sur les lignes intérieures, et moins sur la première et la dernière). Comme tous les sont strictement positifs,

∣ α_{i} ∣ = 2 (h_{i - 1} + h_{i}) > h_{i - 1} + h_{i} \geq j \neq = i \sum ∣ a_{ij} ∣ . (7.16)

Pour le système encastré, la première ligne donne $2 h_{0} > h_{0}$ et la dernière $2 h_{n - 1} > h_{n - 1}$ : la dominance stricte tient aussi.

Une matrice à diagonale strictement dominante est inversible. Supposons $A z = 0$ avec $z \neq = 0$ , et soit $k$ un indice où $∣ z_{k} ∣$ est maximal; alors $∣ z_{k} ∣ > 0$ . La ligne donne , d'où

∣ a_{k k} ∣ ∣ z_{k} ∣ \leq j \neq = k \sum ∣ a_{k j} ∣ ∣ z_{j} ∣ \leq ∣ z_{k} ∣ j \neq = k \sum ∣ a_{k j} ∣ < ∣ a_{k k} ∣ ∣ z_{k} ∣,

ce qui est absurde. Donc $z = 0$ et $A$ est inversible. (C'est le même argument qui, au chapitre 4, donne la convergence de la méthode de Jacobi sous dominance diagonale stricte.)

Conclusion. Le système (7.14) a donc une solution unique $(M_{1}, \dots, M_{n - 1})$ . La fonction $s$ définie morceau par morceau par (7.12) est alors bien une spline cubique interpolante: elle est cubique sur chaque morceau par construction; elle interpole, car (7.12) vaut $y_{i}$ en $x_{i}$ et en ; sa dérivée seconde est l'interpolant affine par morceaux des , donc continue, ce qui est (7.8); et sa dérivée première est continue aux nœuds intérieurs, ce qui est exactement (7.13). Elle est donc de classe . Réciproquement, toute spline naturelle a des moments qui vérifient (7.13) avec , donc coïncide avec celle-ci. Le raisonnement est identique pour la spline encastrée avec les deux lignes (7.15).

Un système tridiagonal se résout en O(n), pas en O(n cube)

La matrice de (7.14) est tridiagonale: chaque ligne a au plus trois coefficients non nuls. L'élimination de Gauss du chapitre 3 y devient l'algorithme de Thomas, qui n'est rien d'autre que l'élimination sans pivot exploitant la structure. Elle coûte $8 n - 7$ opérations flottantes, contre $\frac{2}{3} n^{3}$ pour une élimination complète: pour , contre , un rapport de plus de . Et le pivot partiel est inutile ici, précisément parce que la dominance diagonale stricte garantit que les pivots ne s'effondrent pas. L'argument tient en deux phrases: toute sous-matrice principale dominante d'une matrice à diagonale strictement dominante l'est encore, puisque supprimer des colonnes ne fait qu'enlever des termes positifs à la somme hors diagonale; elle est donc inversible par l'argument du théorème 7.3 ci-dessous, tous les mineurs principaux dominants sont non nuls, et la condition d'existence et d'unicité de la factorisation du théorème 3.1 (chapitre 3) est satisfaite.

def resoudre_tridiagonal(a, b, c, d):
    """Algorithme de Thomas. a: sous-diagonale, b: diagonale, c: sur-diagonale, d: second membre."""
    n = len(b)
    cp = [0.0] * n
    dp = [0.0] * n
    cp[0] = c[0] / b[0]
    dp[0] = d[0] / b[0]

Comptons exactement. La première ligne coûte deux divisions. Chacune des $n - 2$ lignes intérieures en coûte six: une multiplication et une soustraction pour le pivot m, une division pour cp, puis une multiplication, une soustraction et une division pour dp. La dernière ligne n'en coûte que cinq, puisque le garde if i < n - 1 lui épargne la division de cp. La descente vaut donc $2 + 6 (n - 2) + 5 = 6 n - 5$ opérations. La remontée traite $n - 1$ lignes à une multiplication et une soustraction chacune, soit . Au total opérations, linéaire en , ce qui est exactement le compte annoncé au chapitre 3 pour l'algorithme de Thomas. Les listes et ont la longueur de la diagonale par commodité, avec et ignorés.

Question 7.2

Écrivez l'algorithme de Thomas. La fonction resoudre_tridiagonal(a, b, c, d) reçoit la sous-diagonale, la diagonale, la sur-diagonale et le second membre, et renvoie la solution. Le programme la teste sur un système $3 \times 3$ dont la solution est $(1, 1, 1)$ .

thomas.py

Question 7.3

Construisez les moments de la spline cubique naturelle. La fonction moments_naturels(xs, ys) doit assembler le système (7.13), le résoudre avec resoudre_tridiagonal et renvoyer la liste des $n + 1$ moments, extrémités nulles comprises. Le programme l'applique aux données de l'exemple 7.1.

spline.py

L'erreur d'une spline cubique

Ce résultat est admis: sa démonstration complète repose sur une majoration uniforme de l'inverse de la matrice tridiagonale (7.14) en norme infinie, qui relève de la théorie de la perturbation matricielle plutôt que de ce chapitre; on la trouve dans Quarteroni, Sacco et Saleri, chapitre 8, et dans de Boor, A Practical Guide to Splines. Nous nous contenterons ici de la vérifier numériquement, ce qui est d'ailleurs la seule vérification qui protège d'une erreur d'implémentation.

Prenons $f (x) = e^{x}$ sur $[0, 1]$ , pour laquelle $f^{(4)} = e^{x}$ et $max ∣ f^{(4)} ∣ = e$ . La table suivante donne l'erreur maximale mesurée sur une grille de 4001 points, pour la spline encastrée (exacte aux bords) et pour la spline naturelle (qui impose alors que vaut et ).

$n$	$h$	spline encastrée	rapport	spline naturelle	rapport
2	0,5	$4, 008 \cdot 1 0^{- 4}$	—	$2, 978 \cdot 1 0^{- 2}$	—
4	0,25	$2, 636 \cdot 1 0^{- 5}$

La lecture est sans ambiguïté. La colonne des rapports de la spline encastrée tend vers 16: diviser $h$ par deux divise l'erreur par $2^{4}$ , c'est l'ordre 4. Celle de la spline naturelle tend vers 4: ordre 2 seulement. La différence tient uniquement aux deux conditions de bord, et elle est brutale — à $n = 32$ , la spline encastrée est $19000$ fois plus précise que la naturelle sur les mêmes nœuds.

Vérifions aussi la borne (7.18) elle-même: à $h = 0, 5$ elle vaut $\frac{5}{384} \cdot 0, 0625 \cdot e = 2, 212 \cdot 1 0^{- 3}$ , à comparer à l'erreur mesurée ; à , elle vaut contre mesurés. La borne est respectée dans les deux cas, avec un facteur de sécurité d'environ 5 — ce qui est le comportement normal d'une majoration uniforme.

Question 7.4

Remettez dans l'ordre les étapes du calcul d'une spline cubique naturelle à partir d'un tableau de données.

Glissez les éléments pour les mettre dans le bon ordre

Assembler le système tridiagonal (7.14) de taille

n - 1

, de diagonale

2 (h_{i - 1} + h_{i})

Le résoudre par l'algorithme de Thomas, en

8 n - 7

opérations

Déduire les coefficients

a_{i}, b_{i}, c_{i}, d_{i}

de chaque morceau par les formules (7.17)

Calculer les pas

h_{i} = x_{i + 1} - x_{i}

et les différences divisées secondes

δ_{i}

Évaluer la spline en un point

x

après avoir localisé l'intervalle qui le contient

Encadrer la solution des deux moments nuls

M_{0} = M_{n} = 0

B-splines, Bézier, et ce que fait vraiment un logiciel de dessin

La forme (7.4) est commode pour calculer mais mauvaise pour manipuler: modifier une donnée $y_{k}$ change tous les moments, donc tous les morceaux, puisque le système (7.14) est global. Un dessinateur veut l'inverse: bouger un point de contrôle et ne modifier la courbe que localement.

La réponse est un changement de base. On peut montrer que l'espace des splines cubiques sur une subdivision donnée, de dimension $n + 3$ , admet une base de fonctions $B_{i, 3}$ — les B-splines — qui sont elles-mêmes des splines cubiques, positives, de somme égale à 1 en tout point, et à support borné: chaque $B_{i, 3}$ est nulle en dehors de quatre intervalles consécutifs. Écrire $s = \sum_{i} c_{i} B_{i, 3}$ rend la représentation locale: changer ne modifie que quatre morceaux. La matrice du problème d'interpolation devient bande, et le coût reste linéaire.

Une autre famille mérite d'être nommée, parce que c'est celle que les bibliothèques de tracé emploient le plus souvent: la spline cubique d'Hermite par morceaux. On y impose, en chaque nœud, la valeur et la pente; chaque morceau est alors déterminé par ses quatre conditions locales, sans aucun système à résoudre, et la courbe est de classe $C^{1}$ mais pas $C^{2}$ . On échange donc une dérivée seconde continue contre la localité et un coût nul — et, si l'on choisit les pentes de façon à préserver la monotonie des données (schéma de Fritsch–Carlson, pchip dans Matlab et SciPy), on obtient en prime une courbe qui n'invente pas d'oscillation entre deux points voisins, ce qu'une spline cubique ordinaire, elle, peut faire.

En pratique, les logiciels de dessin vectoriel et les formats de police n'interpolent d'ailleurs pas: ils approchent. Le tracé d'un chemin en PostScript, en PDF ou en SVG est une suite de courbes de Bézier cubiques, définies par quatre points de contrôle dont deux seulement sont sur la courbe; les deux autres sont les poignées que l'on tire à la souris. Les polices TrueType utilisent des Bézier quadratiques, les polices PostScript et OpenType/CFF des cubiques. La conception assistée par ordinateur et l'usinage emploient les NURBS, qui sont des B-splines rationnelles à poids, parce qu'elles représentent exactement les coniques — un cercle n'est pas une courbe polynomiale, mais c'est une courbe rationnelle. Toute cette famille descend du même décompte que celui du théorème 7.2; ce qui change est la base choisie et le fait d'imposer, ou non, le passage par les points.

Question 7.5

On interpole 21 points par une spline cubique. Combien de coefficients le problème comporte-t-il au total, en comptant quatre coefficients par morceau?

Passer près plutôt que passer par

Changeons de question. Un capteur renvoie six couples $(m_{i}, U_{i})$ censés être alignés, et ils ne le sont pas, parce qu'aucune mesure n'est exacte. Faire passer une courbe par les six points est absurde deux fois: la courbe reproduirait le bruit, et sa valeur en un point non mesuré serait plus fausse que la droite qu'on cherchait.

Le problème et sa formulation matricielle

On se donne $n$ couples $(x_{i}, y_{i})$ , $i = 1, \dots, n$ , et une famille de $m + 1$ fonctions de base avec — typiquement , mais rien ne l'impose. On cherche les coefficients qui rendent la fonction la plus proche possible des données.

Trois remarques avant la théorie. Le problème est linéaire en $c$ , même si $g$ ne l'est pas en $x$ : ajuster une parabole est un problème linéaire. Les résidus sont verticaux, pas orthogonaux à la courbe: on suppose que $x$ est connu exactement et que toute l'erreur est sur $y$ ; si les deux variables sont bruitées, il faut une autre méthode (la régression orthogonale, ou moindres carrés totaux). Et on minimise des carrés et non des valeurs absolues, pour deux raisons: la fonction $S$ est alors différentiable partout, ce qui donne un système linéaire; et la solution est l'estimateur de vraisemblance maximale si les erreurs sont gaussiennes indépendantes de même variance. Le prix de ce choix est la sensibilité aux points aberrants, et nous le paierons plus loin.

Les équations normales

Démonstration. Posons $r^{*} = y - A c^{*}$ et écrivons un vecteur quelconque sous la forme $c = c^{*} + v$ . Alors $A c - y = A v - r^{*}$ , et en développant le carré de la norme euclidienne,

S (c^{*} + v) = ∥ A v - r^{*} ∥_{2}^{2} = ∥ r^{*} ∥_{2}^{2} - 2 v^{T} A^{T} r^{*} + ∥ A v ∥_{2}^{2} . (7.22)

C'est une identité exacte, valable pour tout $v$ : la fonction $S$ est un polynôme du second degré en $v$ , sans aucun reste.

Condition suffisante. Supposons $A^{T} r^{*} = 0$ , c'est-à-dire $A^{T} A c^{*} = A^{T} y$ . Le terme croisé de (7.22) disparaît et il reste

S (c^{*} + v) = S (c^{*}) + ∥ A v ∥_{2}^{2} \geq S (c^{*})

pour tout $v$ : $c^{*}$ est bien un minimum global.

Condition nécessaire. Supposons réciproquement que $c^{*}$ minimise $S$ . Fixons un $v$ quelconque et considérons la fonction d'une variable réelle

φ (t) = S (c^{*} + t v) = S (c^{*}) - 2 t v^{T} A^{T} r^{*} + t^{2} ∥ A v ∥_{2}^{2},

obtenue en remplaçant $v$ par $t v$ dans (7.22). C'est un trinôme en $t$ , dérivable, qui atteint son minimum en $t = 0$ par hypothèse. Donc $φ^{'} (0) = - 2 v^{T} A^{T} r^{*} = 0$ . Ceci valant pour , en particulier pour les vecteurs de la base canonique, on conclut , soit (7.21).

Lecture géométrique. La $j$ -ième composante de $A^{T} r^{*}$ est le produit scalaire de la $j$ -ième colonne de $A$ avec $r^{*}$ . Dire que $A^{T} r^{*} = 0$ , c'est donc dire que le résidu est orthogonal à chaque colonne, donc à tout leur espace engendré . Autrement dit, est la de sur , et le théorème de la projection garantit qu'elle existe et est unique.

Unicité de $c^{*}$ . Si $A$ est de rang plein en colonnes et si $A^{T} A z = 0$ , alors $z^{T} A^{T} A z = ∥ A z ∥_{2}^{2} = 0$ , donc , donc par injectivité. La matrice , qui est manifestement symétrique, est donc inversible; et comme pour , elle est définie positive. Le minimum est atteint en l'unique .

Pour la droite $g (x) = c_{0} + c_{1} x$ , le système (7.21) est un $2 \times 2$ que l'on peut écrire explicitement:

(n \sum x_{i} \sum x_{i} \sum x_{i}^{2}) (c_{0} c_{1}) = (\sum y_{i} \sum x_{i} y_{i}), (7.23)

de déterminant $D = n \sum x_{i}^{2} - (\sum x_{i})^{2}$ , strictement positif dès que deux abscisses diffèrent (c'est, à un facteur $n$ près, la variance des ).

Exemple 7.3 · Un étalonnage, calculé en fractions exactes

Un capteur de force est étalonné en suspendant des masses connues et en lisant la tension de sortie. Les six mesures — fictives, avec une incertitude de lecture volontairement grossière, de l'ordre de $0, 05$ V, pour que les résidus soient visibles — sont:

$m_{i}$ (kg)	0	2	4	6	8	10
$U_{i}$ (V)	0,06	0,14	0,48	0,55	0,90

Le coefficient de détermination

Comment dire si l'ajustement est bon? La somme $S$ dépend des unités et du nombre de points. On la compare donc à la dispersion des données elles-mêmes.

La décomposition est un théorème de Pythagore: $y - \overset{y}{ˉ} 1 = (y - A c^{*}) + (A c^{*} - \overset{y}{ˉ} 1)$ , et les deux termes de droite sont orthogonaux puisque le premier est le résidu et que le second appartient à — à condition que y appartienne aussi, ce qui exige la présence de la constante dans la base.

Sur l'exemple 7.3: $\overset{ˉ}{U} = \frac{13}{25} = 0, 52$ V, $SCT = \frac{3619}{5000} = 0, 7238$ , et donc exactement. On vérifie , et

R^{2} = 1 - \frac{0 , 0238}{0 , 7238} = \frac{500}{517} = 0, 9671.

Question 7.6

Écrivez droite_moindres_carres(xs, ys) qui résout le système (7.23) et renvoie le couple $(c_{0}, c_{1})$ , ainsi que coefficient_determination qui renvoie $R^{2}$ . Le programme les applique aux six mesures de l'exemple 7.3.

moindres_carres.py

Explorateur 7.1 · Un nuage de points, une droite, et un point qui s'écarte

Réglez d'abord la dispersion du nuage: la droite en trait plein suit, mais elle reste centrée sur le modèle. Déplacez ensuite le dernier point, celui qui est entouré, vers le haut du cadre: la droite en trait plein bascule et la somme des carrés explose, alors que la droite bleue — la même régression calculée sans ce point — ne bouge pas d'un pixel. C'est toute la fragilité des moindres carrés devant une donnée aberrante. Augmentez enfin le nombre de points: le même écart pèse moins lourd, sans jamais devenir inoffensif.

Nombre de points n8 points

Dispersion du nuage σ (V)0,08

Dernier point mesuré, en m = 10 kg (V)1,00

Pente c₁ (tous les points)

0,1083

Ordonnée c₀

−0,0201

Pente c₁ sans le point entouré

0,1145

Somme des carrés S

0,0323

Coefficient R²

0,9689

Résidu du point entouré

−0,063

Régression polynomiale et le mauvais conditionnement de A transposée A

Rien n'oblige à ajuster une droite. Avec $φ_{j} (x) = x^{j}$ , $j = 0, \dots, m$ , on ajuste un polynôme de degré $m$ à $n > m + 1$ points, et les équations normales s'écrivent, coefficient par coefficient,

(A^{T} A)_{j k} = i = 1 \sum n x_{i}^{j + k}, (A^{T} y)_{j} = i = 1 \sum n x_{i}^{j} y_{i} . (7.25)

La matrice $A^{T} A$ est la matrice de Gram de la famille $(φ_{0}, \dots, φ_{m})$ pour le produit scalaire discret $⟨ u, v ⟩ = \sum_{i} u (x_{i}) v (x_{i})$ ; pour la base monomiale, c'est une matrice de de la distribution des abscisses. Elle est symétrique définie positive (théorème 7.5), donc justiciable de la factorisation de Cholesky du chapitre 3, et il est tentant de s'arrêter là: on assemble, on résout par Gauss (chapitre 3), et c'est fini. C'est précisément ce qu'il ne faut pas faire au-delà du degré 3 ou 4, et la raison est le conditionnement du chapitre 4.

Le calcul qui fait peur

Prenons $n = 11$ points d'abscisses $x_{i} = i /10$ , $i = 0, \dots, 10$ , régulièrement réparties sur $[0, 1]$ , et faisons croître le degré. Le nombre de conditionnement $κ_{\infty} (A^{T} A) = ∥ A^{T} A ∥_{\infty} ∥(A^{T} A)^{- 1} ∥_{\infty}$ , calculé pour qu'il ne soit pas lui-même victime de l'arrondi, vaut:

degré $m$	taille de $A^{T} A$	$κ_{\infty} (A^{T} A)$	chiffres perdus
1	$2 \times 2$	22,5	1,4
2	$3 \times 3$

Chaque degré supplémentaire multiplie $κ$ par un facteur voisin de 32, soit une perte d'environ un chiffre et demi par degré. La dernière colonne est $lo g_{10} κ$ : c'est, d'après le chapitre 4, le nombre de chiffres significatifs qu'une perturbation relative des données peut faire disparaître du résultat. Au degré 6, sur seize chiffres disponibles en double précision, il en reste sept.

D'où cela vient-il? De ce que les colonnes $1, x, x^{2}, \dots$ deviennent presque colinéaires sur un intervalle borné: sur $[0, 1]$ , les fonctions $x^{5}$ et $x^{6}$ ont un cosinus d'angle supérieur à $0, 99$ . La matrice (7.25) est une somme de Riemann: pour des abscisses équiréparties sur ,

(A^{T} A)_{j k} = i \sum x_{i}^{j + k} \approx n \int_{0}^{1} x^{j + k} d x = \frac{n}{j + k + 1}, (7.26)

c'est-à-dire $n$ fois la matrice de Hilbert $H_{m + 1}$ — l'exemple canonique de matrice mal conditionnée, celui du chapitre 4 et de l'annexe A. Le conditionnement étant invariant par multiplication par un scalaire, $κ (A^{T} A) \approx κ (H_{m + 1})$ , et l'on connaît ces valeurs exactement: , , . Comparez ligne à ligne avec la table: contre , contre . , et c'est pour cela que le chapitre 4 était nécessaire avant celui-ci.

Les trois remèdes

Centrer et réduire. Le conditionnement dépend de l'intervalle. Sur les mêmes 11 abscisses translatées sur $[1, 2]$ , au degré 5, $κ_{\infty} (A^{T} A)$ passe de $1, 88 \cdot 1 0^{7}$ à $6, 24 \cdot 1 0^{11}$ : quatre ordres de grandeur perdus pour une simple translation, parce que sur toutes les puissances valent à peu près la même chose. À l'inverse, en posant et en ajustant en , on obtient au degré 5 — que sur et deux cent cinquante millions de fois mieux que sur , pour un changement de variable qui coûte deux opérations par point. C'est le remède le moins cher et le plus souvent oublié.

Changer de base. Rien n'oblige à prendre $φ_{j} (x) = x^{j}$ . Si l'on choisit des polynômes $φ_{j}$ orthogonaux pour le produit scalaire discret $⟨ u, v ⟩ = \sum_{i} u (x_{i}) v (x_{i})$ — obtenus par exemple par Gram–Schmidt sur la base monomiale, ou par la relation de récurrence à trois termes des polynômes de Legendre discrets —, alors devient : les équations normales se découplent et chaque coefficient se lit directement,

c_{j} = \frac{\sum _{i} φ _{j} ( x _{i} ) y _{i}}{\sum _{i} φ _{j} ( x _{i} ) ^{2}} . (7.27)

Après normalisation des colonnes, $A^{T} A = I$ et $κ = 1$ : le conditionnement disparaît par construction. Avantage supplémentaire, ajouter un degré ne change aucun des coefficients déjà calculés, ce qui permet d'augmenter le degré progressivement et de s'arrêter quand $SCR$ cesse de diminuer significativement.

Factoriser. La méthode de référence est la factorisation $A = QR$ avec $Q$ à colonnes orthonormées et $R$ triangulaire supérieure: les équations normales deviennent $R c = Q^{T} y$ , une simple remontée triangulaire, et l'on n'a jamais formé $A^{T} A$ . C'est ce que font numpy.linalg.lstsq, Matlab et LAPACK. La construction de et relève des transformations de Householder, hors du programme de ce chapitre; retenez le principe et le fait que le coût, opérations environ, reste du même ordre que celui des équations normales.

Question 7.7

Écrivez equations_normales(xs, ys, m) qui construit la matrice $A^{T} A$ et le vecteur $A^{T} y$ de la régression polynomiale de degré $m$ , **sans construire $A$ ** — les formules (7.25) suffisent. Le solveur de Gauss est fourni. Le programme affiche la somme des carrés pour les degrés 1, 2 et 3 sur les données de l'exemple 7.3.

regression.py

def resoudre_gauss(A, b):
  """Elimination de Gauss avec pivot partiel."""
  n = len(b)
  M = [A[i][:] + [b[i]] for i in range(n)]
  for k in range(n):
      p = max(range(k, n), key=lambda r: abs(M[r][k]))
      M[k], M[p] = M[p], M[k]
      for r in range(k + 1, n):
          f = M[r][k] / M[k][k]
          for j in range(k, n + 1):
              M[r][j] -= f * M[k][j]
  x = [0.0] * n
  for i in range(n - 1, -1, -1):
      s = M[i][n] - sum(M[i][j] * x[j] for j in range(i + 1, n))
      x[i] = s / M[i][i]
  return x

def equations_normales(xs, ys, m):
  """Doit renvoyer (A^T A, A^T y) pour la base 1, x, ..., x^m."""
  n = len(xs)
  G = [[float(n), sum(xs)], [sum(xs), sum(x * x for x in xs)]]
  v = [sum(ys), sum(xs[i] * ys[i] for i in range(n))]
  return G, v

def regression_polynomiale(xs, ys, m):
  G, v = equations_normales(xs, ys, m)
  return resoudre_gauss(G, v)

def somme_carres(xs, ys, c):
  s = 0.0
  for i in range(len(xs)):
      p = sum(c[j] * xs[i] ** j for j in range(len(c)))
      s += (ys[i] - p) ** 2
  return s

masses = [0.0, 2.0, 4.0, 6.0, 8.0, 10.0]
tensions = [0.06, 0.14, 0.48, 0.55, 0.90, 0.99]
for m in [1, 2, 3]:
  c = regression_polynomiale(masses, tensions, m)
  print("degre %d : S = %.6f" % (m, somme_carres(masses, tensions, c)))
print("c (degre 2) = " + " ".join("%.6f" % v for v in regression_polynomiale(masses, tensions, 2)))

Question 7.8

Vous ajustez un polynôme de degré 6 à 11 points d'abscisses régulièrement réparties, en formant $A^{T} A$ et en résolvant par élimination de Gauss en double précision. À quoi devez-vous vous attendre?

Le point aberrant

Les moindres carrés minimisent une somme de carrés: un résidu deux fois plus grand pèse quatre fois plus. Un point faux pèse donc démesurément, et l'effet est facile à chiffrer sur les données de l'exemple 7.3.

Supposons qu'en relevant la dernière mesure l'opérateur ait heurté la table et lu $1, 80$ V au lieu de $0, 99$ V. Un seul nombre change, sur six. Les nouvelles sommes sont $\sum U_{i} = 3, 93$ et $\sum m_{i} U_{i} = 30, 70$ , d'où

c_{0} = - \frac{47}{350} = - 0, 1343 V, c_{1} = \frac{221}{1400} = 0, 1579 V/kg .

Comparons aux valeurs saines $c_{0} = 0, 02$ et $c_{1} = 0, 10$ :

grandeur	données saines	avec le point faux	sans le point du tout
pente $c_{1}$	0,1000	0,1579	0,1045
ordonnée $c_{0}$	0,0200	−0,1343	0,0080
somme $S$	0,0238	0,2876	0,0219
$R^{2}$	0,9671	0,8584	0,9522

La pente augmente de 57,9 %, et l'ordonnée à l'origine change de signe — un capteur qui délivrerait $- 134$ mV à vide. Mais la ligne décisive est la troisième colonne: supprimer purement et simplement la mesure ne déplace la pente que de $0, 1000$ à $0, 1045$ , soit $4, 5 %$ . Autrement dit, un point faux fait treize fois plus de dégâts qu'un point manquant. C'est le résultat à retenir: devant une donnée suspecte, l'écarter est presque toujours moins coûteux que la garder.

Une mise en garde sur la dernière ligne du tableau, qui est un piège tendu par $R^{2}$ lui-même. La troisième colonne affiche $0, 9522$ , moins que les $0, 9671$ de la colonne saine, et l'on serait tenté d'en conclure que retirer le point a dégradé l'ajustement. C'est faux, parce que les trois $R^{2}$ ne portent pas sur les mêmes données: en passant de six à cinq points on retire la mesure la plus éloignée de la moyenne, donc $SCT$ chute de $0, 7238$ à $0, 4587$ alors que bouge à peine, de à . Le quotient monte mécaniquement. ; d'un jeu à l'autre, seuls les coefficients et l'écart-type résiduel sont comparables — ici vaut V sur six points et V sur cinq, deux valeurs du même ordre, ce qui est la bonne lecture.

Deux mécanismes se combinent. Le premier est la pondération quadratique. Le second est le bras de levier: la mesure fautive est à l'extrémité de l'intervalle des abscisses, là où un écart fait pivoter la droite au maximum. La même erreur commise au centre du nuage aurait déplacé l'ordonnée sans presque toucher la pente. L'explorateur ci-dessus permet de vérifier les deux effets, en déplaçant le point marqué et en augmentant $n$ : la droite ajustée sans lui, en bleu, ne bouge jamais, et c'est la mesure exacte du dégât.

Que faire d'un point aberrant

Dans l'ordre. Un: tracer les résidus et regarder. Un résidu supérieur à trois ou quatre écarts-types résiduels mérite une enquête, jamais une suppression automatique. Deux: revenir au cahier de manipulation. Un point aberrant a presque toujours une cause identifiable — une unité, une virgule, un capteur saturé, une saisie. Trois: si le point est réel et qu'on veut néanmoins un ajustement robuste, changer de critère. Minimiser $\sum ∣ r_{i} ∣$ (moindres écarts absolus, ou régression $L^{1}$ ) donne un problème de programmation linéaire, sans formule fermée mais très peu sensible aux aberrants; minimiser la médiane des carrés est plus robuste encore. Quatre: ne jamais supprimer un point au seul motif qu'il dégrade $R^{2}$ . La différence entre écarter une mesure fausse et écarter une mesure gênante est toute la différence entre la science et son contraire.

Question 7.9

Reprenez les six mesures de l'exemple 7.3 en remplaçant la dernière par $U_{6} = 1, 80$ V. Calculez la nouvelle pente $c_{1}$ , en V/kg.

Modèles linéarisables, et le piège du logarithme

Beaucoup de lois ne sont pas linéaires en leurs paramètres, et pourtant un changement de variable les ramène au cas précédent:

y = a e^{b x} ⟹ ln y = ln a + b x, y = a x^{b} ⟹ ln y = ln a + b ln x, y = \frac{1}{a + b x} ⟹ \frac{1}{y} = a + b x . (7.28)

La recette est universelle: on transforme, on ajuste une droite par (7.23), on revient. Elle est correcte du point de vue algébrique, et elle ne résout pas le problème posé. Car ajuster $ln y$ minimise $\sum_{i} (ln y_{i} - ln a - b x_{i})^{2}$ , alors que le problème d'origine demandait de minimiser . Ce sont deux critères différents, et ils n'ont pas le même minimum.

Le logarithme comprime les grandes valeurs et dilate les petites: minimiser l'écart sur $ln y$ revient à peu près à minimiser l'écart relatif sur $y$ , donc à accorder aux mesures les plus petites — souvent les plus bruitées, celles de la fin d'une décroissance — un poids qu'elles ne méritent pas.

Linéariser change le problème, et il faut le dire

Ajuster le logarithme est légitime, rapide et souvent suffisant — et c'est une excellente façon de trouver un point de départ. Mais trois précautions s'imposent. Un: dire dans le rapport ce qui a été minimisé. «Ajustement par moindres carrés sur $ln U$ » et «ajustement par moindres carrés sur $U$ » ne sont pas la même phrase et ne donnent pas les mêmes coefficients. Deux: si le bruit est additif et d'amplitude constante (un instrument à résolution fixe), le critère correct porte sur $y$ ; s'il est multiplicatif, d'amplitude proportionnelle au signal (une mesure à erreur relative constante), alors le critère logarithmique est le bon, et c'est la linéarisation qui est justifiée. La question n'est donc pas «quelle méthode est meilleure» mais «quel est le modèle d'erreur». Trois: une réparation intermédiaire existe, la pondération. Minimiser $\sum_{i} w_{i} (ln y_{i} - ln a - b x_{i})^{2}$ avec compense au premier ordre la distorsion du logarithme, puisque ; on retrouve alors, à l'ordre un, l'ajustement sur sans sortir du linéaire.

Synthèse

L'interpolation par morceaux échange l'ordre contre la robustesse. L'interpolant affine a une erreur $\frac{h ^{2}}{8} max ∣ f^{''} ∣$ qui tend toujours vers zéro, là où l'interpolation globale peut diverger; son défaut est sa dérivée discontinue, inacceptable dès qu'une trajectoire physique est en jeu.
Le décompte est la clé de la spline cubique. Pour $n$ intervalles: $4 n$ coefficients, $2 n$ conditions d'interpolation, de continuité de et de continuité de , soit conditions. Il reste conditions libres, ce qui engendre les familles naturelle ( nulle aux bords), encastrée ( imposée), périodique et .

Exercices

Vous pouvez afficher le corrigé directement sous chaque énoncé après avoir cherché la solution.

Exercice 7.1 · Décompte et conditions de bord

Une spline cubique interpole 15 points. Donner le nombre d'inconnues, le nombre de conditions imposées par la définition, et le nombre de conditions libres.
Écrire les deux équations supplémentaires d'une spline naturelle, puis celles d'une spline encastrée avec $s^{'} (a) = 1$ et $s^{'} (b) = - 2$ .
La condition not-a-knot impose la continuité de $s^{'''}$ en et en . Montrer qu'elle revient à dire que et sont le même polynôme.

Exercice 7.2 · Une spline naturelle sur trois points

On interpole $(0, 0)$ , $(1, 1)$ , $(2, 0)$ par une spline cubique naturelle.

Écrire et résoudre le système des moments.
Donner les deux morceaux $s_{0}$ et $s_{1}$ explicitement.

Exercice 7.3 · Une droite des moindres carrés, à la main

On dispose des cinq mesures suivantes, où $x$ est un temps en minutes et $y$ une concentration en mg/L:

$x_{i}$	1	2	3	4	5
$y_{i}$	3,0	4,6	5,8	7,4	8,2

Calculer les quatre sommes nécessaires et résoudre les équations normales.

Exercice 7.4 · Base orthogonale et conditionnement

On ajuste un polynôme de degré 2 aux cinq abscisses $x_{i} = - 2, - 1, 0, 1, 2$ .

Écrire $A^{T} A$ pour la base monomiale $1, x, x^{2}$ et calculer .

Exercice 7.5 · Le résidu des moindres carrés est orthogonal, et ce qu'on en déduit

Soit $A \in R^{n \times p}$ de rang plein en colonnes, $y \in R^{n}$ , et $c^{*}$ la solution des équations normales. On pose .

Références

Quarteroni, A., Sacco, R. et Saleri, F., Méthodes numériques — algorithmes, analyse et applications, Springer, Milan, chap. 8 (interpolation par morceaux, splines et approximation au sens des moindres carrés).
de Boor, C., A Practical Guide to Splines, éd. révisée, Springer, New York (la référence sur les B-splines et la condition not-a-knot).
Rappaz, J. et Picasso, M., Introduction à l'analyse numérique, Presses polytechniques et universitaires romandes, Lausanne, chap. 5.
Burden, R. L. et Faires, J. D., Numerical Analysis, Cengage, Boston, chap. 3 (splines cubiques) et chap. 8 (approximation discrète).
Trefethen, L. N. et Bau, D., Numerical Linear Algebra, SIAM, Philadelphie, leçons 11 et 18 à 20 (moindres carrés, factorisation $QR$ et conditionnement).
Björck, Å., Numerical Methods for Least Squares Problems, SIAM, Philadelphie (l'ouvrage de référence sur les équations normales et leurs alternatives).

$n$	polynôme unique de degré $n$	spline cubique naturelle
5	0,433	0,423
10	1,92	0,0220
15	2,11	0,0309
20	—	0,00318
40	—	0,000278

critère	ajustement linéarisé	vrais moindres carrés
$\sum (U_{i} - a e^{b t_{i}})^{2}$ , en V $^{2}$	0,033088	0,022451
$\sum (ln U_{i} - ln a - b t_{i})^{2}$	0,049864	0,051782

Objectifs du chapitre

Ce que le chapitre 6 laisse ouvert

L'interpolation affine par morceaux

Le vrai défaut de l'interpolant affine: sa dérivée

La spline cubique interpolante

Le système tridiagonal des moments

L'erreur d'une spline cubique

B-splines, Bézier, et ce que fait vraiment un logiciel de dessin

Passer près plutôt que passer par

Le problème et sa formulation matricielle

Les équations normales

Le coefficient de détermination

Régression polynomiale et le mauvais conditionnement de A transposée A

Le calcul qui fait peur

Les trois remèdes

Le point aberrant

Modèles linéarisables, et le piège du logarithme

Synthèse

Exercices

Références

Compter les conditions

Assembler le second membre

Résoudre le système 2 par 2

Le premier morceau, et sa vérification