Ce tuto traite d'ajustement aux moindres carrés de courbes et de surfaces.
1) Contenu proposé :
  1. une introduction aux modèles linéaires
  2. une présentation de l'opérateur \ de MATLAB (mldivide)
  3. un exemple de code complet pour calculer une droite de régression
  4. un exemple de code complet pour ajuster Z = k + a*X + b*Y + c*X^2 + dX*Y + e*Y^2 en connaissant Z, Y et Y
  5. un aperçu du redressement d'images

Il a trouvé sa place dans cette FAQ car l'opérateur \ (anti-division, encore appelée division à gauche, qui correspond à la fonction MATLAB mldivide) reste méconnu de la plupart des utilisateurs MATLAB parmi ceux qui e ont le plus besoin, or c'est un des outils les plus puissants de MATLAB.

2) Problème posé :
Je connais la valeur de Z en une collection de points X.
Z a en général une colonne (1 variable connue) mais peut en avoir plusieurs (ajustement simultané de plusieurs variables sur les mêmes points).
X a une colonne (ajustement de courbe le long d'un axe), deux colonnes (ajustement d'une surface) ou plus (ajustement d'un modèle sur N variables).
Je veux ajuster une droite, ou une surface, ou (presque) n'importe quelle équation sur les valeurs observées.
je veux ensuite pouvoir calculer :
* les valeurs du modèle aux points observés (pour vérifier que le modèle est bon)
* les valeurs du modèle sur des données quelconques, par exemple régulièrement espacées.

3) Exemples d'applications concernées
Entrent en particulier dans cette catégorie les problèmes suivants :
  • régression linéaire simple
  • calcul du coefficient de corrélation r²
  • ajustement de courbe polynomiales
  • redressement d'image

4) pas à pas dans Z = k + a*X + b*Y + c*X^2 + dX*Y + e*Y^2
La première étape est cruciale. Il s'agit d'écrire ce problème sous forme de matrice.
Code : Sélectionner tout - Visualiser dans une fenêtre à part
1
2
3
4
X = X(:), Y=Y(:), Z=Z(:) ; % X, Y et Z doivent être en colonnes
nx = length(X) ; % le nombre de points
U = ones(nx, 1) ; % une colonne de 1
M = [U X Y X.^2 X.*Y Y.^2] ; % M a nx lignes, et autant de colonnes que de paramètres à ajuster, y-compris une colonne de 1 qui correspond au k
Après cette mise en forme, le problème se résume à l'équation suivante d'inconnue K :
Code : Sélectionner tout - Visualiser dans une fenêtre à part
M * K = Z ;
Avec K, un vecteur colonne qui a autant de lignes que M a de colonnes (et que j'ai e coefficients différents à ajuster). ici 5.

5) la méthode classique:
Un peu de théorie n'est pas inutile pour comprendre la suite. Ceux qui sont pressés peuvent sauter cette section.
Ne soyons pas effrayés par le fait que X, K et Z sont des matrices. Faisons comme si de rien n'était : la solution est alors toute simple : je multiple chaque terme de l'égalité par l'inverse d M et le tour est joué.
inv(M) * M * K = inv(M) * Z ;
donc K = inv(M) * Z ;
Note : si M est un nombre réel, son inverse est 1/M. il n'y a pas de magie ici.

Mais il y a quand même un os : M est loin d'être inversible car M n'est même pas une matrice carré !!!
La beauté de la chose, c'est que la solution n'est pas loin : il suffit au préalable de multiplier M par sa transposée. On construit (M' * M) qui est non seulement carrée, mais, de plus, "symétrique définie positive". "symétrique définie positive". (symétrique : elle est égale à sa transposée. définie positives : toutes ses valeurs propres sont >0. Ces matrices ont des propriétés remarquables, en particulier pour leur décomposition.)

en reprenant notre petite équation à résoudre :
M * K = Z ;
M' * M * K = M' * Z ;
inv(M' * M) * (M' * M) * K = inv(M' * M) * M' * Z ;
K = inv(M' * M) * M' * Z ;

6) la méthode avec \ de Matlab
MATLAB fait cela mieux et plus vite en une seule ligne :
Code : Sélectionner tout - Visualiser dans une fenêtre à part
K = M\Z ;
c'est mieux et plus rapide car, comme le dit la doc, il n'est pas nécessaire d'inverser effectivement la matrice M'M pour résoudre le problème. Je ne rentre pas dans les détails. ce sont les conséquences des propriétés remarquables des matrices symétriques définies positives.

Note sur l'écriture M\Z :
"\" est une "division par la gauche"
la logique de cette écriture typique de MATLAB est la suivante :
équation de départ : M*K = Z
division par la gauche des deux membres de l'égalité : M\M*K = M\Z
en considérant que M divisé par lui-même donne l'identité, on tombe sur l'équation finale : K = M\Z.


7) calcul des valeurs estimées et du coefficient de corrélation:
Le coefficient de corrélation est la partie de la variance des données Z expliquée par le modèle.
Code : Sélectionner tout - Visualiser dans une fenêtre à part
1
2
3
4
5
6
Z_modele = M*K ; 
residus = Z-Z_modele ;
variance_d_origine = var(Z) ;
variance_expliquee = var(Z_modele) ;
variance_residuelle = var(residus) ;
coefficient_correlation_r2 = variance_expliquee/variance_d_origine ;

7) application à la droite de correlation y = ax+b:
Code : Sélectionner tout - Visualiser dans une fenêtre à part
1
2
3
4
5
6
7
8
9
10
11
12
function [a, b, r2, Y_modele] = y_egal_ax_plus_b (X, Y)
if ndims(X)>1, error('X doit être un vecteur'), end
if ndims(Y)>1, error('Y doit être un vecteur'), end
if length(X)~=lenngth(Y), error('X et Y doivent avoir la même taille'), end
s = size(X) ;
M = [ones(max(s),1) X(:)] ;
K = M \ Y(:) ;
a = K(2) ;
b = K(1) ;
Y_modele = reshape(M*K, s) ;
r2 = var(Y_modele)/var(Y) ;
return

