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:
- 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 )
- Optimización global (ej .: basinhopping , diff_evolution )
- Minimización de residuos de mínimos cuadrados ( mínimos cuadrados ) y algoritmos para ajustar curvas a mínimos cuadrados no lineales (curve_fit)
- Minimizar las funciones escalares de una variable (minim_scalar) y encontrar raíces (root_scalar)
- 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
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(xi−x2i−1)( deltai,j−2xi−1,j)−2(1−xi−1) deltai−1,j=
=200(xj−x2j−1)−400xj(xj+1−x2j)−2(1−xj)
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(x1−x20)−2(1−x0),
frac partialf partialxN−1=200(xN−1−x2N−2).$
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− mathbfH−1 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,j−2xi−1 deltai−1,j−400xi( deltai+1,j−2xi deltai,j)−400 deltai,j(xi+1−x2i)+2 deltai,j=
=(202+1200x2i−400xi+1) deltai,j−400xi deltai+1,j−400xi−1 deltai−1,j$
donde i,j in left[1,N−2 right] y i,j in left[0,N−1 right] determinar la matriz N vecesN .
Los restantes elementos distintos de cero de la matriz son iguales a:
frac partial2f partialx20=1200x20−400x1+2
frac partial2f partialx0x1= frac partial2f partialx1x0=−400x0
frac partial2f partialxN−1xN−2= frac partial2f partialxN−2xN−1=−400xN−2
frac partial2f partialx2N−1=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(1200x20−400x1+2 right)p0−400x0p1 vdots−400xi−1pi−1+ left(202+1200x2i−400xi+1 right)pi−400xipi+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/