SciPy (prononcé sai pie) est un package d'application mathématique basé sur l'extension Numpy Python. Avec SciPy, votre session Python interactive devient le même environnement complet de science des données et de prototypage de systèmes complexes que MATLAB, IDL, Octave, R-Lab et SciLab. Aujourd'hui, je souhaite parler brièvement de la façon d'utiliser certains algorithmes d'optimisation bien connus dans le package scipy.optimize. Une aide plus détaillée et à jour sur l'utilisation des fonctions peut toujours être obtenue à l'aide de la commande help() ou en utilisant Shift+Tab.
introduction
Afin de vous épargner, ainsi qu'aux lecteurs, la recherche et la lecture de sources primaires, les liens vers les descriptions des méthodes se trouveront principalement sur Wikipédia. En règle générale, ces informations sont suffisantes pour comprendre les méthodes en termes généraux et les conditions de leur application. Pour comprendre l'essence des méthodes mathématiques, suivez les liens vers des publications plus faisant autorité, que vous trouverez à la fin de chaque article ou dans votre moteur de recherche préféré.
Ainsi, le module scipy.optimize inclut la mise en œuvre des procédures suivantes :
- Minimisation conditionnelle et inconditionnelle de fonctions scalaires de plusieurs variables (minim) à l'aide de divers algorithmes (simplex de Nelder-Mead, BFGS, gradients conjugués de Newton, и )
- Optimisation globale (par exemple : , )
- Minimiser les résidus (least_squares) et algorithmes d'ajustement de courbe utilisant les moindres carrés non linéaires (curve_fit)
- Minimiser les fonctions scalaires d'une variable (minim_scalar) et rechercher des racines (root_scalar)
- Solveurs multidimensionnels de systèmes d'équations (racine) utilisant divers algorithmes (Powell hybride, ou des méthodes à grande échelle telles que ).
Dans cet article, nous ne considérerons que le premier élément de toute cette liste.
Minimisation inconditionnelle d'une fonction scalaire de plusieurs variables
La fonction minim du package scipy.optimize fournit une interface générale pour résoudre les problèmes de minimisation conditionnelle et inconditionnelle des fonctions scalaires de plusieurs variables. Pour démontrer comment cela fonctionne, nous aurons besoin d’une fonction appropriée de plusieurs variables, que nous minimiserons de différentes manières.
À ces fins, la fonction Rosenbrock de N variables est parfaite, qui a la forme :
Malgré le fait que la fonction de Rosenbrock et ses matrices de Jacobi et de Hesse (respectivement les dérivées première et seconde) soient déjà définies dans le package scipy.optimize, nous la définirons nous-mêmes.
import numpy as np
def rosen(x):
"""The Rosenbrock function"""
return np.sum(100.0*(x[1:]-x[:-1]**2.0)**2.0 + (1-x[:-1])**2.0, axis=0)Pour plus de clarté, dessinons en 3D les valeurs de la fonction Rosenbrock de deux variables.
Code de dessin
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
# Настраиваем 3D график
fig = plt.figure(figsize=[15, 10])
ax = fig.gca(projection='3d')
# Задаем угол обзора
ax.view_init(45, 30)
# Создаем данные для графика
X = np.arange(-2, 2, 0.1)
Y = np.arange(-1, 3, 0.1)
X, Y = np.meshgrid(X, Y)
Z = rosen(np.array([X,Y]))
# Рисуем поверхность
surf = ax.plot_surface(X, Y, Z, cmap=cm.coolwarm)
plt.show()
Sachant à l'avance que le minimum est 0 à , regardons des exemples sur la façon de déterminer la valeur minimale de la fonction Rosenbrock à l'aide de diverses procédures scipy.optimize.
Méthode du simplexe de Nelder-Mead
Soit un point initial x0 dans un espace à 5 dimensions. Trouvons le point minimum de la fonction Rosenbrock le plus proche en utilisant l'algorithme (l'algorithme est spécifié comme valeur du paramètre méthode) :
from scipy.optimize import minimize
x0 = np.array([1.3, 0.7, 0.8, 1.9, 1.2])
res = minimize(rosen, x0, method='nelder-mead',
options={'xtol': 1e-8, 'disp': True})
print(res.x)Optimization terminated successfully.
Current function value: 0.000000
Iterations: 339
Function evaluations: 571
[1. 1. 1. 1. 1.]La méthode simplexe est le moyen le plus simple de minimiser une fonction explicitement définie et assez fluide. Elle ne nécessite pas de calculer les dérivées d'une fonction, il suffit de préciser uniquement ses valeurs. La méthode de Nelder-Mead est un bon choix pour les problèmes simples de minimisation. Cependant, comme il n’utilise pas d’estimations de gradient, la recherche du minimum peut prendre plus de temps.
Méthode Powell
Un autre algorithme d'optimisation dans lequel seules les valeurs de fonction sont calculées est . Pour l'utiliser, vous devez définir method = 'powell' dans la fonction minim.
x0 = np.array([1.3, 0.7, 0.8, 1.9, 1.2])
res = minimize(rosen, x0, method='powell',
options={'xtol': 1e-8, 'disp': True})
print(res.x)Optimization terminated successfully.
Current function value: 0.000000
Iterations: 19
Function evaluations: 1622
[1. 1. 1. 1. 1.]Algorithme de Broyden-Fletcher-Goldfarb-Shanno (BFGS)
Pour obtenir une convergence plus rapide vers une solution, la procédure utilise le gradient de la fonction objectif. Le gradient peut être spécifié sous forme de fonction ou calculé à l'aide de différences de premier ordre. Dans tous les cas, la méthode BFGS nécessite généralement moins d’appels de fonction que la méthode simplex.
Trouvons la dérivée de la fonction de Rosenbrock sous forme analytique :
Cette expression est valable pour les dérivées de toutes les variables sauf la première et la dernière, qui sont définies comme :
Regardons la fonction Python qui calcule ce dégradé :
def rosen_der (x):
xm = x [1: -1]
xm_m1 = x [: - 2]
xm_p1 = x [2:]
der = np.zeros_like (x)
der [1: -1] = 200 * (xm-xm_m1 ** 2) - 400 * (xm_p1 - xm ** 2) * xm - 2 * (1-xm)
der [0] = -400 * x [0] * (x [1] -x [0] ** 2) - 2 * (1-x [0])
der [-1] = 200 * (x [-1] -x [-2] ** 2)
return derLa fonction de calcul du gradient est spécifiée comme la valeur du paramètre jac de la fonction minim, comme indiqué ci-dessous.
res = minimize(rosen, x0, method='BFGS', jac=rosen_der, options={'disp': True})
print(res.x)Optimization terminated successfully.
Current function value: 0.000000
Iterations: 25
Function evaluations: 30
Gradient evaluations: 30
[1.00000004 1.0000001 1.00000021 1.00000044 1.00000092]Algorithme de gradient conjugué (Newton)
Algorithme est une méthode de Newton modifiée.
La méthode de Newton est basée sur l'approximation d'une fonction dans une zone locale par un polynôme du deuxième degré :
où est la matrice des dérivées secondes (matrice hessienne, hessienne).
Si le Hessien est défini positif, alors le minimum local de cette fonction peut être trouvé en assimilant le gradient nul de la forme quadratique à zéro. Le résultat sera l'expression :
Le Hessian inverse est calculé à l’aide de la méthode du gradient conjugué. Un exemple d'utilisation de cette méthode pour minimiser la fonction de Rosenbrock est donné ci-dessous. Pour utiliser la méthode Newton-CG, vous devez spécifier une fonction qui calcule le Hessian.
La Hessienne de la fonction Rosenbrock sous forme analytique est égale à :
où и , définissez la matrice .
Les éléments non nuls restants de la matrice sont égaux à :
Par exemple, dans un espace à cinq dimensions N = 5, la matrice hessienne de la fonction Rosenbrock a la forme d'une bande :
Code qui calcule ce Hessian ainsi que le code pour minimiser la fonction de Rosenbrock à l'aide de la méthode du gradient conjugué (Newton) :
def rosen_hess(x):
x = np.asarray(x)
H = np.diag(-400*x[:-1],1) - np.diag(400*x[:-1],-1)
diagonal = np.zeros_like(x)
diagonal[0] = 1200*x[0]**2-400*x[1]+2
diagonal[-1] = 200
diagonal[1:-1] = 202 + 1200*x[1:-1]**2 - 400*x[2:]
H = H + np.diag(diagonal)
return H
res = minimize(rosen, x0, method='Newton-CG',
jac=rosen_der, hess=rosen_hess,
options={'xtol': 1e-8, 'disp': True})
print(res.x)Optimization terminated successfully.
Current function value: 0.000000
Iterations: 24
Function evaluations: 33
Gradient evaluations: 56
Hessian evaluations: 24
[1. 1. 1. 0.99999999 0.99999999]Un exemple avec la définition de la fonction produit du Hessian et d'un vecteur arbitraire
Dans les problèmes du monde réel, le calcul et le stockage de l’intégralité de la matrice hessienne peuvent nécessiter beaucoup de temps et de ressources mémoire. Dans ce cas, il n’est en fait pas nécessaire de spécifier la matrice hessienne elle-même, car la procédure de minimisation ne nécessite qu'un vecteur égal au produit du Hessian avec un autre vecteur arbitraire. Ainsi, d'un point de vue informatique, il est de loin préférable de définir immédiatement une fonction qui renvoie le résultat du produit du Hessien avec un vecteur arbitraire.
Considérons la fonction de Hess, qui prend le vecteur de minimisation comme premier argument et un vecteur arbitraire comme deuxième argument (avec d'autres arguments de la fonction à minimiser). Dans notre cas, calculer le produit de la fonction Hessienne de la fonction Rosenbrock avec un vecteur arbitraire n'est pas très difficile. Si p est un vecteur arbitraire, alors le produit ressemble à:
La fonction qui calcule le produit du Hessian et d'un vecteur arbitraire est transmise comme valeur de l'argument hessp à la fonction minimiser :
def rosen_hess_p(x, p):
x = np.asarray(x)
Hp = np.zeros_like(x)
Hp[0] = (1200*x[0]**2 - 400*x[1] + 2)*p[0] - 400*x[0]*p[1]
Hp[1:-1] = -400*x[:-2]*p[:-2]+(202+1200*x[1:-1]**2-400*x[2:])*p[1:-1]
-400*x[1:-1]*p[2:]
Hp[-1] = -400*x[-2]*p[-2] + 200*p[-1]
return Hp
res = minimize(rosen, x0, method='Newton-CG',
jac=rosen_der, hessp=rosen_hess_p,
options={'xtol': 1e-8, 'disp': True})
Optimization terminated successfully.
Current function value: 0.000000
Iterations: 24
Function evaluations: 33
Gradient evaluations: 56
Hessian evaluations: 66Algorithme de région de confiance à gradient conjugué (Newton)
Un mauvais conditionnement de la matrice hessienne et des directions de recherche incorrectes peuvent rendre l'algorithme de gradient conjugué de Newton inefficace. Dans de tels cas, la préférence est donnée à (région de confiance) conjugue les gradients de Newton.
Exemple avec la définition de la matrice Hessienne :
res = minimize(rosen, x0, method='trust-ncg',
jac=rosen_der, hess=rosen_hess,
options={'gtol': 1e-8, 'disp': True})
print(res.x)Optimization terminated successfully.
Current function value: 0.000000
Iterations: 20
Function evaluations: 21
Gradient evaluations: 20
Hessian evaluations: 19
[1. 1. 1. 1. 1.]Exemple avec la fonction produit du Hessien et d'un vecteur arbitraire :
res = minimize(rosen, x0, method='trust-ncg',
jac=rosen_der, hessp=rosen_hess_p,
options={'gtol': 1e-8, 'disp': True})
print(res.x)Optimization terminated successfully.
Current function value: 0.000000
Iterations: 20
Function evaluations: 21
Gradient evaluations: 20
Hessian evaluations: 0
[1. 1. 1. 1. 1.]Méthodes de type Krylov
Comme la méthode trust-ncg, les méthodes de type Krylov sont bien adaptées à la résolution de problèmes à grande échelle car elles utilisent uniquement des produits matrice-vecteur. Leur essence est de résoudre un problème dans une région de confiance limitée par un sous-espace de Krylov tronqué. Pour les problèmes incertains, il est préférable d'utiliser cette méthode, car elle utilise un plus petit nombre d'itérations non linéaires en raison du plus petit nombre de produits matrice-vecteur par sous-problème, par rapport à la méthode trust-ncg. De plus, la solution au sous-problème quadratique est trouvée avec plus de précision qu’en utilisant la méthode trust-ncg.
Exemple avec la définition de la matrice Hessienne :
res = minimize(rosen, x0, method='trust-krylov',
jac=rosen_der, hess=rosen_hess,
options={'gtol': 1e-8, 'disp': True})
Optimization terminated successfully.
Current function value: 0.000000
Iterations: 19
Function evaluations: 20
Gradient evaluations: 20
Hessian evaluations: 18
print(res.x)
[1. 1. 1. 1. 1.]
Exemple avec la fonction produit du Hessien et d'un vecteur arbitraire :
res = minimize(rosen, x0, method='trust-krylov',
jac=rosen_der, hessp=rosen_hess_p,
options={'gtol': 1e-8, 'disp': True})
Optimization terminated successfully.
Current function value: 0.000000
Iterations: 19
Function evaluations: 20
Gradient evaluations: 20
Hessian evaluations: 0
print(res.x)
[1. 1. 1. 1. 1.]
Algorithme de solution approximative dans la région de confiance
Toutes les méthodes (Newton-CG, trust-ncg et trust-krylov) sont bien adaptées à la résolution de problèmes à grande échelle (avec des milliers de variables). Cela est dû au fait que l’algorithme du gradient conjugué sous-jacent implique une détermination approximative de la matrice de Hesse inverse. La solution est trouvée de manière itérative, sans développement explicite du Hessien. Puisqu'il vous suffit de définir une fonction pour le produit d'un hessien et d'un vecteur arbitraire, cet algorithme est particulièrement adapté au travail avec des matrices clairsemées (diagonales de bande). Cela permet de réduire les coûts de mémoire et de gagner beaucoup de temps.
Pour les problèmes de taille moyenne, le coût de stockage et d’affacturage de la toile de jute n’est pas critique. Cela signifie qu’il est possible d’obtenir une solution en moins d’itérations, résolvant presque exactement les sous-problèmes de la région de confiance. Pour ce faire, certaines équations non linéaires sont résolues de manière itérative pour chaque sous-problème quadratique. Une telle solution nécessite généralement 3 ou 4 décompositions de Cholesky de la matrice hessienne. En conséquence, la méthode converge en moins d’itérations et nécessite moins de calculs de fonctions objectives que les autres méthodes de région de confiance mises en œuvre. Cet algorithme implique uniquement la détermination de la matrice hessienne complète et ne prend pas en charge la possibilité d'utiliser la fonction produit de la hessienne et un vecteur arbitraire.
Exemple avec minimisation de la fonction Rosenbrock :
res = minimize(rosen, x0, method='trust-exact',
jac=rosen_der, hess=rosen_hess,
options={'gtol': 1e-8, 'disp': True})
res.xOptimization terminated successfully.
Current function value: 0.000000
Iterations: 13
Function evaluations: 14
Gradient evaluations: 13
Hessian evaluations: 14
array([1., 1., 1., 1., 1.])On s'arrêtera probablement là. Dans le prochain article, j'essaierai de raconter les choses les plus intéressantes sur la minimisation conditionnelle, l'application de la minimisation dans la résolution de problèmes d'approximation, la minimisation d'une fonction d'une variable, les minimisateurs arbitraires et la recherche des racines d'un système d'équations à l'aide de scipy.optimize. emballer.
Source:
Source: habr.com