8) application à l'ajustement d'une courbe polynomiale de degré quelconque sur un plan:
l'ajustement d'une courbe sur un plan se fait exactement de la même manière que la droite de corrélation. On va juste généraliser la construction de la matrice M à un polynôme quelconque de X et de Y. le code ci-dessous est directement fonctionnel

Code : Sélectionner tout - Visualiser dans une fenêtre à part
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
% 1/ préparation des données à simuler
% 1.a) X et Y sont tirés dans [0 1[
taille=6000;
X=rand(taille,1) ;
Y=rand(taille,1) ;
 
% 1.b) on simule Z = X^a * Y^b après avoir ajouté du bruit à X et Y
bruit = 0.1 ;
a=1.75 ;
b=-0.9 ;
Z=(X+bruit*rand(taille,1)).^a .* (Y+bruit*rand(taille,1)).^b ;
 
% 2/ définition du modèle à résoudre
px=13 ; % degré max en X
py=13 ; % degré max en Y
ptot=13 ; % degré max total
[i,j]=meshgrid(0:px,0:py) ; % génération de toutes les possibilités
u=find(i+j<=ptot) ; % écarter les degrés trop élevés
i=i(u) ; j=j(u) ;
nk=length(u) ; % Nombre de paramètres à calculer ?
 
% 3/ Résolution du problème
% 3.a) on calcule 4 matrices de dimension [taille * nk]
nx = length(X) ; % combien de données ?
X=repmat(X,1,nk) ;
Y=repmat(Y,1,nk) ;
i=repmat(i',nx,1) ;
j=repmat(j',nx,1) ;
 
% 3.b) on les assemble dans la matrice M
M=X.^i.*Y.^j ;
 
% 3.c1) résolution normale par le \ de Matlab
k_matlab = M\Z ;
 
