# Python pour le calcul scientifique/Régression et optimisation

Rappelons que dorénavant les programmes commencent tous par :

## Régression
### Régression linéaire

La fonction permettant la régression linéaire par la méthode des moindres carrés s'appelle est np.linalg.lstsq (linear algebra — least squares).
Nous disposons d'un nuage de n points (x$_i$, y$_i$)$_{0≤i≤n–1}$. Nous voulons les modéliser par une loi affine :
    y = a⋅x + b
les paramètres a et b étant à déterminer. On peut représenter le problème sous la forme d'un système d'équations :
    ${\begin{cases}y_{0}=ax_{0}+b+\varepsilon _{0}\\y_{1}=ax_{1}+b+\varepsilon _{1}\\\vdots \\y_{n-1}=ax_{n-1}+b+\varepsilon _{n-1}\end{cases}}$
ce qui peut également s'écrire sous la forme d'une équation matricielle
    Y = XA + E
avec
    $\mathrm {Y} ={\begin{pmatrix}y_{0}\\y_{1}\\\vdots \\y_{n-1}\end{pmatrix}}{\text{ ; }}\mathrm {X} ={\begin{pmatrix}x_{0}&1\\x_{1}&1\\\vdots &\vdots \\x_{n-1}&1\end{pmatrix}}{\text{ ; }}\mathrm {A} ={\begin{pmatrix}a\\b\end{pmatrix}}{\text{ ; }}\mathrm {E} ={\begin{pmatrix}\varepsilon _{0}\\\varepsilon _{1}\\\vdots \\\varepsilon _{n-1}\end{pmatrix}}$
La syntaxe est alors de la forme : resultat = np.linalg.lstsq(X, Y). La fonction renvoie une liste de valeurs :
* resultat[0] contient les valeurs (a, b) ;
* resultat[1] contient la somme des carrés des résidus.
Par exemple
Si nous voulons que la loi passe par l'origine (b = 0), le système d'équations devient :
    ${\begin{cases}y_{0}=ax_{0}+\varepsilon _{0}\\y_{1}=ax_{1}+\varepsilon _{1}\\\vdots \\y_{n-1}=ax_{n-1}+\varepsilon _{n-1}\end{cases}}$
soit
    Y = XA + E
avec
    $\mathrm {X} ={\begin{pmatrix}x_{0}\\x_{1}\\\vdots \\x_{n-1}\end{pmatrix}}{\text{, }}\mathrm {A} =(a){\text{.}}$
Il suffit donc de modifier le code source comme suit :

#### Ressource

    (en) « numpy.linalg.lstsq », sur Numpy and Scipy Documentation (consulté le 17 avril 2019)

### Régression polynomiale

NumPy propose un module polynomial.polynomial. Une des manières de décrire le polynôme P(x) = a$_0$ + a$_1$x + a$_2$x$^2$ + … + a$_n$x$^n$ est d'utiliser un vecteur [a0, a1, …, an]. Pour déterminer le polynôme de degré n le plus proche d'un jeu de données (x, y), nous disposons de la fonction numpy.polynomial.polynomial.polyfit(x, y, n). Pour avoir le résidus, il faut ajouter le paramètre full = True.
Par exemple, avec cet outil, la régression linéaire consiste à faire une régression par un polynôme de degré 1. En reprenant l'exemple précédent, nous pouvons conserver les sections de code * Paramètres *, * Fonctions *, * Nuage de points * et * Affichage des résultats *, et nous ajoutons le code suivant :
Code complet
Si nous voulons avoir une loi qui passe par l'origine, il faut indiquer non pas le degré maximum du polynôme mais la liste des degrés. Comme nous ne voulons que le terme en x$^1$ et pas le terme en x$^0$, nous indiquons donc [1]. Le code devient donc :
    Pour plus de détails voir : Python pour le calcul scientifique/Polynômes.

### Régression non-linéaire

