SciPy, optimisation

SciPy (prononcé sai pie) est un ensemble de procédures mathématiques appliquées basées sur l'extension Numpy Python. Avec SciPy, une session Python interactive se transforme en le même environnement complet de traitement de données et de prototypage pour des systèmes complexes comme MATLAB, IDL, Octave, R-Lab et SciLab. Aujourd'hui, je veux 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 en utilisant la commande help () ou en utilisant Shift + Tab.

Présentation

Afin de vous éviter, ainsi que les lecteurs, de rechercher et de lire la source, les liens vers les descriptions de méthodes seront principalement vers Wikipedia. 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, nous suivons les liens vers des publications plus fiables, qui se trouvent à la fin de chaque article ou dans votre moteur de recherche préféré.

Ainsi, le module scipy.optimize inclut l'implémentation des procédures suivantes:

1. Minimisation conditionnelle et inconditionnelle des fonctions scalaires de plusieurs variables (minim) à l'aide de différents algorithmes (Nelder-Mead simplex, BFGS, conjugué des gradients de Newton, COBYLA et SLSQP )
2. Optimisation globale (ex: bassinhopping , diff_evolution )
3. Minimisation des résidus des moindres carrés ( moindres carrés ) et algorithmes d'ajustement des courbes aux moindres carrés non linéaires (curve_fit)
4. Minimiser les fonctions scalaires d'une variable (minim_scalar) et trouver des racines (root_scalar)
5. Solveurs multidimensionnels du système d'équations (racine) utilisant divers algorithmes (Powell hybride, Levenberg-Marquardt ou méthodes à grande échelle, comme Newton-Krylov ).

Dans cet article, nous ne considérerons que le premier élément de cette liste entière.

Minimisation inconditionnelle de la fonction scalaire de plusieurs variables

La fonction minim du package scipy.optimize fournit une interface commune pour résoudre les problèmes de minimisation conditionnelle et inconditionnelle des fonctions scalaires de plusieurs variables. Pour démontrer son travail, 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 (les première et deuxième dérivées, respectivement) sont déjà définies dans le paquet scipy.optimize, nous la définissons 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é, nous dessinons en 3D les valeurs de la fonction de Rosenbrock de deux variables.

 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 pour $$$x_i = 1$ , examinez des exemples de détermination de la valeur minimale de la fonction Rosenbrock à l'aide de diverses procédures scipy.optimize.

La méthode Nelder-Mead Simplex (Nelder-Mead)

Soit un point initial x0 dans un espace à 5 dimensions. Trouvez le point minimum de la fonction Rosenbrock le plus proche d'elle en utilisant l'algorithme Nelder-Mead simplex (l'algorithme est spécifié comme la valeur du paramètre de 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 simplex est le moyen le plus simple de minimiser une fonction clairement définie et assez fluide. Il ne nécessite pas le calcul des dérivées d'une fonction, il suffit de ne spécifier que ses valeurs. La méthode Nelder-Mead est un bon choix pour les problèmes de minimisation simples. Cependant, comme il n'utilise pas d'estimations de gradient, il peut prendre plus de temps pour trouver le minimum.

Méthode Powell

Un autre algorithme d'optimisation dans lequel seules les valeurs de fonction sont calculées est la méthode de Powell . 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-Channo (BFGS)

Pour obtenir une convergence plus rapide vers la solution, la procédure BFGS utilise le gradient de la fonction objectif. Le gradient peut être spécifié en tant que fonction ou calculé en utilisant des 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.

Nous trouvons la dérivée de la fonction de Rosenbrock sous la 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 gradient:

 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 der 

La 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)

L'algorithme de gradient conjugué de Newton est une méthode de Newton modifiée.
La méthode de Newton est basée sur l'approximation d'une fonction dans une région locale par un polynôme du deuxième degré:

où $$$\ mathbf {H} \ left (\ mathbf {x} _ {0} \ right)$ est une matrice de dérivées secondes (matrice de Hesse, Hesse).
Si la Hesse est définie positive, alors le minimum local de cette fonction peut être trouvé en égalant le gradient nul de la forme quadratique à zéro. Le résultat est une expression:

La toile de jute inverse est calculée en utilisant la méthode du gradient conjugué. Un exemple d'utilisation de cette méthode pour minimiser la fonction Rosenbrock est donné ci-dessous. Pour utiliser la méthode Newton-CG, vous devez définir une fonction qui évalue la Hesse.
La Hesse de la fonction de Rosenbrock sous forme analytique est égale à:

où $$$i, j \ in \ left [1, N-2 \ right]$ et $$$i, j \ in \ left [0, N-1 \ right]$ déterminer la matrice $$$N \ fois N$ .

Les autres éléments non nuls de la matrice sont égaux à:

Par exemple, dans l'espace à cinq dimensions N = 5, la matrice de Hesse pour la fonction de Rosenbrock a la forme de bande:

Le code qui calcule cette toile de jute avec le code pour minimiser la fonction de Rosenbrock en utilisant 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 du produit de la Hesse et un vecteur arbitraire

Dans les problèmes du monde réel, le calcul et le stockage de toute la matrice de Hesse peuvent nécessiter des ressources importantes en temps et en mémoire. De plus, en fait, il n'est pas nécessaire de spécifier la matrice de Hesse elle-même, car la procédure de minimisation ne nécessite qu'un vecteur égal au produit de la Hesse avec un autre vecteur arbitraire. Ainsi, d'un point de vue informatique, il est bien préférable de déterminer immédiatement la fonction qui renvoie le résultat du produit de la Hesse avec un vecteur arbitraire.

Considérez la fonction hess, qui prend un vecteur de minimisation comme premier argument, et un vecteur arbitraire comme deuxième argument (avec d'autres arguments de la fonction minimisée). Dans notre cas, il n'est pas très difficile de calculer le produit de la fonction Hessienne de Rosenbrock avec un vecteur arbitraire. Si p est un vecteur arbitraire, alors le produit $$$H (x) \ cdot p$ a la forme:

Une fonction qui calcule le produit de la Hesse et 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: 66 

L'algorithme de région de confiance des gradients conjugués (Newton)

Une mauvaise conditionnalité de la matrice de Hesse et des directions de recherche incorrectes peuvent conduire au fait que l'algorithme des gradients de Newton conjugués peut être inefficace. Dans de tels cas, la préférence est donnée à la méthode de la région de confiance des gradients de Newton conjugués.

Exemple de définition de matrice de Hesse:

 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.] 

Un exemple avec la fonction produit de la Hesse et 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 Krylovsky

Comme la méthode trust-ncg, les méthodes de type Krylovsky sont bien adaptées pour résoudre des problèmes à grande échelle, car elles n'utilisent que des produits à matrice vectorielle. Leur essence est de résoudre le problème dans le domaine confidentiel limité par le sous-espace tronqué de Krylov. Pour les tâches incertaines, il est préférable d'utiliser cette méthode, car elle utilise moins d'itérations non linéaires en raison de moins de produits vectoriels matriciels par sous-tâche, par rapport à la méthode trust-ncg. De plus, la solution de la sous-tâche quadratique est plus précise que la méthode trust-ncg.
Exemple de définition de matrice de Hesse:

 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.] 

Un exemple avec la fonction produit de la Hesse et 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 approximatif basé sur la confiance

Toutes les méthodes (Newton-CG, trust-ncg et trust-krylov) sont bien adaptées pour résoudre des tâches à grande échelle (avec des milliers de variables). Cela est dû au fait que l'algorithme sous-jacent des gradients conjugués implique une détermination approximative de la matrice de Hesse inverse. La solution est itérative, sans décomposition explicite de la Hesse. Puisqu'il suffit de déterminer la fonction pour le produit de la Hesse et un vecteur arbitraire, cet algorithme est particulièrement bon pour travailler avec des matrices éparses (bande diagonale). Cela permet de réduire les coûts de mémoire et de gagner du temps.

Dans les problèmes de taille moyenne, les coûts de stockage et de factorisation de la toile de jute ne sont pas critiques. Cela signifie qu'une solution peut être obtenue en moins d'itérations, résolvant les sous-tâches de la zone de confiance presque exactement. Pour cela, certaines équations non linéaires sont résolues de manière itérative pour chaque sous-tâche quadratique. Une telle solution nécessite généralement 3 ou 4 décompositions de la matrice Holets Hessian. En conséquence, la méthode converge en moins d'itérations et nécessite moins de calcul de la fonction objectif que d'autres méthodes implémentées du domaine de confiance. Cet algorithme n'implique que la détermination de la matrice hessienne complète et ne prend pas en charge la possibilité d'utiliser la fonction produit de la Hesse et un vecteur arbitraire.

Un exemple de 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.x 

 Optimization terminated successfully. Current function value: 0.000000 Iterations: 13 Function evaluations: 14 Gradient evaluations: 13 Hessian evaluations: 14 array([1., 1., 1., 1., 1.]) 

Cela, peut-être, demeure. Dans le prochain article, je vais essayer de raconter les plus intéressants sur la minimisation conditionnelle, l'application de la minimisation dans la résolution des problèmes d'approximation, la minimisation de la fonction d'une variable, les minimiseurs arbitraires et la recherche des racines du système d'équations à l'aide du package scipy.optimize

Source: https://docs.scipy.org/doc/scipy/reference/