% 3.c2) pour comparaison, résolution par inv(X'X) 
k_manuel = inv(M' * M) * M' * Z ;
 
% 3.d) comparaison des coefficients du polynôme par les deux méthodes
% C'est quasiment la même chose, mais pas tout à fait.
% Dans certains cas dégénérés, l'inversion peut être fausse et l'anti-divison continuer à trouver le bon résultat.
disp('Réultats de chaque calcul') ;
disp ([k_matlab k_manuel]) ;
 
disp('différence relative entre les deux méthodes') ;
disp(k_matlab-k_manuel) ;
 
% 4. Calcul des Z simulés
Z_matlab = M*k_matlab ;
Z_manuel = M*k_manuel ;
differences = (Z_matlab-Z_manuel) ;
 
disp ('l''écart type des différences en Z entre chaque version est') ;
disp (std(differences)) ;
 
% 5. calcul du coefficient de corrélation
residus = Z-Z_matlab ;
r2 = var(Z_matlab)/var(Z) ;
disp(sprintf('coefficient de corrélation r2 = %f', r2)) ;
9) application au redressement d'image:
Problème : une collection de points [Xi,Yi] connus sur l'image et je connais leurs coordonnées dans l'image redressée que je veux calculer.
Je veux ajuster un polynôme du deuxième degré en XY pour recalculer X, et un autre polynôme pour ajuster Y.
L'idée est de calculer la matrice M exactement comme sur l'exemple précédent, puis d'ajuster K avec deux colonnes, une pour les X, et une pour les Y.
Regardez bien les petites différences avec le code précédent : elles sont infimes.
Code : Sélectionner tout - Visualiser dans une fenêtre à part
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
% 1/ préparation des données à simuler
% 1.a) X et Y sont tirés dans [1 npix[
npix = 400 ; % taille de l'image à redresser
taille = 10 ; % nombre de points de recalage
 
% [Xw, Yw] sont les coordonnées des points de recalage observés sur l'image
Xw=1+round(npix*rand(taille,1)) ; % les points vus sur l'image
Yw=1+round(npix*rand(taille,1)) ;
imagew = sqrt(((0:1/npix:1)' * (0:1/npix:1))) ; % prendre ici les vraies données
figure(1) ; imshow(imagew) ;
 
% [Xi Yi] sont les coordonnées théoriques de ces points à obtenir sur l'image recalée 
% pour la démo on applique une formule quelconque qui donne une "petite" déformation
force = 1/5 ; % augmenter la force pour augmenter la déformation
a = 1+(rand(1,7)-0.5)*force ;
b = 1+(rand(1,7)-0.5)*force ;
Xi = (a(1) + a(2)*Xw.^a(2) + a(3)*Yw.^a(4) + sqrt(a(5)*Xw.^a(6).*Yw.^a(7)))/3 ; % une petite déformation 
Yi = (b(1) + b(2)*Xw.^b(2) + b(3)*Yw.^b(4) + sqrt(b(5)*Xw.^b(6).*Yw.^b(7)))/3 ; % une petite déformation 
 
% 2/ définition du modèle à résoudre
px=2 ; % degré max en X
py=2 ; % degré max en Y
ptot=2 ; % degré max total
[i,j]=meshgrid(0:px,0:py) ; % génération de toutes les possibilités
u=find(i+j<=ptot) ; % écarter les degrés trop élevés
i=i(u) ; j=j(u) ;
nk=length(u) ; % Nombre de paramètres à calculer
 
% 3/ Résolution du problème
% 3.a) on calcule 4 matrices
nx = length(Xw) ; 
X=repmat(Xw,1,nk) ;
Y=repmat(Yw,1,nk) ;
im=repmat(i',nx,1) ;
jm=repmat(j',nx,1) ;
 
% 3.b) on les assemble dans la matrice M
M=X.^im.*Y.^jm ;
 
% 3.c1) résolution normale par \
% noter que K a deux colonnes. une pour les coefficients e recalage en X, l'autre pour Y
K = M \ [Xi Yi] ;
disp (sprintf ('noter que les dimensions de K sont (%d lignes, %d colonnes)', size(K))) ;
 
XY_model = M*K ; % calculer les points recalés
r2 = (var(XY_model)) ./ var([Xi Yi]) ; % r2 se calcule indépendamment en X et en Y
disp(sprintf('coefficient de correlation r2 = %f en X et %f en Y', r2)) ;
 
% calcul de la nouvelle image i à partir de l'image connue mais déformée w
[Xni, Yni] = meshgrid (1:npix, 1:npix) ;
sn = prod(size(Xni)) ;
Xni = reshape(Xni,sn,1) ;
Yni = reshape(Yni,sn,1) ;
nx = length(Xw) ; % combien de données ?
X=repmat(Xni,1,nk) ;
Y=repmat(Yni,1,nk) ;
im=repmat(i',sn,1) ;
jm=repmat(j',sn,1) ;
M=X.^im.*Y.^jm ;
 
XYnw = M * K ;
XYnw = max(min(round(XYnw),npix),1) ;
imagei = zeros(npix) ;
disp('...'); tic ;
for i=1:sn
imagei(Xni(i), Yni(i)) = imagew(XYnw(i,1), XYnw(i,2)) ;
end
toc ;
figure(2) ; imshow(imagei) ;
Note : Ce programme m'a servi à vous montrer un modèle linéaire avec deux colonnes pour les coefficients K. C'est finalement très simple.
le programme marche bien dans pas mal de cas si on a assez de points de recalage et qu'on prend un polynôme de degré suffisant. je l'ai moi-même utilisé 'pour de vrai'.
Toutefois, si on veut spécifiquement faire quelque chose comme une translation-rotation-mise à l'échelle, ce n'est pas l'outil qu'il faut.

Conclusion:
Pour en savoir plus : opérateur mldivide
Et pour aller plus loin ?

Ça devient tellement facile d'ajuster n'importe quoi sur autre chose qu'à un moment, on a besoin de limites. "Je suis déjà au degré 13 avec 105 paramètres. Je n'ai qu'une ligne à modifier pour passer au degré 14. Jusqu’où s'arrêtera-t-il ?

Il y a pour cela une méthode rigoureuse qui passe par le calcul d'une grandeur appelée F. F représente la quantité de variance expliquée par un modèle meilleur, par rapport à la variance déjà expliquée par un modèle moins bon. Une fois calculé, on regarde dans une table la valeur de F en fonction du nombre de paramètres du modèle et du nombre de données pour savoir si l'éventuel accroissement de r2 a des chances d'être le fruit du hasard.

J'expliquerai ça une autre fois, dans un autre mini tuto...

OL