Approximation au sens des moindres carrés

Grégory Vial

5 mars 2003

Une manière simple d’approcher une fonction f par une suite de polynômes p n sur
un intervalle I est d’utiliser l’interpolation de Lagrange. Cependant, il n’y a pas nécessai-
rement convergence de la suite (pn ) vers la fonction d’origine quand le degré des polynômes
tend vers l’infini (c’est ce que l’on nomme le phénomène de Runge, voir [2]). Une solu-
tion consiste à fixer le degré et à effectuer une interpolation par morceaux (il y a alors
convergence de (pn ) vers f quand la taille des sous-intervalles tend vers 0).
Il arrive cependant qu’on préfère utiliser une représentation sur l’intervalle I tout entier.
On adopte alors un point de vue différent de l’interpolation : plutôt que d’imposer que p n
et f coı̈ncident en certains nœuds, on demande qu’ils soient proches d’une manière plus
globale. Précisément, soit f une fonction continue sur l’intervalle I, on s’intéressera aux
deux cas suivants.
Exemple 1. On recherche un polynôme p n de degré inférieur ou égal à n tel que la quantité
Z
|f (x) − pn (x)|2 dx soit minimale.
I

Exemple 2. On recherche un polynôme p n de degré inférieur ou égal à n tel que la quantité

N
X
|f (xi ) − pn (xi )|2 soit minimale,
i=0

où (xi )0≤i≤N est une subdivision de I et N ≥ n.

Point de vue théorique

Soit F un espace vectoriel muni d’une forme bilinéaire symétrique positive h·, ·i ; la
semi-norme associée est notée k·k. On considère P un sous-espace de F de dimension finie
et on suppose que k·k est une norme sur P.
Pour f ∈ F , on recherche p ∈ P tel que la quantité kf − pk soit minimale. Les deux
exemples cités plus haut rentrent dans ce cadre :
Z
Exemple 1. F = C (I, R), P = Rn [X] et hf, gic = f (x)g(x) dx.
I
N
X
Exemple 2. F = C (I, R), P = Rn [X] et hf, gid = f (xi )g(xi ), où (xi ) est une subdivi-
i=0
sion de I et N ≥ n.

Proposition 1 Soit f ∈ F . Il existe un unique élément p ∈ P tel que la quantité kf − pk

soit minimale.

1
Démonstration. Introduisons l’application Φ : q 7→ kf − qk 2 = kqk2 − 2hf, qi + kf k2 . La
fonctionelle Φ est continue et vérifie

lim Φ(q) = +∞.

kqk→+∞

Donc Φ admet un minimum sur P (qui est de dimension finie). Tout polynôme p, en lequel
il est atteint, est un point critique de Φ :

∀q ∈ P, hf − p, qi = 0.

Si on suppose que Φ atteint son minimum en deux polynômes p 1 et p2 , alors pour le choix
q = p1 − p2 , on obtient kp1 − p2 k = 0, d’où p1 = p2 , car h·, ·i est définie sur P.

Remarque. L’application Φ est quadratique, cela explique le terme moindres carrés. Si

on ne s’intéresse pas à l’approximation au sens des moindres cubes, par exemple, c’est
parce-que dans ce cas la différentielle de Φ n’est plus linéaire.

On peut préciser la convergence du polynôme d’approximation pour les deux exemples

considérés.

Proposition 2 Soit pn ∈ Rn [X] le polynôme obtenu dans le cadre de l’exemple 1 (resp. 2).
Alors pn converge vers f dans la norme k·k c (resp. semi-norme k·k d ) quand n tend vers
l’infini.

Démonstration. Les deux (semi-)normes k·k c et k·kd sont dominées par la norme uni-
forme : p √
∀g ∈ F, kgkc ≤ |I| kgk∞ et kgkd ≤ N kgk∞ .
Pour toute fonction continue f , il existe une suite q n ∈ Rn [X] convergeant uniformément
vers f . Par définition de pn , pour la norme continue k·kc ,
p
kf − pn kc ≤ kf − qn kc ≤ |I| kf − qn k∞ −→ 0.

Il en est de même pour la norme discrète k·k d .

Mise en œuvre pratique

On s’intéresse dans ce paragraphe au calcul explicite du polynôme d’approximation au
sens des moindres carrés.

Norme continue
Soit f une fonction continue sur l’intervalle I ; on recherche un polynôme p n , de degré
inférieur ou égal à n, tel que
Z
|f − pn |2 dx soit minimale.
I

On a prouvé, dans la démonstration de la proposition 1, la relation de Pythagore

Z Z
(1) ∀q ∈ Rn [X], pn q dx = f q dx.
I I

2
 Il s’agit d’un système
P linéaire carré : soit (e 0 , . . . , en ) une base de Rn [X], si on recherche
pn sous la forme pn = αj ej , alors (1) devient
n Z
X  Z
∀i = 0, . . . n, ei ej dx αj = f ei dx,
I
j=0 | {z } | I {z }
aij bi

qui n’est autre que le système Aα = b, avec des notations naturelles.

Le choix de la base canonique ej = X j conduit à la matrice de Hilbert, qui est très mal
conditionnée (voir [3]), et n’est donc jamais utilisé. Il est préférable de choisir les polynômes
ej tels que la matrice A soit l’identité, c’est-à-dire la famille orthogonale pour le produit
scalaire de L2 (I). Le calcul des composantes du vecteur b nécessite la plupart du temps
une intégration numérique.

Norme discrète
Soit f une fonction continue sur l’intervalle I ; on recherche un polynôme p n , de degré
inférieur ou égal à n, tel que
N
X
(2) |f (xi ) − pn (xi )|2 soit minimale.
i=0

On peut écrire matriciellement de problème de minimisation : soient A la matrice de Van-

