Systèmes non déterminés

On rencontre souvent des problèmes de moindres-carrés. Mathématiquement, ces problèmes consistent à minimiser le carré de la norme euclidienne d’une fonction à valeurs vectorielles, qui peut être linéaire (on parle alors de moindres-carrés linéaire, voir la section 19.1) ou non linéaire (on parle alors de moindres-carrés non linéaire, voir la section 19.3). Ce chapitre leur est consacré. Étant donné la grande variété de problèmes qui peuvent être modélisés comme problème de moindres-carrés, on les rencontre sous des appellations différentes : estimation de paramètres, problème d’identification, de calibration ou de régression en statistiques [199, 333]. 19.1 Moindres-carrés linéaire Presque chaque soir, je fais une nouvelle édition du tableau, qu’il est facile d’améliorer n’importe où. Contre la monotonie du travail d’arpentage, c’est toujours une plaisante distraction ; on peut aussi voir immédiatement si quelque chose de douteux s’est glissé, ce qui reste à obtenir, etc. Je vous recommande cette méthode comme modèle. Il ne vous arrivera presque plus jamais de pratiquer une élimination directe, du moins quand vous avez plus de deux inconnues. La procédure indirecte peut se faire à moitié endormi, ou en pensant à autre chose. C.F. Gauss, extrait d’une lettre à G.L. Gerling, datée du 26 décembre 1823, dans laquelle il loue les mérites de sa méthode de relaxation, par rapport aux méthodes directes, pour résoudre l’équation normale d’un problème de moindres-carrés linéaire. Traduction de A. Michel-Pajus.

Définition du problème

 Un problème de moindres-carrés linéaire (MCL) est un problème d’optimisation qui s’écrit de la manière suivante : min x∈Rn  f(x) = 1 2 kAx − bk 2 2  , (19.1) où k · k2 est la norme ℓ2, A est une matrice de type m × n et b ∈ R m. Lorsque m = n, on parle de problème d’interpolation linéaire. Ce problème peut se voir de différentes manières.Le premier point de vue est « concret ». On cherche à déterminer des paramètres x ∈ R n d’un « système » au moyen de mesures b ∈ R m réalisées sur celui-ci. La loi qui relie les paramètres x aux quantités mesurées Ax est supposée linéaire. Le problème de moindres-carrés linéaire permet de déterminer les paramètres x qui donnent des quantités mesurées Ax au plus proche des mesures b. De ce point de vue, le problème consiste à projeter b sur l’image de A au moyen du produit scalaire euclidien (voir la section 2.5.1, on pourrait prendre d’autres produits scalaires d’ailleurs). Le second point de vue est « abstrait ». On s’intéresse à la résolution du système linéaire Ax = b. On n’impose pas que b ∈ R(A), si bien que ce système n’a peut-être pas de solution. Le problème de moindres-carrés linéaire cherche alors à résoudre ce système linéaire « au mieux », en minimisant le résidu Ax − b. Il s’agit donc d’un problème « fondamental », auquel sont rattachés divers concepts bien connus en algèbre linéaire. On note r = rg A le rang de A

Exemple : régression linéaire

En statistiques, un problème de régression linéaire consiste à déterminer au mieux la dépendance affine supposée entre m mesures et n−1 variables. En général m ≫ n. Pour prendre les notations de la section 19.1.1, la i-ième mesure bi résulte des valeurs (ai,1, . . . , ai,n−1) données au vecteur dont elle dépend. Si la mesure dépend affinement du vecteur, il doit y avoir un vecteur de paramètres x ∈ R n , inconnu mais indépendant des mesures prises, tel que le i-ième résidu ri = bi −   nX−1 j=1 ai,jxj + xn   soit nul (xn est le coefficient indépendant des mesures rendant la dépendance affine plutôt que linéaire). Si l’on note b =   b1 . . . bm   ∈ R m et A =   a1,1 · · · a1,n−1 1 . . . . . . . . . am,1 · · · a1,n−1 1   ∈ R m×n , le résidu r = b − Ax ∈ R m doit être nul. En présence d’erreurs de mesure ou d’une relation mesures-données non linéaire, on peut chercher à annuler r au mieux, ce qui conduit à résoudre le problème de moindres-carrés linéaire (19.1). Dans le cas où n = 2, les couples {(bi , ai,1)} m i=1 représentent un nuage de m points dans le plan. Pour la solution x ∈ R 2 du problème de moindres-carrés, l’application affine a ∈ R 7→ ax1 + x2 donne la relation affine « la plus proche » des couples {(bi , ai,1)} m i=1. La figure 19.1 reprise de [619 ; 2013, figure 2.2] donne ainsi la droite passant au plus près, au sens des moindres-carrés, des mesures de l’indice des problèmes sanitaires et sociaux en fonction des inégalités de revenus ; il y a une mesure (un point) pour chacun des 21 pays considérés dans l’étude.

L’ensemble des solutions 

