SciPy (pronunciato come sai pai) Γ¨ un pacchetto di procedure matematiche applicate, basato sull'estensione Numpy di Python. Con SciPy, una sessione interattiva di Python diventa un ambiente completo per l'elaborazione dei dati e per la prototipazione di sistemi complessi, simile a MATLAB, IDL, Octave, R-Lab e SciLab. Oggi voglio brevemente illustrare come applicare alcuni noti algoritmi di ottimizzazione nel pacchetto scipy.optimize. Puoi sempre ottenere una documentazione piΓΉ dettagliata e aggiornata sulle funzioni utilizzando il comando help() o premendo Shift+Tab.
Introduzione
Per evitare di costringere me stesso e i lettori a cercare e leggere le fonti originali, i riferimenti alle descrizioni dei metodi saranno principalmente su Wikipedia. In genere, queste informazioni sono sufficienti per comprendere i metodi in termini generali e le condizioni della loro applicazione. Per comprendere il fondamento dei metodi matematici, faremo riferimento a pubblicazioni piΓΉ autorevoli, che puoi trovare alla fine di ogni articolo o nel tuo motore di ricerca preferito.
Pertanto, il modulo scipy.optimize include l'implementazione delle seguenti procedure:
- Minimizzazione condizionale e incondizionata di funzioni scalari di piΓΉ variabili (minim) utilizzando diversi algoritmi (simplex di Nelder-Mead, BFGS, gradienti coniugati di Newton, e )
- Ottimizzazione globale (ad esempio: , )
- Minimizzazione dei residui (least_squares) e algoritmi per l'adattamento di curve con MNC non lineare (curve_fit)
- Minimizzazione di funzioni scalari di una variabile (minim_scalar) e ricerca di radici (root_scalar)
- Risolutori multidimensionali per sistemi di equazioni (root) utilizzando diversi algoritmi (ibrido di Powell, oppure metodi su larga scala, come ).
In questo articolo ci concentreremo solo sul primo punto di tutto questo elenco.
Minimizzazione incondizionata di una funzione scalare di piΓΉ variabili
La funzione minim del pacchetto scipy.optimize fornisce un'interfaccia generale per risolvere problemi di minimizzazione condizionale e incondizionata di funzioni scalari di piΓΉ variabili. Per dimostrarne il funzionamento, avremo bisogno di una funzione di piΓΉ variabili adeguata che andremo a minimizzare in vari modi.
Per questo scopo Γ¨ ideale la funzione di Rosenbrock in N variabili, che ha la forma:
Anche se la funzione di Rosenbrock e le sue matrici jacobiane e hessiane (rispettivamente prima e seconda derivata) sono giΓ definite nel pacchetto scipy.optimize, definiamola da noi.
import numpy as np
def rosen(x):
"""La funzione di Rosenbrock"""
return np.sum(100.0*(x[1:]-x[:-1]**2.0)**2.0 + (1-x[:-1])**2.0, axis=0)Per chiarezza, tracciamo in 3D i valori della funzione di Rosenbrock in due variabili.
Codice per la visualizzazione
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
# Configuriamo il grafico 3D
fig = plt.figure(figsize=[15, 10])
ax = fig.gca(projection='3d')
# Impostiamo l'angolo di visualizzazione
ax.view_init(45, 30)
# Creiamo i dati per il grafico
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]))
# Disegniamo la superficie
surf = ax.plot_surface(X, Y, Z, cmap=cm.coolwarm)
plt.show()
Sapendo in anticipo che il minimo Γ¨ 0 per 
, consideriamo esempi su come determinare il valore minimo della funzione di Rosenbrock utilizzando diverse procedure di scipy.optimize.
Metodo del simplesso di Nelder-Mead
Supponiamo di avere un punto iniziale x0 in uno spazio a 5 dimensioni. Troviamo il punto minimo piΓΉ vicino della funzione di Rosenbrock utilizzando l'algoritmo. (l'algoritmo Γ¨ specificato come valore del parametro method):
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)Ottimizzazione terminata con successo.
Valore attuale della funzione: 0.000000
Iterazioni: 339
Valutazioni della funzione: 571
[1. 1. 1. 1. 1.]Il metodo del simplex Γ¨ il modo piΓΉ semplice per minimizzare una funzione esplicitamente definita e piuttosto liscia. Non richiede il calcolo delle derivate della funzione, Γ¨ sufficiente fornire solo i suoi valori. Il metodo di Nelder-Mead Γ¨ una buona scelta per semplici problemi di minimizzazione. Tuttavia, poichΓ© non utilizza le stime del gradiente, potrebbe richiedere piΓΉ tempo per trovare il minimo.
Metodo di Powell
Un altro algoritmo di ottimizzazione che calcola solo i valori delle funzioni Γ¨ . Per usarlo, Γ¨ necessario impostare method = βpowellβ nella funzione 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)Ottimizzazione terminata con successo.
Valore attuale della funzione: 0.000000
Iterazioni: 19
Valutazioni della funzione: 1622
[1. 1. 1. 1. 1.]Algoritmo di Broyden-Fletcher-Goldfarb-Shanno (BFGS)
Per ottenere una convergenza piΓΉ rapida verso la soluzione, la procedura utilizza il gradiente della funzione obiettivo. Il gradiente puΓ² essere fornito come funzione o calcolato tramite differenze di primo ordine. In ogni caso, solitamente il metodo BFGS richiede meno chiamate alla funzione rispetto al metodo simplex.
Calcoliamo analiticamente la derivata della funzione di Rosenbrock:
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 valore del parametro jac della funzione minim, come mostrato di seguito.
res = minimize(rosen, x0, method='BFGS', jac=rosen_der, options={'disp': True})
print(res.x)Ottimizzazione terminata con successo.
Valore attuale della funzione: 0.000000
Iterazioni: 25
Valutazioni della funzione: 30
Valutazioni del gradiente: 30
[1.00000004 1.0000001 1.00000021 1.00000044 1.00000092]Algoritmo dei gradienti coniugati (Newton)
Algoritmo Γ¨ un metodo modificato di Newton.
Il metodo di Newton si basa sull'approssimazione di una funzione in una regione locale mediante un polinomio di secondo grado:
dove 
Γ¨ la matrice delle seconde derivate (matrice di Hess, hessiana).
Se la hessiana Γ¨ definita positiva, si puΓ² trovare il minimo locale di questa funzione ponendo a zero il gradiente della forma quadratica. Ne risulta l'espressione:
L'inverso della hessiana Γ¨ calcolato utilizzando il metodo dei gradienti coniugati. Un esempio di utilizzo di questo metodo per minimizzare la funzione di Rosenbrock Γ¨ riportato di seguito. Per utilizzare il metodo Newton-CG, Γ¨ necessario definire una funzione che calcola la hessiana.
La hessiana della funzione di Rosenbrock in forma analitica Γ¨ uguale a:
dove 
e 
, definisce la matrice 
.
Gli altri elementi non nulli della matrice sono uguali a:
Ad esempio, in uno spazio a cinque dimensioni N = 5, la matrice di Hess per la funzione di Rosenbrock ha una forma a nastro:
Codice che calcola questa hessiana insieme al codice per minimizzare la funzione di Rosenbrock utilizzando il metodo dei gradienti coniugati (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)Ottimizzazione terminata con successo.
Valore attuale della funzione: 0.000000
Iterazioni: 24
Valutazioni della funzione: 33
Valutazioni del gradiente: 56
Valutazioni della Hessiana: 24
[1. 1. 1. 0.99999999 0.99999999]Esempio di definizione della funzione del prodotto della Hessiana e di un vettore arbitrario
Nelle applicazioni pratiche, il calcolo e la memorizzazione dell'intera matrice Hessiana possono richiedere risorse significative in termini di tempo e memoria. Tuttavia, non Γ¨ realmente necessario specificare la matrice Hessiana stessa, poichΓ© per la procedura di minimizzazione Γ¨ sufficiente avere un vettore che rappresenti il prodotto della Hessiana con un altro vettore arbitrario. Pertanto, dal punto di vista computazionale, Γ¨ molto piΓΉ conveniente definire immediatamente una funzione che restituisca il risultato del prodotto della Hessiana con un vettore arbitrario.
Consideriamo la funzione hess, che accetta un vettore di minimizzazione come primo argomento e un vettore arbitrario come secondo argomento (insieme ad altri argomenti della funzione minimizzata). Nel nostro caso, calcolare il prodotto dell' Hessiano della funzione di Rosenbrock con un vettore arbitrario non Γ¨ molto complicato. Se p Γ¨ un vettore arbitrario, il prodotto 
ha la forma:
La funzione che calcola il prodotto dell'Hessiano e di un vettore arbitrario viene passata come valore dell'argomento hessp alla funzione 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})
Ottimizzazione terminata con successo.
Valore attuale della funzione: 0.000000
Iterazioni: 24
Valutazioni della funzione: 33
Valutazioni del gradiente: 56
Valutazioni dell'Hessiano: 66L'algoritmo delle regioni di fiducia (trust region) dei gradienti coniugati (Newton)
Una scarsa definizione della matrice di Hess e direzioni di ricerca errate possono rendere inefficace l'algoritmo dei gradienti coniugati di Newton. In tali casi, si preferisce (trust-region) dei gradienti coniugati di Newton.
Esempio con la definizione della matrice di Hess:
res = minimize(rosen, x0, method='trust-ncg',
jac=rosen_der, hess=rosen_hess,
options={'gtol': 1e-8, 'disp': True})
print(res.x)Ottimizzazione terminata con successo.
Valore attuale della funzione: 0.000000
Iterazioni: 20
Valutazioni della funzione: 21
Valutazioni del gradiente: 20
Valutazioni della Hessiana: 19
[1. 1. 1. 1. 1.]Esempio con la funzione prodotto della Hessiana 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)Ottimizzazione terminata con successo.
Valore attuale della funzione: 0.000000
Iterazioni: 20
Valutazioni della funzione: 21
Valutazioni del gradiente: 20
Valutazioni della Hessiana: 0
[1. 1. 1. 1. 1.]Metodi di tipo Krylov
Come il metodo trust-ncg, i metodi di tipo Krylov sono ben adatti per risolvere problemi su larga scala, poichΓ© utilizzano solo prodotti matrice-vettore. La loro essenza risiede nel risolvere il problema in un'area fidata, limitata da uno spazio ridotto di Krylov. Per problemi indeterminati, Γ¨ meglio usare questo metodo, poichΓ© richiede un numero inferiore di iterazioni non lineari grazie a un minor numero di prodotti matrice-vettore per ogni sottoproblema, rispetto al metodo trust-ncg. Inoltre, la soluzione del sottoproblema quadratico Γ¨ piΓΉ precisa rispetto al metodo trust-ncg.
Esempio con la definizione della matrice di Hess:
res = minimize(rosen, x0, method='trust-krylov',
jac=rosen_der, hess=rosen_hess,
options={'gtol': 1e-8, 'disp': True})
Ottimizzazione terminata con successo.
Valore attuale della funzione: 0.000000
Iterazioni: 19
Valutazioni della funzione: 20
Valutazioni del gradiente: 20
Valutazioni dell' Hessiana: 18
print(res.x)
[1. 1. 1. 1. 1.]
Esempio con la funzione prodotto della Hessiana e un vettore arbitrario:
res = minimize(rosen, x0, method='trust-krylov',
jac=rosen_der, hessp=rosen_hess_p,
options={'gtol': 1e-8, 'disp': True})
Ottimizzazione terminata con successo.
Valore attuale della funzione: 0.000000
Iterazioni: 19
Valutazioni della funzione: 20
Valutazioni del gradiente: 20
Valutazioni dell' Hessiana: 0
print(res.x)
[1. 1. 1. 1. 1.]
Algoritmo per la risoluzione approssimativa nell'area fidata
Tutti i metodi (Newton-CG, trust-ncg e trust-krylov) sono ben adattati per risolvere problemi su larga scala (con migliaia di variabili). Questo Γ¨ dovuto al fatto che l'algoritmo sottostante del gradiente coniugato implica la determinazione approssimativa della matrice inverse di Hess. La soluzione viene trovata in modo iterativo, senza la decomposizione esplicita dell'Hessiano. PoichΓ© Γ¨ necessario solo determinare la funzione per il prodotto dell'Hessiano e di un vettore arbitrario, questo algoritmo Γ¨ particolarmente efficace nel trattare matrici sparse (diagonali a nastro). CiΓ² garantisce un basso consumo di memoria e un notevole risparmio di tempo.
Nei problemi di dimensioni medio-piccole, i costi di memorizzazione e la fattorizzazione dell'ematossidina non sono determinanti. CiΓ² significa che si puΓ² ottenere una soluzione con un numero ridotto di iterazioni, risolvendo le sottoattivitΓ quasi esattamente nell'area di fiducia. Per fare questo, alcune equazioni non lineari vengono risolte in modo iterativo per ogni sottoattivitΓ quadratica. Di solito, questa soluzione richiede 3 o 4 decomposizioni della matrice Hessiana di Cholesky. Di conseguenza, il metodo converge con un minor numero di iterazioni e richiede meno calcoli della funzione obiettivo rispetto ad altri metodi di area di fiducia implementati. Questo algoritmo implica solo la definizione della matrice Hessiana completa e non supporta l'uso della funzione del prodotto dell'ematossidina e di un vettore arbitrario.
Esempio di minimizzazione della funzione di Rosenbrock:
res = minimize(rosen, x0, method='trust-exact',
jac=rosen_der, hess=rosen_hess,
options={'gtol': 1e-8, 'disp': True})
res.xOttimizzazione terminata con successo.
Valore attuale della funzione: 0.000000
Iterazioni: 13
Valutazioni della funzione: 14
Valutazioni del gradiente: 13
Valutazioni dell'Hessiano: 14
array([1., 1., 1., 1., 1.])Su questo, direi basta. Nel prossimo articolo cercherΓ² di raccontare le cose piΓΉ interessanti sulla minimizzazione vincolata, sull'applicazione della minimizzazione nella risoluzione di problemi di approssimazione, sulla minimizzazione di funzioni di una variabile, sui minimizzatori generali e sulla ricerca delle radici di un sistema di equazioni con il pacchetto scipy.optimize.
Fonte:
Fonte: habr.com