La régression non linéaire nécessite de charger le module optimize de SciPy. Ce module propose alors plusieurs fonctions.
La première est scipy.optimize.least_squares. Cette méthode consiste à minimiser la fonction de résidus, qui est donc une fonction de coûts, par la méthode des moindres carrés ; il faut donc lui donner la fonction de résidus en paramètre. C'est une méthode itérative. Nous définissons donc :
* une fonction residus(A, x, y) où A est la liste des paramètres de la fonction ;
* un jeu de paramètres initial A0 (point de départ des itérations) ;
* un nuage de points (x, y).
La syntaxe est alors resultat = scipy.optimize.least_squares(residus, A0, arg=(x, y)). La fonction renvoie un dictionnaire contenant entre autres :
* resultat.x : paramètres optimisés Aopt ;
* resultat.cost : valeur de la fonction de coûts.
Reprenons l'exemple précédent de la régression linéaire ; nous reprenons les sections * Paramètres *, * Fonctions *, * Nuage de points * et * Affichage des résultats *, et nous ajoutons le code suivant :
Notez que :
* du point de vue de l'utilisation finale, la fonction que l'on cherche est une fonction de x, définie par un jeu de paramètres A ; on peut la noter ƒ$_A$(x) ;
* du point de vue de l'optimisation, on cherche le minimum de la fonction résidus, une fonction de A, définie par un jeu de paramètres (x$_i$, y$_i$)$_{i∈[0;n–1]}$ qui sont les points expérimentaux ; on peut l'écrire résidus$_{(xi,yi)i∈[0;n–1]}$(A) = ∑$_i$[y$_i$ – ƒ$_A$(x$_i$)].
Pour cette raison, lorsque l'on définit la fonction de résidus, c'est impérativement A qui vient en premier : residus(A, *args, **kwargs).
Code complet
Si nous voulons que la loi passe par l'origine (b = 0), il suffit de définir une loi purement linéaire :
Le module optimize propose également la fonction scipy.optimize.curve_fit qui permet d'éviter de définir la fonction de résidus. On lui indique directement la fonction ƒ$_A$(x). Les deux principales différences par rapport à .least_squares sont :
* lorsque l'on définit la fonction ƒ$_A$(x), on indique d'abord x puis ensuite la liste séparée des paramètres : f(x, a0, a1, …, am) ;
* la fonction .curve_fit n'indique pas la somme quadratique des résidus (le khi carré, χ²), il faut la calculer à part.
La syntaxe est Aopt, Acov = scipy.optimize.curve_fit(f, x, y, p0=(a00, a01, …, a0m)) où A$_0$ = (a$^0$$_0$, a$^0$$_1$, …, a$^0$$_m$) est le jeu de paramètres initial.
La variable de sortie A$_{cov}$ est la matrice de covariance du jeu de paramètres ; elle permet de connaître l'erreur standard sur la détermination des paramètres A.
Reprenons l'exemple précédent de la régression linéaire ; nous reprenons les sections * Paramètres *, * Nuage de points * et * Affichage des résultats *, et nous ajoutons le code suivant :
Code complet
Là encore, si nous voulons une loi passant par l'origine, le code devient :

#### Ressources

* (en) « scipy.optimize.least_squares », sur Numpy and Scipy Documentation (consulté le 17 avril 2019)
* (en) « scipy.optimize.curve_fit », sur Numpy and Scipy Documentation (consulté le 17 avril 2019)

### Régression sur un espace de dimension supérieure

Nous considérons maintenant une loi y = ƒ$_A$(x) où x est un vecteur de dimension m et y est un scalaire :
    ƒ : ℝ$^n$ → ℝ ou bien ℂ$^n$ → ℂ.

### Régression linéaire

