SciPy, otimização

SciPy (pronuncia-se sai pie) é um pacote de aplicativos matemáticos baseado na extensão Numpy Python. Com o SciPy, uma sessão interativa em Python se transforma no mesmo ambiente completo de processamento e prototipagem de dados para sistemas complexos como MATLAB, IDL, Octave, R-Lab e SciLab. Hoje, quero falar brevemente sobre como usar alguns algoritmos de otimização conhecidos no pacote scipy.optimize. Ajuda mais detalhada e atualizada sobre o uso de funções sempre pode ser obtida usando o comando help () ou Shift + Tab.

1. Introdução

Para evitar que você e os leitores pesquisem e leiam a fonte, os links para as descrições dos métodos serão principalmente para a Wikipedia. Como regra, essas informações são suficientes para entender os métodos em termos gerais e as condições para sua aplicação. Para entender a essência dos métodos matemáticos, seguimos os links para publicações mais autoritativas, que podem ser encontradas no final de cada artigo ou em seu mecanismo de pesquisa favorito.

Portanto, o módulo scipy.optimize inclui a implementação dos seguintes procedimentos:

1. Minimização condicional e incondicional de funções escalares de várias variáveis ​​(minim) usando vários algoritmos (Nelder-Mead simplex, BFGS, gradientes de Newton conjugados, COBYLA e SLSQP )
2. Otimização global (por exemplo: basinhopping , diff_evolution )
3. Minimização de resíduos de mínimos quadrados (least_squares) e algoritmos para ajustar curvas a mínimos quadrados não lineares (curve_fit)
4. Minimizar funções escalares de uma variável (minim_scalar) e encontrar raízes (root_scalar)
5. Solucionadores multidimensionais do sistema de equações (raiz) usando vários algoritmos (híbrido Powell, Levenberg-Marquardt ou métodos de larga escala, como Newton-Krylov ).

Neste artigo, consideraremos apenas o primeiro item de toda a lista.

Minimização incondicional da função escalar de várias variáveis

A função minim do pacote scipy.optimize fornece uma interface comum para resolver os problemas de minimização condicional e incondicional das funções escalares de várias variáveis. Para demonstrar seu trabalho, precisaremos de uma função adequada de várias variáveis, que minimizaremos de maneiras diferentes.

Para esses propósitos, a função Rosenbrock de N variáveis ​​é perfeita, que tem a forma:

Apesar do fato de que a função Rosenbrock e suas matrizes Jacobi e Hessian (a primeira e a segunda derivadas, respectivamente) já estão definidas no pacote scipy.optimize, nós a definimos.

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 maior clareza, desenhamos em 3D os valores da função Rosenbrock de duas variáveis.

 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 antecipadamente que o mínimo é 0 para $$$x_i = 1$ , considere exemplos de como determinar o valor mínimo da função Rosenbrock usando vários procedimentos scipy.optimize.

O método simples Nelder-Mead (Nelder-Mead)

Haja um ponto inicial x0 no espaço 5-dimensional. Encontre o ponto mínimo da função Rosenbrock mais próximo usando o algoritmo simplex Nelder-Mead (o algoritmo é especificado 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 é a maneira mais fácil de minimizar uma função claramente definida e razoavelmente suave. Não requer o cálculo de derivadas de uma função, basta especificar apenas seus valores. O método Nelder-Mead é uma boa opção para problemas simples de minimização. No entanto, como não usa estimativas de gradiente, pode levar mais tempo para encontrar o mínimo.

Método Powell

Outro algoritmo de otimização no qual apenas os valores das funções são calculados é o método Powell . Para usá-lo, você precisa definir method = 'powell' na função 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 obter uma convergência mais rápida para a solução, o procedimento BFGS usa o gradiente da função objetivo. O gradiente pode ser especificado como uma função ou calculado usando diferenças de primeira ordem. Em qualquer caso, geralmente o método BFGS requer menos chamadas de função que o método simplex.

Encontramos a derivada da função Rosenbrock na forma analítica:

Essa expressão é válida para derivadas de todas as variáveis, exceto a primeira e a última, que são definidas como:

Vejamos a função Python que calcula esse 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 

A função de cálculo do gradiente é especificada como o valor do parâmetro jac da função minim, conforme mostrado abaixo.

 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)

O algoritmo de gradiente conjugado de Newton é um método de Newton modificado.
O método de Newton baseia-se na aproximação de uma função em uma região local por um polinômio de segundo grau:

onde $$$\ mathbf {H} \ left (\ mathbf {x} _ {0} \ right)$ é uma matriz de segundas derivadas (matriz de Hessian, Hessian).
Se o Hessiano é definido positivamente, então o mínimo local dessa função pode ser encontrado equiparando o gradiente zero da forma quadrática a zero. O resultado é uma expressão:

O Hessian inverso é calculado usando o método do gradiente conjugado. Um exemplo do uso desse método para minimizar a função Rosenbrock é dado abaixo. Para usar o método Newton-CG, você deve definir uma função que avalie o Hessian.
A função hessiana da Rosenbrock na forma analítica é igual a:

onde $$$i, j \ in \ esquerda [1, N-2 \ direita]$ e $$$i, j \ in \ esquerda [0, N-1 \ direita]$ determinar a matriz $$$N \ vezes N$ .

Os demais elementos diferentes de zero da matriz são iguais a:

Por exemplo, no espaço tridimensional N = 5, a matriz Hessiana da função Rosenbrock tem a forma de fita:

O código que calcula esse Hessian junto com o código para minimizar a função Rosenbrock usando o método do 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] 

Um exemplo com a definição da função do produto do Hessian e um vetor arbitrário

Em problemas do mundo real, calcular e armazenar toda a matriz hessiana pode exigir recursos significativos de tempo e memória. Além disso, de fato, não há necessidade de especificar a própria matriz de Hessian, pois o procedimento de minimização requer apenas um vetor igual ao produto do Hessian com outro vetor arbitrário. Assim, do ponto de vista computacional, é muito preferível determinar imediatamente a função que retorna o resultado do produto do Hessian com um vetor arbitrário.

Considere a função hess, que usa um vetor de minimização como primeiro argumento e um vetor arbitrário como o segundo argumento (junto com outros argumentos da função minimizada). No nosso caso, não é muito difícil calcular o produto da função hessiana de Rosenbrock com um vetor arbitrário. Se p é um vetor arbitrário, então o produto $$$H (x) \ cdot p$ tem a forma:

Uma função que calcula o produto do Hessian e um vetor arbitrário é passada como o valor do argumento hessp para a função minimize:

 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 

O algoritmo da região de confiança dos gradientes conjugados (Newton)

A baixa condicionalidade da matriz Hessiana e as direções de busca incorretas podem levar ao fato de que o algoritmo dos gradientes conjugados de Newton pode ser ineficiente. Nesses casos, é dada preferência ao método da região de confiança dos gradientes de Newton conjugados.

Exemplo de definição de 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.] 

Um exemplo com a função de produto do Hessian e um vetor arbitrário:

 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 de tipo Krylovsky

Assim como o método trust-ncg, os métodos do tipo Krylovsky são adequados para resolver problemas de larga escala, pois usam apenas produtos vetoriais matriciais. A essência deles é resolver o problema no campo confidencial limitado pelo subespaço de Krylov truncado. Para tarefas incertas, é melhor usar esse método, pois ele usa menos iterações não lineares devido a menos produtos vetoriais de matriz por subtarefa, em comparação com o método trust-ncg. Além disso, a solução para a subtarefa quadrática é mais precisa que o método trust-ncg.
Exemplo de definição de 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.] 

Um exemplo com a função de produto do Hessian e um vetor arbitrário:

 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 Baseado em Confiança

Todos os métodos (Newton-CG, trust-ncg e trust-krylov) são adequados para resolver tarefas de larga escala (com milhares de variáveis). Isto é devido ao fato de que o algoritmo subjacente de gradientes conjugados implica uma determinação aproximada da matriz Hessiana inversa. A solução é iterativa, sem decomposição explícita do hessiano. Como é necessário apenas determinar a função para o produto do Hessian e um vetor arbitrário, esse algoritmo é especialmente bom para trabalhar com matrizes esparsas (diagonal da fita). Isso fornece baixos custos de memória e economia de tempo significativa.

Em problemas de tamanho médio, os custos de armazenamento e fatoração do Hessian não são críticos. Isso significa que uma solução pode ser obtida em menos iterações, resolvendo as subtarefas da área de confiança quase exatamente. Para isso, algumas equações não lineares são resolvidas iterativamente para cada subtarefa quadrática. Essa solução geralmente requer 3 ou 4 decomposições da matriz de Holets Hessian. Como resultado, o método converge em menos iterações e requer menos computação da função objetivo do que outros métodos implementados do domínio de confiança. Este algoritmo implica apenas a determinação da matriz completa do Hessian e não suporta a capacidade de usar a função do produto do Hessian e um vetor arbitrário.

Um exemplo de como minimizar a função 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.]) 

Isso talvez habite. No próximo artigo, tentarei dizer o mais interessante sobre a minimização condicional, a aplicação da minimização na resolução de problemas de aproximação, minimizando a função de uma variável, minimizadores arbitrários e encontrando as raízes do sistema de equações usando o pacote scipy.optimize.

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