Der-Monde de coefficient générique a ij = xji et b le vecteur de composantes bi = f (xi ).
Alors (2) devient

(3) Chercher α ∈ Rn tel que kAα − bk2 soit minimale,

les composantes du vecteur α sont celles du polynôme p n dans la base canonique de Rn [X].
Notons que la matrice A comporte N lignes et n colonnes ; comme N ≥ n, on parle de
matrice sur-déterminée.

Proposition 3 Le problème (3) est équivalent au système linéaire suivant, dit équations
normales :
ATAα = AT b.

Démonstration. Il s’agit simplement d’exprimer que α annule la différentielle de la

fonctionnelle quadratique kAα − bk 22 .
La proposition précédente fournit un procédé de construction du polynôme d’approxi-
mation pn : la matrice ATA étant symétrique définie positive, on peut utiliser les méthodes
de Choleski ou du gradient conjugué. Cependant, les équations normales sont souvent as-
sez mal conditionnées ; on leur préfère une méthode basée sur la décomposition QR de la
matrice A (voir [1]).

Proposition 4 Soit Q, orthogonale, et R, triangulaire supérieure, deux matrices telles

que A = QR. On écrit, par blocs,
   
R1 T d1
R= et Q b = .
0 d2

Alors α est donné par R1 α = d1 .

3
Démonstration. Il suffit de remarquer que la norme euclidienne est stable par transfor-
mation orthogonale.
Remarque. En général, le nombre d’équations N + 1 est beaucoup plus grand que la
dimension n + 1 de l’espace Rn [X] ; le cas particulier n = 1 correspond à la régression
linéaire.

La fonction matlab suivante calcule le polynôme d’approximation de f :

function valp=approx_mc(strf,x,n,t) d1=d(1:n+1);

% Calcule l’évaluation du polyn^
ome R1=R(1:n+1,1:n+1);
% d’approximation de degré n de f p=R1\d1;
% sur l’intervalle x aux abscisses t. valp=polyval(flipud(p),t);
A=[];
for j=0:n -------------
A=[A,x.^j];
end function y=f(x)
[Q,R]=qr(A); y=1./(1+x.^2);
d=Q’*feval(strf,x);

Noter l’utilisation de la commande feval pour l’appel de la fonction f. Il est possible

d’utiliser les équations normales :

p=(A’*A)\(A’*feval(strf,x));

ou encore la résolution de matlab pour le calcul du vecteur p :

p=A\feval(strf,x);

qui résout le système au sens des moindres carrés quand il est sur-déterminé. Il existe aussi
une fonction matlab qui calcule directement les coefficients du polynôme d’approximation
au sens des moindres carrés : polyfit.
À l’aide du script main.m, on a tracé la fonction f et son polynôme d’approximation
au sens des moindes carrés pour n = 10 et on a mis en évidence la convergence

kf − pn kd −−−−−−→ 0.
n→+∞

N=100; % Convergence en n
M=5;
x=linspace(-M,M,N)’; nn=1:30;
t=linspace(-M,M,1000); err=[];
for n=nn
% Tracé pour n=10 px=approx_mc(’f’,x,n,x);
err=[err,norm(f(x)-px,’inf’)];
n=10; end
pt=approx_mc(’f’,x,n,t); figure
close all plot(nn,err)
plot(t,f(t),t,pt) title(’Erreur en norme discrète
title(’f et p pour n=10’) en fonction du degré’)
xlabel(’x’); xlabel(’degré n’)
legend(’f’,’p’); ylabel(’erreur ||f-p||_d’)

4
 Les graphes construits par le programme main.m sont consignés dans la figure 1. On a
choisi la fonction f définie sur [−5, 5] par
1
f (x) = ,
1 + x2
qui est l’exemple habituellement choisi pour illustrer le phénomène de Runge dans l’inter-
polation de Lagrange. En revanche, l’approximation au sens des moindres carrés converge
quand le degré augmente.
f et p pour n=10 Erreur en norme discrète en fonction du degré
1.2 0.8

f
p
0.7
1

0.6
0.8

0.5

erreur ||f-p|| d
0.6

0.4

0.4
0.3

0.2
0.2

0
0.1

-0.2 0
-5 -4 -3 -2 -1 0 1 2 3 4 5 0 5 10 15 20 25 30
x degré n

Fig. 1 – Résultats du script main.m.

Régression linéaire
Le cas particulier n = 1 est particulièrement utile quand on cherche à mettre en
évidence une vitesse de convergence : il permet d’approcher un nuage de points par une
droite “moyenne”. Le script reglin.m illustre son utilisation :

% On clique avec le bouton

% de droite le dernier point. 1
Regression linéaire

figure
0.9
plot([0,1,1,0],[0,0,1,1],’k’)
hold on 0.8

button=1; 0.7

xx=[];yy=[]; 0.6

while (button~=3) 0.5

[x,y,button]=ginput(1);
0.4
plot(x,y,’ro’)
xx=[xx;x]; 0.3

yy=[yy;y]; 0.2

end 0.1

p=polyfit(xx,yy,1); 0
t=0:.01:1; 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

plot(t,p(1)*t+p(2));
title(’Regression linéaire’) Fig. 2 – Résultat du script reglin.m.

5
Références
[1] P. G. Ciarlet. Introduction à l’analyse numérique matricielle et à l’optimisation.
Collection Mathématiques Appliquées pour la Maı̂trise. [Collection of Applied Mathe-
matics for the Master’s Degree]. Masson, Paris 1982.
[2] J.-P. Demailly. Analyse numérique et équations différentielles. PUG, Grenoble 1996.
[3] M. Schatzman. Analyse numérique. InterEditions, Paris 1991. Cours et exercices
pour la licence. [Course and exercises for the bachelor’s degree].