Python pour le calcul scientifique/Régression et optimisation

Rappelons que dorénavant les programmes commencent tous par :

#!/usr/bin/python3

import numpy as np
import matplotlib.pyplot as plt

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

# **************
# * Paramètres *
# **************

nbpts = 20 # nombre de points

a = 5. # paramètres de la droite
b = 2.
var = 0.5 # variance du bruit

# *************
# * Fonctions *
# *************

def droite(A, x):
    return A[0]*x + A[1]

# *******************
# * Nuage de points *
# *******************

x = np.linspace(0, 1, nbpts)

y = droite((a, b), x) + var*np.random.randn(nbpts)

# **************
# * Résolution *
# **************

X = np.vstack([x, np.ones(len(x))]).T
# X = [[x1, 1],
#      [x2, 1],
#      ...
#      [xn, 1]]

resultat = np.linalg.lstsq(X, y)

aopt, bopt = resultat[0]
erreur = resultat[1][0]/nbpts

# ***************************
# * Affichage des résultats *
# ***************************

plt.plot(x, y, "b+") # nuage de points
plt.plot(x, droite((aopt, bopt), x)) # droite de régression

print("a =", aopt.round(2), ", b =", bopt.round(2),
      "; χ² =", erreur.round(2))

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 :

# **************
# * Résolution *
# **************

X = x.reshape(nbpts, 1)
# X = [[x1],
#      [x2],
#      ...
#      [xn]]

resultat = np.linalg.lstsq(X, y)

[…]

print("a =", aopt.round(2), "; χ² =", erreur.round(2))

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₀ + a₁x + a₂x² + … + a_nxⁿ 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 :

import numpy.polynomial.polynomial as nppol

# **************
# * Résolution *
# **************

A, param = nppol.polyfit(x, y, 1, full=True)

aopt = A[1]
bopt = A[0]
erreur = param[0]/nbpts

Code complet

import numpy as np
import numpy.polynomial.polynomial as nppol
import matplotlib.pyplot as plt

# **************
# * Paramètres *
# **************

nbpts = 20 # nombre de points

a = 5. # paramètres de la droite
b = 2.
var = 0.5 # variance du bruit

# *************
# * Fonctions *
# *************

def droite(A, x):
    return A[0]*x + A[1]

# *******************
# * Nuage de points *
# *******************

x = np.linspace(0, 1, nbpts)

y = droite((a, b), x) + var*np.random.randn(nbpts)

# **************
# * Résolution *
# **************

A, param = nppol.polyfit(x, y, 1, full=True)

aopt = A[1]
bopt = A[0]
erreur = param[0]/nbpts

# ***************************
# * Affichage des résultats *
# ***************************

plt.plot(x, y, "b+") # nuage de points
plt.plot(x, droite((aopt, bopt), x)) # droite de régression

print("a =", aopt.round(2), ", b =", bopt.round(2),
      "; χ² =", erreur.round(2))

 

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¹ et pas le terme en x⁰, nous indiquons donc [1]. Le code devient donc :

A, param = nppol.polyfit(x, y, [1], full=True)

Pour plus de détails voir : Python pour le calcul scientifique/Polynômes#Régression polynomiale.

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 :

import scipy.optimize as opt

# **************
# * Paramètres *
# **************

A0 = (1, 0) # paramètres initiaux

# *************
# * Fonctions *
# *************

def residus(A, x, y):
    return y - droite(A, x)

# **************
# * Résolution *
# **************

resultat = opt.least_squares(residus, A0, args=(x, y))

Aopt = resultat.x # paramètres à l'optimum
aopt = Aopt[0]
bopt = Aopt[1]
erreur = 2*resultat.cost/nbpts

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

import numpy as np
import matplotlib.pyplot as plt
import scipy.optimize as opt

# **************
# * Paramètres *
# **************

nbpts = 20 # nombre de points

a = 5. # paramètres de la droite
b = 2.
var = 0.5 # variance

A0 = (1, 0) # paramètres initiaux

# *************
# * Fonctions *
# *************

def droite(A, x):
    return A[0]*x + A[1]

def residus(A, x, y):
    return y - droite(A, x)

# *******************
# * Nuage de points *
# *******************

x = np.linspace(0, 1, nbpts)

y = droite((a, b), x) + var*np.random.randn(nbpts)

# **************
# * Résolution *
# **************

resultat = opt.least_squares(residus, A0, args=(x, y))

Aopt = resultat.x # paramètres à l'optimum
aopt = Aopt[0]
bopt = Aopt[1]
erreur = 2*resultat.cost/nbpts

# ***************************
# * Affichage des résultats *
# ***************************

plt.plot(x, y, "b+") # nuage de points
plt.plot(x, droite((aopt, bopt), x)) # droite de régression

print("a =", aopt.round(2), "; b =", bopt.round(2),
      "; χ² =", erreur.round(2))

 

Si nous voulons que la loi passe par l'origine (b = 0), il suffit de définir une loi purement linéaire :

