SciPy (pronunciado sai pie) é un paquete de aplicacións matemáticas baseado na extensión Numpy Python. Con SciPy, a túa sesión interactiva de Python convértese no mesmo ambiente completo de ciencia de datos e de prototipado de sistemas complexos que MATLAB, IDL, Octave, R-Lab e SciLab. Hoxe quero falar brevemente de como usar algúns algoritmos de optimización coñecidos no paquete scipy.optimize. Sempre se pode obter axuda máis detallada e actualizada sobre o uso de funcións usando o comando help() ou usando Maiús+Tab.
Introdución
Para aforrar a vostede e aos lectores de buscar e ler fontes primarias, as ligazóns ás descricións dos métodos estarán principalmente na Wikipedia. Como regra xeral, esta información é suficiente para comprender os métodos en termos xerais e as condicións para a súa aplicación. Para comprender a esencia dos métodos matemáticos, siga as ligazóns a publicacións máis autorizadas, que se poden atopar ao final de cada artigo ou no seu buscador favorito.
Así, o módulo scipy.optimize inclúe a implementación dos seguintes procedementos:
- Minimización condicional e incondicional de funcións escalares de varias variables (mínimo) mediante varios algoritmos (Nelder-Mead simplex, BFGS, gradientes conxugados de Newton, и )
- Optimización global (por exemplo: , )
- Minimización de residuos (mínimos_cadrados) e algoritmos de axuste de curvas usando mínimos cadrados non lineais (curve_fit)
- Minimización das funcións escalares dunha variable (minim_scalar) e busca de raíces (root_scalar)
- Resolvedores multidimensionais de sistemas de ecuacións (raíz) usando varios algoritmos (Powell híbrido, ou métodos a gran escala como ).
Neste artigo consideraremos só o primeiro elemento de toda esta lista.
Minimización incondicional dunha función escalar de varias variables
A función minim do paquete scipy.optimize proporciona unha interface xeral para resolver problemas de minimización condicional e incondicional de funcións escalares de varias variables. Para demostrar como funciona, necesitaremos unha función adecuada de varias variables, que minimizaremos de diferentes xeitos.
Para estes efectos, a función de Rosenbrock de N variables é perfecta, que ten a forma:
A pesar de que a función de Rosenbrock e as súas matrices de Jacobi e Hessian (as derivadas primeira e segunda, respectivamente) xa están definidas no paquete scipy.optimize, definirémola nós mesmos.
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 máis claridade, debuxemos en 3D os valores da función de Rosenbrock de dúas variables.
Código de debuxo
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()
Sabendo de antemán que o mínimo é 0 at , vexamos exemplos de como determinar o valor mínimo da función de Rosenbrock usando varios procedementos scipy.optimize.
Método simple de Nelder-Mead
Sexa un punto inicial x0 no espazo de 5 dimensións. Imos atopar o punto mínimo da función de Rosenbrock máis próximo a el usando o algoritmo (o algoritmo especifícase como o valor do parámetro do 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.]O método simplex é o xeito máis sinxelo de minimizar unha función definida explícitamente e bastante suave. Non require calcular as derivadas dunha función; basta con especificar só os seus valores. O método Nelder-Mead é unha boa opción para problemas de minimización sinxelos. Non obstante, como non utiliza estimacións de gradiente, pode tardar máis tempo en atopar o mínimo.
Método de Powell
Outro algoritmo de optimización no que só se calculan os valores das funcións é . Para usalo, cómpre configurar method = 'powell' na 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-Shanno (BFGS).
Para obter unha converxencia máis rápida a unha solución, o procedemento utiliza o gradiente da función obxectivo. O gradiente pódese especificar como función ou calcularse usando diferenzas de primeira orde. En calquera caso, o método BFGS normalmente require menos chamadas de función que o método simplex.
Atopemos a derivada da función de Rosenbrock en forma analítica:
Esta expresión é válida para as derivadas de todas as variables excepto a primeira e a última, que se definen como:
Vexamos a 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 derA función de cálculo do gradiente especifícase como o valor do parámetro jac da función minim, como se mostra 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 conxugado (Newton)
Algoritmo é un método de Newton modificado.
O método de Newton baséase na aproximación dunha función nunha área local mediante un polinomio de segundo grao:
onde é a matriz de segundas derivadas (matriz de Hessian, Hessian).
Se o Hessiano é definido positivo, entón o mínimo local desta función pódese atopar igualando o gradiente cero da forma cuadrática a cero. O resultado será a expresión:
O Hessian inverso calcúlase mediante o método do gradiente conxugado. A continuación dáse un exemplo de uso deste método para minimizar a función de Rosenbrock. Para usar o método Newton-CG, debes especificar unha función que calcule o Hessiano.
O Hessiano da función de Rosenbrock en forma analítica é igual a:
onde и , define a matriz .
Os restantes elementos distintos de cero da matriz son iguais a:
Por exemplo, nun espazo de cinco dimensións N = 5, a matriz de Hesse para a función de Rosenbrock ten a forma dunha banda:
Código que calcula este Hessian xunto co código para minimizar a función de Rosenbrock usando o método de gradiente conxugado (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 exemplo coa definición da función produto do Hesse e dun vector arbitrario
Nos problemas do mundo real, a computación e o almacenamento de toda a matriz de Hesse poden requirir recursos de memoria e tempo significativos. Neste caso, en realidade non hai necesidade de especificar a propia matriz de Hesse, porque o procedemento de minimización require só un vector igual ao produto do Hesse con outro vector arbitrario. Así, desde o punto de vista computacional, é moi preferible definir inmediatamente unha función que devolva o resultado do produto do Hessiano cun vector arbitrario.
Considere a función hess, que toma o vector de minimización como primeiro argumento e un vector arbitrario como segundo argumento (xunto con outros argumentos da función a minimizar). No noso caso, calcular o produto do Hessian da función de Rosenbrock cun vector arbitrario non é moi difícil. Se p é un vector arbitrario, entón o produto parece:
A función que calcula o produto do Hessiano e un vector arbitrario pásase como valor do argumento hessp á 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 rexión de confianza de gradiente conxugado (Newton)
Un mal acondicionamento da matriz de Hesse e as direccións de busca incorrectas poden facer que o algoritmo de gradiente conxugado de Newton sexa ineficaz. Nestes casos, dáselles preferencia (rexión de confianza) conxuga os gradientes de Newton.
Exemplo coa definición da matriz hessiana:
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.]Exemplo coa función produto do Hessiano e 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
Do mesmo xeito que o método trust-ncg, os métodos de tipo Krylov son moi axeitados para resolver problemas a gran escala porque só usan produtos matriciales. A súa esencia é resolver un problema nunha rexión de confianza limitada por un subespazo de Krylov truncado. Para problemas incertos, é mellor utilizar este método, xa que utiliza un número menor de iteracións non lineais debido ao menor número de produtos matriz-vector por subproblema, en comparación co método trust-ncg. Ademais, a solución ao subproblema cuadrático atópase con máis precisión que usando o método trust-ncg.
Exemplo coa definición da matriz hessiana:
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.]
Exemplo coa función produto do Hessiano e 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 a solución aproximada na rexión de confianza
Todos os métodos (Newton-CG, trust-ncg e trust-krylov) son moi axeitados para resolver problemas a gran escala (con miles de variables). Isto débese ao feito de que o algoritmo de gradiente conxugado subxacente implica unha determinación aproximada da matriz de Hesse inversa. A solución encóntrase de forma iterativa, sen expansión explícita do hessiano. Dado que só precisa definir unha función para o produto dun vector Hessiano e un vector arbitrario, este algoritmo é especialmente bo para traballar con matrices dispersas (diagonal de banda). Isto proporciona baixos custos de memoria e un importante aforro de tempo.
Para problemas de tamaño medio, o custo de almacenamento e factorización do Hessian non é crítico. Isto significa que é posible obter unha solución en menos iteracións, resolvendo case exactamente os subproblemas da rexión de confianza. Para iso, algunhas ecuacións non lineais resólvense iterativamente para cada subproblema cuadrático. Tal solución require normalmente 3 ou 4 descomposicións de Cholesky da matriz de Hesse. Como resultado, o método converxe en menos iteracións e require menos cálculos de funcións obxectivas que outros métodos de rexión de confianza implementados. Este algoritmo só implica determinar a matriz hessiana completa e non admite a capacidade de usar a función produto do hessiano e un vector arbitrario.
Exemplo coa minimización da 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 pararemos alí. No seguinte artigo intentarei contar as cousas máis interesantes sobre a minimización condicional, a aplicación da minimización na resolución de problemas de aproximación, a minimización dunha función dunha variable, os minimizadores arbitrarios e a busca das raíces dun sistema de ecuacións mediante o scipy.optimize. paquete.
Fonte:
Fonte: www.habr.com