Le principe est le même, la matrice X est de dimension n × (m + 1) :
    ${\begin{cases}y_{0}=a_{0}x_{0,0}+a_{1}x_{1,0}+\cdots +a_{m-1}x_{m-1,0}+a_{m}+\varepsilon _{0}\\y_{1}=a_{0}x_{0,1}+a_{1}x_{1,1}+\cdots +a_{m-1}x_{m-1,1}+a_{m}+\varepsilon _{1}\\\vdots \\y_{n-1}=a_{0}x_{0,n-1}+a_{1}x_{1,n-1}+\cdots +a_{m-1}x_{m-1,n-1}+a_{m}+\varepsilon _{n-1}\end{cases}}$
    Y = XA + E
avec
    $\mathrm {Y} ={\begin{pmatrix}y_{0}\\y_{1}\\\vdots \\y_{n-1}\end{pmatrix}}{\text{ ; }}\mathrm {X} ={\begin{pmatrix}x_{0,0}&x_{1,0}&\cdots &x_{m-1,0}&1\\x_{0,1}&x_{1,1}&\cdots &x_{m-1,1}&1\\\vdots &\vdots &&\vdots &\vdots \\x_{0,n-1}&x_{1,n-1}&\cdots &x_{m-1,n-1}&1\end{pmatrix}}{\text{ ; }}\mathrm {A} ={\begin{pmatrix}a_{0}\\a_{1}\\\vdots \\a_{m-1}\\a_{m}\end{pmatrix}}{\text{.}}$
Par exemple, si nous considérons un espace de dimension 2 :
    ƒ$_A$ : ℝ$^2$ → ℝ
       (x, y) ↦ z = ƒ$_A$(x, y)
une loi bilinéaire z = ax + by + c décrit un plan, les matrices sont :
    $\mathrm {Y} ={\begin{pmatrix}z_{0}\\z_{1}\\\vdots \\z_{n-1}\end{pmatrix}}{\text{ ; }}\mathrm {X} ={\begin{pmatrix}x_{0}&y_{0}&1\\x_{1}&y_{1}&1\\\vdots &\vdots &\vdots \\x_{n-1}&y_{n-1}&1\end{pmatrix}}{\text{ ; }}\mathrm {A} ={\begin{pmatrix}a\\b\\c\end{pmatrix}}{\text{.}}$
et le code devient :

## Optimisation

L'optimisation consiste à trouver le minimum d'une fonction ƒ — le minimum ou le maximum, il suffit le cas échéant de changer de signe. Si la fonction ƒ est définie sur un ensemble E, nous cherchons donc
    X$_{opt}$ ∈ E tel que ∀ X ∈ E, ƒ(X) ≥ ƒ(X$_{opt}$).
Notez qu'un tel X$_{opt}$ n'existe pas forcément, et qu'il n'est pas forcément unique. Par exemple, la fonction ƒ(x) = x ne possède pas de minimum sur ℝ, et la fonction ƒ(x) = constante possède une infinité de minima.
On peut restreindre E, en particulier par des inégalités, par exemple trouver le minimum de ƒ sur ℝ avec x > 0.

### Cas général

Dans le cas d'une fonction réelle, le module scipy.optimize propose la fonction minimize_scalar. Par défaut, cette fonction utilise la méthode de Brent qui est une méthode itérative. La syntaxe est simplement resultat = optimize.minimise_scalar(f).
Le résultat est un objet de la classeOptimizeResult : c'est un dictionnaire contenant entre autres :
* x : la valeur de x donnant le minimum ;
* success : un booléen valant True si la méthode a convergé, False sinon ;
* fun : valeur de la fonction au minimum ;
* jac, hes : valeur du jacobien et du hessien ;
* nit : nombre d'itérations.
Pour les fonction de plusieurs variables, nous disposons de la fonction minimize. Par défaut, cette fonction utilise la méthode BFGS (méthode de Broyden-Fletcher-Goldfarb-Shanno) qui est également une méthode itérative. Pour cette fonction, il faut indiquer un point de départ x0 :
    optimize.minimize(f, x0)
Comme il s'agit d'un algorithme de descente, l'algorithme peut être piégé dans un minimum local.
Avec l'exemple précédent :
code source

## Notes et références

Interpolation, extrapolation et lissage < ↑ > Algèbre linéaire