SciPy (pronunciato sai pie) è un pacchetto di applicazioni matematiche basato sull'estensione Numpy Python. Con SciPy, la tua sessione interattiva di Python diventa lo stesso ambiente completo di data science e prototipazione di sistemi complessi di MATLAB, IDL, Octave, R-Lab e SciLab. Oggi voglio parlare brevemente di come utilizzare alcuni noti algoritmi di ottimizzazione presenti nel pacchetto scipy.optimize. Un aiuto più dettagliato e aggiornato sull'uso delle funzioni può sempre essere ottenuto utilizzando il comando help() o utilizzando Shift+Tab.
Introduzione
Per evitare a te stesso e ai lettori di cercare e leggere fonti primarie, i collegamenti alle descrizioni dei metodi si troveranno principalmente su Wikipedia. Di norma, queste informazioni sono sufficienti per comprendere i metodi in termini generali e le condizioni per la loro applicazione. Per comprendere l'essenza dei metodi matematici, segui i collegamenti alle pubblicazioni più autorevoli, che puoi trovare alla fine di ogni articolo o nel tuo motore di ricerca preferito.
Quindi, il modulo scipy.optimize prevede l'implementazione delle seguenti procedure:
- Minimizzazione condizionale e incondizionata di funzioni scalari di più variabili (minima) utilizzando vari algoritmi (simplex di Nelder-Mead, BFGS, gradienti coniugati di Newton, и )
- Ottimizzazione globale (ad esempio: , )
- Minimizzazione dei residui (least_squares) e algoritmi di adattamento della curva che utilizzano i minimi quadrati non lineari (curve_fit)
- Minimizzazione delle funzioni scalari di una variabile (minim_scalar) e ricerca delle radici (root_scalar)
- Risolutori multidimensionali di sistemi di equazioni (radice) che utilizzano vari algoritmi (ibrido Powell, o metodi su larga scala come ).
In questo articolo considereremo solo il primo elemento dell'intero elenco.
Minimizzazione incondizionata di una funzione scalare di più variabili
La funzione minimi del pacchetto scipy.optimize fornisce un'interfaccia generale per risolvere problemi di minimizzazione condizionale e incondizionata di funzioni scalari di più variabili. Per dimostrare come funziona, avremo bisogno di una funzione adeguata di più variabili, che minimizzeremo in diversi modi.
Per questi scopi è perfetta la funzione di Rosenbrock di N variabili, che ha la forma:
Nonostante il fatto che la funzione Rosenbrock e le sue matrici Jacobi e Hessiane (rispettivamente la prima e la seconda derivata) siano già definite nel pacchetto scipy.optimize, la definiremo noi stessi.
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)Per chiarezza disegniamo in 3D i valori della funzione di Rosenbrock di due variabili.
Codice disegno
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()
Sapendo in anticipo che il minimo è 0 a , diamo un'occhiata ad esempi su come determinare il valore minimo della funzione Rosenbrock utilizzando varie procedure scipy.optimize.
Metodo del simplesso di Nelder-Mead
Sia presente un punto iniziale x0 nello spazio a 5 dimensioni. Troviamo il punto minimo della funzione Rosenbrock più vicino ad essa utilizzando l'algoritmo (l'algoritmo è specificato come valore del parametro del metodo):
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.]Il metodo del simplesso è il modo più semplice per minimizzare una funzione definita esplicitamente e abbastanza fluida. Non è necessario calcolare le derivate di una funzione; è sufficiente specificarne solo i valori. Il metodo Nelder-Mead è una buona scelta per semplici problemi di minimizzazione. Tuttavia, poiché non utilizza stime del gradiente, potrebbe essere necessario più tempo per trovare il minimo.
Metodo Powell
Un altro algoritmo di ottimizzazione in cui vengono calcolati solo i valori delle funzioni è . Per usarlo, è necessario impostare metodo = 'powell' nella funzione minimi.
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 di Broyden-Fletcher-Goldfarb-Shanno (BFGS).
Per ottenere una convergenza più rapida verso una soluzione, la procedura utilizza il gradiente della funzione obiettivo. Il gradiente può essere specificato come funzione o calcolato utilizzando differenze del primo ordine. In ogni caso, il metodo BFGS richiede in genere meno chiamate di funzione rispetto al metodo simplex.
Troviamo la derivata della funzione di Rosenbrock in forma analitica:
Questa espressione è valida per le derivate di tutte le variabili tranne la prima e l'ultima, che sono definite come:
Diamo un'occhiata alla funzione Python che calcola questo 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 funzione di calcolo del gradiente è specificata come il valore del parametro jac della funzione minimi, come mostrato di seguito.
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 del gradiente coniugato (Newton)
Algoritmo è un metodo di Newton modificato.
Il metodo di Newton si basa sull'approssimazione di una funzione in un'area locale mediante un polinomio di secondo grado:
dove è la matrice delle derivate seconde (matrice Hessiana, Hessiana).
Se l'Assia è definita positiva, il minimo locale di questa funzione può essere trovato eguagliando a zero il gradiente zero della forma quadratica. Il risultato sarà l'espressione:
L'Hessiano inverso viene calcolato utilizzando il metodo del gradiente coniugato. Di seguito è riportato un esempio dell'utilizzo di questo metodo per ridurre al minimo la funzione Rosenbrock. Per utilizzare il metodo Newton-CG, è necessario specificare una funzione che calcoli l'Hessiana.
La funzione dell'Assia della Rosenbrock in forma analitica è uguale a:
dove и , definire la matrice .
I restanti elementi diversi da zero della matrice sono uguali a:
Ad esempio, in uno spazio pentadimensionale N = 5, la matrice dell'Assia per la funzione Rosenbrock ha la forma di una banda:
Codice che calcola questa iuta insieme al codice per minimizzare la funzione Rosenbrock utilizzando il metodo del gradiente coniugato (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 esempio con la definizione della funzione prodotto dell'Hessiano e un vettore arbitrario
Nei problemi del mondo reale, il calcolo e l'archiviazione dell'intera matrice dell'Assia può richiedere notevoli risorse di tempo e di memoria. In questo caso, in realtà non è necessario specificare la matrice dell'Assia stessa, perché la procedura di minimizzazione richiede solo un vettore uguale al prodotto dell'Hessiano con un altro vettore arbitrario. Pertanto, da un punto di vista computazionale, è di gran lunga preferibile definire immediatamente una funzione che restituisca il risultato del prodotto dell'Hessiano con un vettore arbitrario.
Consideriamo la funzione di Hess, che accetta il vettore di minimizzazione come primo argomento e un vettore arbitrario come secondo argomento (insieme ad altri argomenti della funzione da minimizzare). Nel nostro caso, calcolare il prodotto della funzione dell'Assia della Rosenbrock con un vettore arbitrario non è molto difficile. Se p è un vettore arbitrario, quindi il prodotto sembra:
La funzione che calcola il prodotto dell'Hessiano e di un vettore arbitrario viene passata come valore dell'argomento hessp alla funzione di minimizzazione:
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 della regione di fiducia del gradiente coniugato (Newton)
Uno scarso condizionamento della matrice assiana e direzioni di ricerca errate possono rendere inefficace l'algoritmo del gradiente coniugato di Newton. In questi casi viene data preferenza a (regione di fiducia) coniuga i gradienti di Newton.
Esempio con la definizione della matrice 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.]Esempio con la funzione prodotto dell'Assia e un vettore 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.]Metodi di tipo Krylov
Come il metodo trust-ncg, i metodi di tipo Krylov sono adatti per risolvere problemi su larga scala perché utilizzano solo prodotti di matrice vettoriale. La loro essenza è risolvere un problema in una regione di confidenza limitata da un sottospazio di Krylov troncato. Per problemi incerti, è meglio utilizzare questo metodo, poiché utilizza un numero minore di iterazioni non lineari a causa del minor numero di prodotti matrice-vettore per sottoproblema, rispetto al metodo trust-ncg. Inoltre, la soluzione al sottoproblema quadratico viene trovata in modo più accurato rispetto all'utilizzo del metodo trust-ncg.
Esempio con la definizione della matrice 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.]
Esempio con la funzione prodotto dell'Assia e un vettore 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 per la soluzione approssimata nella regione di confidenza
Tutti i metodi (Newton-CG, trust-ncg e trust-krylov) sono adatti per risolvere problemi su larga scala (con migliaia di variabili). Ciò è dovuto al fatto che l'algoritmo del gradiente coniugato sottostante implica una determinazione approssimativa della matrice Hessiana inversa. La soluzione viene trovata iterativamente, senza espansione esplicita dell'Assia. Poiché è necessario definire solo una funzione per il prodotto di un vettore dell'Assia e di un vettore arbitrario, questo algoritmo è particolarmente utile per lavorare con matrici sparse (diagonale di banda). Ciò garantisce bassi costi di memoria e un notevole risparmio di tempo.
Per problemi di medie dimensioni, il costo di stoccaggio e factoring dell'Assia non è critico. Ciò significa che è possibile ottenere una soluzione in meno iterazioni, risolvendo quasi esattamente i sottoproblemi della trust regione. Per fare ciò, alcune equazioni non lineari vengono risolte iterativamente per ciascun sottoproblema quadratico. Tale soluzione richiede solitamente 3 o 4 decomposizioni di Cholesky della matrice dell'Assia. Di conseguenza, il metodo converge in meno iterazioni e richiede meno calcoli della funzione obiettivo rispetto ad altri metodi implementati per la regione di confidenza. Questo algoritmo implica solo la determinazione della matrice dell'Assia completa e non supporta la possibilità di utilizzare la funzione prodotto dell'Assia e un vettore arbitrario.
Esempio con minimizzazione della funzione 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.])Probabilmente ci fermeremo qui. Nel prossimo articolo cercherò di raccontare le cose più interessanti sulla minimizzazione condizionale, sull'applicazione della minimizzazione nella risoluzione di problemi di approssimazione, sulla minimizzazione di una funzione di una variabile, sui minimizzatori arbitrari e sulla ricerca delle radici di un sistema di equazioni utilizzando scipy.optimize pacchetto.
Fonte:
Fonte: habr.com