def droite(x, A):
    return A*x

def residus(A, x, y):
    return y - droite(x, A)

[…]

Aopt = opt.least_squares(residus, A0, args=(x, y))

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₀ = (a⁰₀, a⁰₁, …, a⁰_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 :

import scipy.optimize as opt

# **************
# * Paramètres *
# **************

a00 = 1 # paramètres initiaux
a01 = 0

# *************
# * Fonctions *
# *************

def droite(x, a, b):
    return a*x + b

# **************
# * Résolution *
# **************

Aopt, Acov = opt.curve_fit(droite, x, y, p0=(a00, a01))

aopt = Aopt[0]
bopt = Aopt[1]

res = y - droite(x, aopt, bopt)
erreur = (res*res).sum()/nbpts

Code complet

import numpy as np
import matplotlib.pyplot as plt
import scipy.optimize as opt

# **************
# * Paramètres *
# **************

nbpts = 20 # nombre de points

a = 5. # paramètres de la droite
b = 2.
var = 0.5 # variance

a00 = 1 # paramètres initiaux
a01 = 0

# *************
# * Fonctions *
# *************

def droite(x, a, b):
    return a*x + b

# *******************
# * Nuage de points *
# *******************

x = np.linspace(0, 1, nbpts)

y = droite(x, a, b) + var*np.random.randn(nbpts)

# **************
# * Résolution *
# **************

Aopt, Acov = opt.curve_fit(droite, x, y, p0=(a00, a01))

aopt = Aopt[0]
bopt = Aopt[1]

res = y - droite(x, aopt, bopt)
erreur = (res*res).sum()/nbpts

# ***************************
# * Affichage des résultats *
# ***************************

plt.plot(x, y, "b+") # nuage de points
plt.plot(x, droite(x, aopt, bopt)) # droite de régression

print("a =", aopt.round(2), "; b =", bopt.round(2),
      "; χ² =", erreur.round(2))

 

Là encore, si nous voulons une loi passant par l'origine, le code devient :

def droite(x, A):
    return A*x

[…]

Aopt, Acov = opt.curve_fit(droite, x, y, p0=A0)

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 :

ƒ : ℝⁿ → ℝ ou bien ℂⁿ → ℂ.

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 : ℝ² → ℝ

(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 :

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# **************
# * Paramètres *
# **************

nbpts = 5 # nombre de points par axe

a = 5. # paramètres du plan
b = 3.
c = 2.
sigma = 0.5 # écart type

# *************
# * Fonctions *
# *************

def plan(x, y, A):
    return A[0]*x + A[1]*y + A[2]

# *******************
# * Nuage de points *
# *******************

x = np.linspace(0, 1, nbpts)
y = np.linspace(0, 1, nbpts)

grilleX, grilleY = np.meshgrid(x, y)
X = grilleX.flatten().reshape(nbpts*nbpts, 1)
Y = grilleY.flatten().reshape(nbpts*nbpts, 1)

Z = plan(X, Y, (a, b, c)) + sigma*np.random.randn(nbpts*nbpts, 1)

# **************
# * Résolution *
# **************

points = np.hstack((X, Y, np.ones_like(X)))
# points = [[x1, y1, 1]
#           [x2, y2, 1]
#            ...
#           [xn, yn, 1]]

resultat = np.linalg.lstsq(points, Z)

aopt, bopt, copt = resultat[0]
erreur = resultat[1][0]/(nbpts*nbpts)

fig = plt.figure()
ax = fig.add_subplot(1, 1, 1, projection="3d")
ax.set_title("np.linalg.lstsq()")

ax.scatter(X, Y, Z, color="r")
ax.plot_surface(grilleX, grilleY, plan(grilleX, grilleY, (aopt, bopt, copt)), alpha=0.8)
print("""*****************
np.linalg.lstsq()""")
print("a =", aopt.round(2), ", b =", bopt.round(2), "; c =", copt.round(2),
      "; χ² =", erreur.round(2))

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.

import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize

def fonction(x):
    '''fonction présentant un minimum proche de zéro'''
    
    return x*x + 5*np.sin(x)

x = np.linspace(-2*np.pi, 2*np.pi, 40)
y = fonction(x)

minimum = optimize.minimize_scalar(fonction)
xopt = minimum.x
ymin = minimum.fun

print(minimum)
print(round(xopt, 2), round(ymin, 2))

plt.plot(x, y)
plt.plot([xopt], [ymin], "+k")

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

import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize

def fonction(x):
    '''fonction présentant un minimum proche de zéro'''
    
    return x*x + 5*np.sin(x)

x = np.linspace(-2*np.pi, 2*np.pi, 40)
y = fonction(x)

minimum = optimize.minimize(fonction, 0)
xopt = minimum.x[0]
ymin = minimum.fun

print(minimum)
print(round(xopt, 2), round(ymin, 2))

plt.plot(x, y)
plt.plot([xopt], [ymin], "+k")

 

Notes et références

---

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