La condition d’optimalité du premier ordre de ce problème s’écrit ∇f(x) = 0 ou encore A TAx = A T b. (19.2) Celle-ci porte le nom d’équation normale de (19.1). La proposition suivante règle la question de l’existence et de l’unicité des solutions de (19.1). Proposition 19.1 (ensemble des solutions) Le problème (19.1) est convexe et admet toujours une solution. Celle-ci est unique si, et seulement si, A est injective. L’ensemble des solutions de (19.1) s’écrit xp + N (A), où xp est une solution particulière de (19.1) et la valeur optimale est 1 2 kP bk 2 2 , où P est le projecteur orthogonal (pour le produit scalaire euclidien) sur R(A) ⊥. Démonstration. Le problème (19.1) consiste à projeter b sur l’image de A (voir la section 2.5.1), qui est un convexe fermé non vide. Il existe donc un unique élément y ∈ R(A) qui est le plus proche de b (proposition 2.25). Cet élément est de la forme y = Ax, où x est une solution de (19.1). Si A est injective f est strictement convexe (ATA est définie positive) et donc (19.1) a une solution unique. Si A n’est pas injective et x est solution, tous les points de x + N (A) (6= {x}) sont aussi solutions ; par ailleurs, si x et x ′ sont deux solutions de (19.1), elles vérifient l’équation normale et donc x − x ′ ∈ N (ATA) = N (A). Enfin, Q = I − P étant le projecteur orthogonal sur R(A), la valeur optimale s’écrit 1 2 kQb − bk 2 2 = 1 2 kP bk 2 2 . ✷ Il existe de nombreuses démonstrations de l’existence d’une solution de (19.1). Celle donnée ci-dessus considère directement le problème d’optimisation. On peut aussi s’intéresser à son système d’optimalité (19.2) (du fait de la convexité du critère — sa hessienne ATA est semi-défini positif — il y a équivalence entre les solutions (19.1) et celles de son système d’optimalité ; théorème 4.9). Pour montrer que ce dernier a toujours une solution, il suffit d’observer que ATb ∈ R(ATA). L’ensemble des solutions de (19.1) sont donc les solutions de l’équation normale (19.2). On peut s’intéresser à la solution de norme minimale, qui est donc définie par  min 1 2 kxk 2 2 ATAx = ATb. (19.3) On peut bien parler de « la » solution de norme minimale, car le critère de ce problème étant strictement convexe, il y a exactement une unique solution de norme minimale. On peut caractériser la solution de ce problème. Comme celui-ci est convexe, ses conditions d’optimalité du premier ordre sont nécessaires et suffisantes (on notera en effet que les contraintes sont qualifiées, car affines). Donc xˆ est la solution de norme minimale si, et seulement si, il existe un multiplicateur λ ∈ R n (non nécessairement unique) tel que  xˆ + ATAλ = 0 ATAxˆ = ATb. (19.4) La première condition s’écrit aussi xˆ ∈ R(ATA) = R(AT). Dès lors, xˆ est solution de norme minimale de (19.1), c’est-à-dire solution de (19.3), si, et seulement si,  xˆ ∈ R(AT) ATAxˆ = ATb. (19.5) Comme xˆ est univoquement déterminé par le système linéaire d’optimalité (19.4) (il peut cependant avoir de multiples solutions en λ si A n’est pas injective), qui se récrit  I ATA ATA 0  xˆ λ  =  0 ATb  , l’application qui à b ∈ R m fait correspondre xˆ est linéaire et la matrice qui la représente est appelée le pseudo-inverse (de Moore-Penrose) de A. On note cette matrice A† et la solution de norme minimale s’écrit xˆ = A † b. On peut voir xˆ comme le résultat d’une tentative de trouver une solution au système Ax = b, qui n’en n’a pas nécessairement (on se rappelle, en particulier, que A n’est pas carrée), au moyen de deux opérations : la première force l’existence d’une solution en relaxant « Ax = b » en un problème de minimisation, celui que l’on trouve dans (19.1), la seconde force l’unicité de la solution par le problème d’optimisation (19.3)

Résolution numérique

On distingue les algorithmes qui utilisent l’équation normale (19.2) (résolution par factorisation de Cholesky ou par gradient conjugué) et ceux qui s’attaquent directement à la formulation originale (19.1) (résolution par factorisation QR ou par factorisation en valeurs singulières). Résolution de l’équation normale par factorisation L’approche la plus simple est de faire la factorisation de Cholesky de ATA et de calculer la solution de l’équation normale en résolvant les deux systèmes triangulaires qui en découlent. Cette approche n’est pas sans inconvénients. D’une part elle oblige de former la matrice ATA, ce qui demande O(mn2 ) opérations (ce n’est souvent pas négligeable). Ensuite, on perd aussi en précision, du fait de l’annulation de l’influence des petits éléments de A (en particulier dans les termes diagonaux (ATA)ii = P j A2 ji). Enfin, cette approche peut éventuellement détruire la creusité éventuelle de A (si A a une ligne pleine, ATA sera en général une matrice pleine). Une autre possibilité, bien adaptée aux matrices A creuses est de résoudre ce que l’on appelle le système linéaire augmenté, équivalent à l’équation normale, que l’on obtient à partir de celle-ci en posant y = −Ax :  I A AT 0  y x  =  0 −ATb  . La matrice K de ce système linéaire est symétrique, mais n’est pas définie positive. D’autre part, l’ordre n + m de K peut être beaucoup plus important que l’ordre n de l’équation normale ; mais si A est creuse, K l’est aussi. On peut alors la factoriser par des méthodes pouvant prendre en compte son caractère creux (comme les solveurs MA27/MA47 de Duff et Reid [181, 182, 183]). On se rappellera que K n’est en général pas définie positive (si A est injective, K a n valeurs propres strictement négatives et m valeurs propres strictement positives), si bien qu’une factorisation de Cholesky ne convient pas.