SciPy (pronunciado sai pie) es un paquete de aplicación matemática basado en la extensión Numpy Python. Con SciPy, su sesión interactiva de Python se convierte en el mismo entorno completo de ciencia de datos y creación de prototipos de sistemas complejos que MATLAB, IDL, Octave, R-Lab y SciLab. Hoy quiero hablar brevemente sobre cómo utilizar algunos algoritmos de optimización conocidos en el paquete scipy.optimize. Siempre se puede obtener ayuda más detallada y actualizada sobre el uso de funciones usando el comando help() o usando Shift+Tab.
introducción
Para evitar que usted y sus lectores tengan que buscar y leer fuentes primarias, los enlaces a las descripciones de los métodos estarán principalmente en Wikipedia. Por 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, siga los enlaces a publicaciones más autorizadas, que puede 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 funciones escalares de varias variables (minim) utilizando varios algoritmos (Nelder-Mead simplex, BFGS, gradientes conjugados de Newton, и )
- Optimización global (por ejemplo: , )
- Minimizar residuos (least_squares) y algoritmos de ajuste de curvas que utilizan mínimos cuadrados no lineales (curve_fit)
- Minimizar funciones escalares de una variable (minim_scalar) y buscar raíces (root_scalar)
- Solucionadores multidimensionales de sistemas de ecuaciones (raíz) utilizando varios algoritmos (híbrido Powell, o métodos a gran escala como ).
En este artículo consideraremos sólo el primer elemento de toda esta lista.
Minimización incondicional de una función escalar de varias variables.
La función minim del paquete scipy.optimize proporciona una interfaz general para resolver problemas de minimización condicional e incondicional de funciones escalares de varias variables. Para demostrar cómo funciona, necesitaremos una función adecuada de varias variables, que minimizaremos de diferentes maneras.
Para estos efectos es perfecta la función de Rosenbrock de N variables, la cual tiene la forma:
A pesar de que la función de Rosenbrock y sus matrices de Jacobi y Hesse (la primera y segunda derivada, respectivamente) ya están definidas en el paquete scipy.optimize, la definiremos 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, dibujemos en 3D los valores de la función de Rosenbrock de dos variables.
Código de dibujo
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 en , veamos ejemplos de cómo determinar el valor mínimo de la función Rosenbrock usando varios procedimientos scipy.optimize.
Método simplex de Nelder-Mead
Sea un punto inicial x0 en un espacio de 5 dimensiones. Encontremos el punto mínimo de la función de Rosenbrock más cercano a él usando el algoritmo (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 sencilla de minimizar una función definida explícitamente y bastante fluida. No requiere calcular las derivadas de una función, basta con especificar sólo sus valores. El método Nelder-Mead es una buena opción para problemas de minimización simples. 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 la función es . Para usarlo, debe establecer método = 'powell' en la función mínima.
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 Broyden-Fletcher-Goldfarb-Shanno (BFGS)
Para obtener una convergencia más rápida a una solución, el procedimiento utiliza el gradiente de la función objetivo. El gradiente se puede especificar como una función o calcularse utilizando diferencias de primer orden. En cualquier caso, el método BFGS normalmente requiere menos llamadas a funciones que el método simplex.
Encontremos la derivada de la función de Rosenbrock en forma analítica:
Esta expresión es válida para las derivadas de todas las variables excepto la primera y la última, que se definen como:
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 derLa función de cálculo de gradiente se especifica como el valor del parámetro jac de la función mínima, 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)
Algoritmo es un método de Newton modificado.
El método de Newton se basa en aproximar una función en un área local mediante un polinomio de segundo grado:
donde es la matriz de segundas derivadas (matriz de Hesse, Hesse).
Si el hessiano es definido 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 será la expresión:
El hessiano inverso se calcula utilizando el método del gradiente conjugado. A continuación se proporciona un ejemplo del uso de este método para minimizar la función de Rosenbrock. Para utilizar el método Newton-CG, debe especificar una función que calcule el hessiano.
El hessiano de la función de Rosenbrock en forma analítica es igual a:
donde и , define la matriz .
Los elementos restantes distintos de cero de la matriz son iguales a:
Por ejemplo, en un espacio de cinco dimensiones N = 5, la matriz de Hesse para la función de Rosenbrock tiene la forma de una banda:
Código que calcula este Hessiano junto con código para minimizar la función de Rosenbrock usando el método del 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 producto del hessiano y un vector arbitrario.
En problemas del mundo real, calcular y almacenar toda la matriz de Hesse puede requerir importantes recursos de tiempo y memoria. En este caso, en realidad no hay necesidad de especificar la matriz de Hesse en sí, porque el procedimiento de minimización requiere sólo un vector igual al producto del hessiano por otro vector arbitrario. Por lo tanto, desde un punto de vista computacional, es mucho preferible definir inmediatamente una función que devuelva el resultado del producto de Hesse con un vector arbitrario.
Considere la función Hess, que toma el vector de minimización como primer argumento y un vector arbitrario como segundo argumento (junto con otros argumentos de la función a minimizar). En nuestro caso, calcular el producto de la función de Hesse de Rosenbrock con un vector arbitrario no es muy difícil. Si p es un vector arbitrario, entonces el producto tiene la forma:
La función que calcula el producto del hessiano y un vector arbitrario se pasa como el valor del argumento hessp a la función de 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: 66Algoritmo de región de confianza de gradiente conjugado (Newton)
Un mal acondicionamiento de la matriz de Hesse y direcciones de búsqueda incorrectas pueden hacer que el algoritmo de gradiente conjugado de Newton sea ineficaz. En tales casos, se da preferencia a (región de confianza) conjugar gradientes de Newton.
Ejemplo con la definición de la matriz 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.]Ejemplo con la función producto del hessiano 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 Krylov
Al igual que el método Trust-ncg, los métodos de tipo Krylov son muy adecuados para resolver problemas a gran escala porque utilizan únicamente productos matriz-vectorial. Su esencia es resolver un problema en una región de confianza limitada por un subespacio de Krylov truncado. Para problemas inciertos, es mejor usar este método, ya que utiliza un número menor de iteraciones no lineales debido al menor número de productos matriz-vector por subproblema, en comparación con el método trust-ncg. Además, la solución al subproblema cuadrático se encuentra con mayor precisión que utilizando el método trust-ncg.
Ejemplo con la definición de la matriz 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.]
Ejemplo con la función producto del hessiano 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 para solución aproximada en la región de confianza.
Todos los métodos (Newton-CG, Trust-ncg y Trust-krylov) son adecuados para resolver problemas a gran escala (con miles de variables). Esto se debe al hecho de que el algoritmo de gradiente conjugado subyacente implica una determinación aproximada de la matriz de Hesse inversa. La solución se encuentra de forma iterativa, sin una expansión explícita del hessiano. Dado que solo necesita definir una función para el producto de un vector de Hesse y un vector arbitrario, este algoritmo es especialmente bueno para trabajar con matrices dispersas (diagonales de banda). Esto proporciona bajos costos de memoria y importantes ahorros de tiempo.
Para problemas de tamaño mediano, el costo de almacenar y factorizar el hessiano no es crítico. Esto significa que es posible obtener una solución en menos iteraciones, resolviendo los subproblemas de la región de confianza casi exactamente. Para ello, se resuelven iterativamente algunas ecuaciones no lineales para cada subproblema cuadrático. Una solución de este tipo suele requerir 3 o 4 descomposiciones de Cholesky de la matriz de Hesse. Como resultado, el método converge en menos iteraciones y requiere menos cálculos de funciones objetivo que otros métodos de regiones de confianza implementados. Este algoritmo solo implica determinar la matriz de Hesse completa y no admite la capacidad de utilizar la función producto de Hesse y un vector arbitrario.
Ejemplo con minimización de 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.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.])Probablemente nos detendremos ahí. En el próximo artículo intentaré contar las cosas más interesantes sobre la minimización condicional, la aplicación de la minimización en la resolución de problemas de aproximación, minimizar una función de una variable, minimizadores arbitrarios y encontrar las raíces de un sistema de ecuaciones usando scipy.optimize. paquete.
Fuente:
Fuente: habr.com