SciPy, optimizaci贸n


SciPy (pronunciado sai pie) es un paquete de aplicaci贸n matem谩tica basado en la extensi贸n Numpy Python. Con SciPy, una sesi贸n interactiva de Python se convierte en el mismo entorno completo de procesamiento de datos y creaci贸n de prototipos para sistemas complejos como MATLAB, IDL, Octave, R-Lab y SciLab. Hoy quiero hablar brevemente sobre c贸mo usar algunos algoritmos de optimizaci贸n bien conocidos en el paquete scipy.optimize. Siempre se puede obtener ayuda m谩s detallada y actualizada sobre el uso de funciones utilizando el comando help () o Shift + Tab.


Introduccion


Para evitar que usted y los lectores busquen y lean la fuente, los enlaces a las descripciones de los m茅todos se dirigir谩n principalmente a Wikipedia. Como regla general, esta informaci贸n es suficiente para comprender los m茅todos en t茅rminos generales y las condiciones para su aplicaci贸n. Para comprender la esencia de los m茅todos matem谩ticos, seguimos los enlaces a publicaciones m谩s autorizadas, que se pueden encontrar al final de cada art铆culo o en su motor de b煤squeda favorito.


Entonces, el m贸dulo scipy.optimize incluye la implementaci贸n de los siguientes procedimientos:


  1. Minimizaci贸n condicional e incondicional de las funciones escalares de varias variables (minim) utilizando varios algoritmos (Nelder-Mead simplex, BFGS, gradientes de Newton conjugados, COBYLA y SLSQP )
  2. Optimizaci贸n global (ej .: basinhopping , diff_evolution )
  3. Minimizaci贸n de residuos de m铆nimos cuadrados ( m铆nimos cuadrados ) y algoritmos para ajustar curvas a m铆nimos cuadrados no lineales (curve_fit)
  4. Minimizar las funciones escalares de una variable (minim_scalar) y encontrar ra铆ces (root_scalar)
  5. Solucionadores multidimensionales del sistema de ecuaciones (ra铆z) utilizando varios algoritmos (Powell h铆brido, Levenberg-Marquardt o m茅todos a gran escala, como Newton-Krylov ).

En este art铆culo consideraremos solo el primer elemento de esta lista completa.


Minimizaci贸n incondicional de la funci贸n escalar de varias variables.


La funci贸n minim del paquete scipy.optimize proporciona una interfaz com煤n para resolver los problemas de minimizaci贸n condicional e incondicional de las funciones escalares de varias variables. Para demostrar su trabajo, necesitaremos una funci贸n adecuada de varias variables, que minimizaremos de diferentes maneras.


Para estos fines, la funci贸n de Rosenbrock de N variables es perfecta, que tiene la forma:


f l e f t ( m a t h b f x r i g h t ) = s u m N - 1 i = 1 [ 100 l e f t ( x i + 1 - x 2 i r i g h t ) 2 + l e f t ( 1 - x i        right)2]


A pesar de que la funci贸n de Rosenbrock y sus matrices de Jacobi y Hessian (la primera y la segunda derivada, respectivamente) ya est谩n definidas en el paquete scipy.optimize, la definimos nosotros mismos.


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) 

Para mayor claridad, dibujamos en 3D los valores de la funci贸n Rosenbrock de dos variables.


C贸digo para renderizar
 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() 


Sabiendo de antemano que el m铆nimo es 0 para xi=1 , considere ejemplos de c贸mo determinar el valor m铆nimo de la funci贸n Rosenbrock usando varios procedimientos scipy.optimize.


El m茅todo simplex de Nelder-Mead (Nelder-Mead)


Deje que haya un punto inicial x0 en el espacio de 5 dimensiones. Encuentre el punto m铆nimo de la funci贸n de Rosenbrock m谩s cercana a ella utilizando el algoritmo simplex de Nelder-Mead (el algoritmo se especifica como el valor del par谩metro del m茅todo):


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

El m茅todo simplex es la forma m谩s f谩cil de minimizar una funci贸n claramente definida y bastante fluida. No requiere el c谩lculo de derivadas de una funci贸n; es suficiente para especificar solo sus valores. El m茅todo Nelder-Mead es una buena opci贸n para problemas simples de minimizaci贸n. Sin embargo, dado que no utiliza estimaciones de gradiente, puede llevar m谩s tiempo encontrar el m铆nimo.


M茅todo Powell


Otro algoritmo de optimizaci贸n en el que solo se calculan los valores de funci贸n es el m茅todo Powell . Para usarlo, debe establecer method = 'powell' en la funci贸n 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.] 

