Approximation Non-linéaire

On cherche maintenant à fitter des points expérimentaux \((x_i, y_i)\) par une fonction \(y=f(x)\) dépendant non-linéairement de \(m\) coefficients :

\[y = f(a_1,a_2,\dots, a_m;x)\]

Il s'agit ensuite de minimiser la distance entre cette fonction et les points expérimentaux, soit

\[\min_{a_1\dots a_m} F(a_1\dots a_m)\]

avec

\[F(a_1\dots a_m) = \sum_{i=1}^n \left[ y_i - f \left( a_1,a_2,\dots a_m;x_i \right)\right]^2\]

La fonction \(F\) dépend non-linéairement des paramètres \(a_1 ~ \dots ~ a_m\) et un algorithme d'optimisation doit être utilisé. Nous renvoyons le lecteur à un cours sur le sujet pour en acquérir les bases théoriques. Il faut simplement savoir que MATLAB possède des routines d'optimisation utiles pour ce genre de problème.

Si vous ne disposez pas de la toolbox «Optimisation», il vous faudra utiliser la routine de minimisation fournie dans la version MATLAB de base qui s'appelle fminsearch.

Pour cela il faut simplement programmer la fonction \(F\) à minimiser[1]. Nous conseillons d'écrire la fonction de fitting[2] \(f\) dans une sous-fonction séparée.

Exemple

Voici comment procéder dans l'exemple suivant (très classique) : des analyses spectrométriques fournissent des pics que l'on cherche à fitter par des gaussiennes.

Par exemple dans le cas de deux pics :

\[y = B_1 \exp{\left[-\left(\frac{x-x_1}{\sigma_1}\right)^2 \right]}+ B_2 \exp{\left[-\left(\frac{x-x_2}{\sigma_2}\right)^2 \right]}\]

Il y a 6 paramètres \(B_1\), \(B_2\), \(x_1\), \(x_2\), \(\sigma_1\) et \(\sigma_2\) que nous regroupons dans un tableau a avec la correspondance suivante :

\[\begin{array}{cc} a_1 & B_1 \\ a_2 & B_2 \\ a_3 & \sigma_1 \\ a_4 & \sigma_2 \\ a_5 & x_1 \\ a_6 & x_2 \\ \end{array}\]

Notons que dans ce genre de problème, les centres des gaussiennes \(x_1\) et \(x_2\) sont souvent connus et correspondent à des longueurs d'ondes de raies de base, de telle sorte que le nombre de paramètres dans cet exemple serait réduit à 4. De façon générale, il est toujours souhaitable de diminuer le nombre de paramètre \(a_1 ~ \dots ~ a_m\)

Écrivons tout d'abord la fonction de fitting :

function y = fGaussiennes(a, x)
    y = a(1)*exp(- ( (x-a(5))/a(3) ).^2) + a(2)*exp(- ( (x-a(6))/a(4) ).^2) ;
end

Il est important que cette fonction puisse accepter un tableau en entrée, d'où l'usage des points devant les opérateurs puissance. Ensuite on écrit la fonction funResidu représentant le résidu[3] à optimiser :

function res = funResidu(a, xi, yi)
    res = sum ( (yi - fGaussiennes(a,xi)).^2 );
end

Il reste à appeler l'optimiseur dans un programme principal :

% Définition des points expérimentaux
% (en général chargés depuis un fichier)
 
xi = ...
yi = ...
 
% Estimation de la solution. Bien souvent le succès de l'opération
% dépend de cette initialisation. Attention à ne pas rentrer des valeurs
% farfelues. Notez de plus que les paramètres doivent être entrés
% dans l'ordre défini plus haut.
 
a0 = [1 400 5 2 800 10] % correspond à B1, sigma1, x1, B2, sigma2, x2
 
% Appel de l'optimiseur.
% Utilisation d'une fonction anonyme pour rechercher le minimum en a
 
a_sol = fminsearch ( @(a) funResidu(a,xi,yi), a0);
 
% (Facultatif mais en général utile)
% Superposition points expérimentaux et fittés.
 
plot (xi, yi, '*', xi, fGaussiennes(a_sol, xi))

Voilà l'utilisation de base. Notez que l'algorithme utilisé par fminsearch est assez basique et est parfois insuffisant.

Si vous disposez de la toolbox «optimisation» (ce que nous conseillons si vous avez de nombreux problèmes de ce type), un grand nombre de routines d'optimisation vous est fourni. Toutes utilisent le même moteur d'optimisation qui propose de nombreux algorithmes possibles. Ce moteur est décliné sous forme de plusieurs fonctions, chacune étant dédiée à un problème bien particulier. Pour le problème d'estimation de paramètres proposé ci-dessus, la routine dédiée est lsqcurvefit (qui signifie «least squares curve fit»). Elle permet de ne programmer que la fonction \(f\) donnée par \(y= f(a_1, a_2, ~\dots, a_m ; x)\) et forme elle-même la somme des carrés.

Voici ce qu'il faut écrire pour résoudre le problème précédent :

function y = fGaussiennes(a, x)
    y = a(1)*exp(- ( (x-a(5))/a(3) ).^2) + a(2)*exp(- ( (x-a(6))/a(4) ).^2) ;
end

C'est en fait la même que précédemment. Maintenant on peut appeler directement lsqcurvefit :

% Définition des points expérimentaux
% (en général chargés depuis un fichier)
 
xi = ...
yi = ...
 
% Estimation de la solution. Bien souvent le succès de l'opération
% dépend de cette initialisation. Attention à ne pas rentrer des valeurs
% farfelues. Notez de plus que les paramètres doivent être entrés
% dans l'ordre défini plus haut.
 
a0 = [1 400 5 2 800 10] % correspond à B1, sigma1, x1, B2, sigma2, x2
 
% Appel de l'optimiseur.
% On remarque que la fonction à passer est maintenant fGaussiennes et plus funResidu
% et qu'on passe les points expérimentaux directement au solveur.
 
asol = lsqcurvefit ( @fGaussiennes, a0, xi, yi);
 
% (Facultatif mais en général utile)
% Superposition points expérimentaux et fittés.
 
plot (xi, yi, '*', xi, fGaussiennes(asol, xi))

Notons enfin que les caractéristiques du moteur d'optimisation (utilisation du jacobien ou non, algorithme, nombre max d'itérations, tolérance ...) peuvent être modifiées grâce à l'appel de la fonction optimset. Nous renvoyons le lecteur intéressé à la documentation.