Algoritmo de Broyden-Fletcher-Goldfarb-Channo (BFGS)


Para obtener una convergencia m谩s r谩pida a la soluci贸n, el procedimiento BFGS utiliza el gradiente de la funci贸n objetivo. El gradiente puede especificarse como una funci贸n o calcularse utilizando diferencias de primer orden. En cualquier caso, generalmente el m茅todo BFGS requiere menos llamadas a funciones que el m茅todo simplex.


Encontramos la derivada de la funci贸n de Rosenbrock en la forma anal铆tica:


 frac partialf partialxj= sum limitsNi=1200(xix2i1)( deltai,j2xi1,j)2(1xi1) deltai1,j=


=200(xjx2j1)400xj(xj+1x2j)2(1xj)


Esta expresi贸n es v谩lida para derivadas de todas las variables, excepto la primera y la 煤ltima, que se definen como:


 frac partialf partialx0=400x0(x1x20)2(1x0),


 frac partialf partialxN1=200(xN1x2N2).$


Veamos la funci贸n de Python que calcula este gradiente:


 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 funci贸n de c谩lculo de gradiente se especifica como el valor del par谩metro jac de la funci贸n minim, como se muestra a continuaci贸n.


 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] 

Algoritmo de gradiente conjugado (Newton)


El algoritmo de gradiente conjugado de Newton es un m茅todo de Newton modificado.
El m茅todo de Newton se basa en aproximar una funci贸n en una regi贸n local por un polinomio de segundo grado:


f left( mathbfx right) approxf left( mathbfx0 right)+ nablaf left( mathbfx0 right) cdot left( mathbfx mathbfx0 right)+ frac12 left( mathbfx mathbfx0 right)T mathbfH left( mathbfx0 right) left( mathbfx mathbfx0 right)


donde  mathbfH left( mathbfx0 right) es una matriz de segundas derivadas (matriz de Hesse, Hesse).
Si el hessiano es definitivo positivo, entonces el m铆nimo local de esta funci贸n se puede encontrar igualando el gradiente cero de la forma cuadr谩tica a cero. El resultado es una expresi贸n:


 mathbfx textrmopt= mathbfx0 mathbfH1 nablaf


La arpillera inversa se calcula utilizando el m茅todo de gradiente conjugado. A continuaci贸n se muestra un ejemplo del uso de este m茅todo para minimizar la funci贸n de Rosenbrock. Para usar el m茅todo Newton-CG, debe definir una funci贸n que eval煤e a Hesse.
La funci贸n de Hesse de Rosenbrock en forma anal铆tica es igual a:


Hi,j= frac partial2f partialxixj=200( deltai,j2xi1 deltai1,j400xi( deltai+1,j2xi deltai,j)400 deltai,j(xi+1x2i)+2 deltai,j=


=(202+1200x2i400xi+1) deltai,j400xi deltai+1,j400xi1 deltai1,j$


donde i,j in left[1,N2 right] y i,j in left[0,N1 right] determinar la matriz N vecesN .


Los restantes elementos distintos de cero de la matriz son iguales a:


 frac partial2f partialx20=1200x20400x1+2


 frac partial2f partialx0x1= frac partial2f partialx1x0=400x0


 frac partial2f partialxN1xN2= frac partial2f partialxN2xN1=400xN2


 frac partial2f partialx2N1=200x


Por ejemplo, en el espacio de cinco dimensiones N = 5, la matriz de Hesse para la funci贸n de Rosenbrock tiene la forma de cinta:


\ tiny \ mathbf {H} = \ begin {bmatrix} 1200x_ {0} ^ {2} -400x_ {1} +2 & -400x_ {0} & 0 & 0 & 0 \\ -400x_ {0} & 202 + 1200x_ {1} ^ {2} -400x_ {2} y -400x_ {1} y 0 y 0 \\ 0 y -400x_ {1} y 202 + 1200x_ {2} ^ {2} -400x_ {3} & -400x_ {2} y 0 \\ 0 y & -400x_ {2} y 202 + 1200x_ {3} ^ {2} -400x_ {4} y -400x_ {3} \\ 0 y 0 y 0 y -400x_ { 3} y 200 \ end {bmatrix}


El c贸digo que calcula este Hessian junto con el c贸digo para minimizar la funci贸n de Rosenbrock usando el m茅todo de gradiente conjugado (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 ejemplo con la definici贸n de la funci贸n del producto de la arpillera y un vector arbitrario


En problemas del mundo real, la computaci贸n y el almacenamiento de toda la matriz de Hesse pueden requerir importantes recursos de tiempo y memoria. Adem谩s, de hecho, no hay necesidad de especificar la matriz de Hesse en s铆, ya que El procedimiento de minimizaci贸n requiere solo un vector igual al producto de la arpillera con otro vector arbitrario. Por lo tanto, desde un punto de vista computacional, es mucho mejor determinar de inmediato la funci贸n que devuelve el resultado del producto de la arpillera con un vector arbitrario.


Considere la funci贸n hess, que toma un vector de minimizaci贸n como primer argumento, y un vector arbitrario como segundo argumento (junto con otros argumentos de la funci贸n minimizada). En nuestro caso, no es muy dif铆cil calcular el producto de la funci贸n Hesse de Rosenbrock con un vector arbitrario. Si p es un vector arbitrario, entonces el producto H(x) cdotp tiene la forma:


 mathbfH left( mathbfx right) mathbfp= beginbmatrix left(1200x20400x1+2 right)p0400x0p1 vdots400xi1pi1+ left(202+1200x2i400xi+1 right)pi400xipi+1 vdots- 400 x N - 2 p N - 2 + 200 p N - 1 e n d b m a t r i x . $ $ 


Una funci贸n que calcula el producto de la arpillera y un vector arbitrario se pasa como el valor del argumento hessp a la funci贸n minimizar:


 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 

El algoritmo de regi贸n de confianza de gradientes conjugados (Newton)


La condicionalidad deficiente de la matriz de Hesse y las direcciones de b煤squeda incorrectas pueden conducir al hecho de que el algoritmo de los gradientes de Newton conjugados puede ser ineficiente. En tales casos, se da preferencia al m茅todo de la regi贸n de confianza de los gradientes de Newton conjugados.


Ejemplo de definici贸n de matriz de arpillera:


 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 ejemplo con la funci贸n del producto de la arpillera y un vector arbitrario:


 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茅todos tipo Krylovsky


Al igual que el m茅todo trust-ncg, los m茅todos de tipo Krylovsky son muy adecuados para resolver problemas a gran escala, ya que solo utilizan productos de vectores de matriz. Su esencia est谩 en resolver el problema en el campo confidencial limitado por el subespacio truncado de Krylov. Para tareas inciertas, es mejor usar este m茅todo, ya que usa menos iteraciones no lineales debido a menos productos de matriz de vectores por subtarea, en comparaci贸n con el m茅todo trust-ncg. Adem谩s, la soluci贸n a la subtarea cuadr谩tica es m谩s precisa que el m茅todo trust-ncg.
Ejemplo de definici贸n de matriz de arpillera:


 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 ejemplo con la funci贸n del producto de la arpillera y un vector arbitrario:


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

Algoritmo aproximado basado en la confianza


Todos los m茅todos (Newton-CG, trust-ncg y trust-krylov) son muy adecuados para resolver tareas a gran escala (con miles de variables). Esto se debe al hecho de que el algoritmo subyacente de los gradientes conjugados implica una determinaci贸n aproximada de la matriz de Hesse inversa. La soluci贸n es iterativa, sin descomposici贸n expl铆cita del hessiano. Dado que solo es necesario determinar la funci贸n para el producto de la arpillera y un vector arbitrario, este algoritmo es especialmente bueno para trabajar con matrices dispersas (diagonal de cinta). Esto proporciona bajos costos de memoria y un importante ahorro de tiempo.


En problemas medianos, los costos de almacenar y factorizar el Hessian no son cr铆ticos. Esto significa que se puede obtener una soluci贸n en menos iteraciones, resolviendo las subtareas del 谩rea de confianza casi exactamente. Para esto, algunas ecuaciones no lineales se resuelven iterativamente para cada subtarea cuadr谩tica. Tal soluci贸n generalmente requiere 3 o 4 descomposiciones de la matriz Holets Hessian. Como resultado, el m茅todo converge en menos iteraciones y requiere menos c谩lculo de la funci贸n objetivo que otros m茅todos implementados del dominio de confianza. Este algoritmo implica solo la determinaci贸n de la matriz hessiana completa y no admite la capacidad de utilizar la funci贸n del producto de la arpillera y un vector arbitrario.


Un ejemplo de minimizar la funci贸n de 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.]) 

Esto, quiz谩s, mora. En el pr贸ximo art铆culo tratar茅 de contar lo m谩s interesante sobre la minimizaci贸n condicional, la aplicaci贸n de la minimizaci贸n para resolver problemas de aproximaci贸n, minimizar la funci贸n de una variable, minimizadores arbitrarios y encontrar las ra铆ces del sistema de ecuaciones usando el paquete scipy.optimize.


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

Source: https://habr.com/ru/post/439288/


All Articles