SciPy (pronunciato come sai pai) è un pacchetto di procedure matematiche applicate basato sull'estensione di Numpy Python. Con SciPy, una sessione interattiva di Python si trasforma in un ambiente completo per l'elaborazione dei dati e il prototipaggio di sistemi complessi, simile a MATLAB, IDL, Octave, R-Lab e SciLab. Oggi voglio raccontarvi brevemente come applicare alcuni noti algoritmi di ottimizzazione nel pacchetto scipy.optimize. Ulteriori informazioni e aggiornamenti sulle funzioni possono essere ottenuti utilizzando il comando help() o premendo Shift+Tab.
Introduzione
Per risparmiare tempo a me stesso e ai lettori nella ricerca e lettura delle fonti primarie, i riferimenti alle descrizioni dei metodi saranno principalmente verso Wikipedia. Di norma, queste informazioni sono sufficienti per comprendere i metodi in modo generale e le condizioni del loro utilizzo. Per comprendere il nocciolo dei metodi matematici, seguiamo i link a pubblicazioni più autorevoli, che possono essere trovate alla fine di ogni articolo o nel motore di ricerca preferito.
Quindi, il modulo scipy.optimize include l'implementazione delle seguenti procedure:
- Minimizzazione condizionata e non condizionata di funzioni scalari di più variabili (minim) utilizzando vari algoritmi (simplex di Nelder-Mead, BFGS, gradienti coniugati di Newton, e )
- Ottimizzazione globale (ad esempio: , )
- Minimizzazione dei residui (least_squares) e algoritmi di fitting di curve non lineari (curve_fit)
- Minimizzazione di funzioni scalari di una variabile (minim_scalar) e ricerca di zeri (root_scalar)
- Risolutori multidimensionali per sistemi di equazioni (root) utilizzando vari algoritmi (ibrido di Powell, o metodi su larga scala come ).
In questo articolo considereremo solo il primo punto di questo elenco.
Minimizzazione non condizionata di una funzione scalare di più variabili
La funzione minim del pacchetto scipy.optimize offre un'interfaccia generale per risolvere problemi di minimizzazione condizionata e non condizionata di funzioni scalari di più variabili. Per dimostrare il suo funzionamento, avremo bisogno di una funzione adeguata di più variabili che minimizzeremo in modo diverso.
A tal fine, la funzione di Rosenbrock di N variabili, che ha la seguente forma:
Nonostante la funzione di Rosenbrock e le sue matrici di Jacobi e Hess (rispettivamente prima e seconda derivata) siano già definite nel pacchetto scipy.optimize, definiamola autonomamente.
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 base a due variabili.
Codice per il tracciamento
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 vista
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 a 
, consideriamo esempi di come determinare il valore minimo della funzione di Rosenbrock utilizzando diverse procedure di scipy.optimize.
Il metodo del simplesso di Nelder-Mead
Supponiamo di avere un punto iniziale x0 nello spazio a 5 dimensioni. Troviamo il punto di 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 simplesso è il modo più semplice per minimizzare una funzione ben definita e abbastanza liscia. Non richiede il calcolo delle derivate della funzione, basta fornire solo i suoi valori. Il metodo di Nelder-Mead è una buona scelta per compiti di minimizzazione semplici. Tuttavia, poiché non utilizza stime del gradiente, potrebbero essere necessari più tempo per trovare il minimo.
Il metodo di Powell
Un altro algoritmo di ottimizzazione, in cui vengono calcolati solo i valori delle funzioni, è . Per utilizzarlo, è 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.]L'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 espresso come una funzione o calcolato tramite differenze di primo ordine. In ogni caso, di solito il metodo BFGS richiede meno chiamate rispetto al metodo del simplesso.
Calcoliamo 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 per calcolare il 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 della funzione in un'area locale tramite un polinomio di secondo grado:
dove 
è la matrice delle derivate seconde (matrice di Hessian, hessiano).
Se l'hessiano è definito positivo, il minimo locale di questa funzione può essere trovato ponendo a zero il gradiente della forma quadratica. Di conseguenza, si ottiene l'espressione:
L'inverso dell'hessiano viene 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 specificare una funzione che calcoli l'hessiano.
L'hessiano della funzione di Rosenbrock in forma analitica è:
dove 
e 
, definiscono la matrice 
.
Gli altri elementi non nulli della matrice sono:
Ad esempio, in uno spazio a cinque dimensioni N = 5, la matrice di Hessian per la funzione di Rosenbrock ha una forma a nastro:
Codice che calcola questo hessiano 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 dell'Hessiana: 24
[1. 1. 1. 0.99999999 0.99999999]Esempio di definizione della funzione del prodotto dell'Hessiana e di un vettore arbitrario
Nei problemi reali, il calcolo e la memorizzazione dell'intera matrice di Hessiana possono richiedere risorse significative in termini di tempo e memoria. Tuttavia, non è necessario definire l'intera matrice di Hessiana, poiché per la procedura di minimizzazione è sufficiente avere un vettore che rappresenti il prodotto dell'Hessiana con un altro vettore arbitrario. Pertanto, dal punto di vista computazionale, è molto più vantaggioso definire direttamente una funzione che restituisce il risultato del prodotto dell'Hessiana con un vettore arbitrario.
Consideriamo la funzione 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 dell'Hessiana della funzione di Rosenbrock con un vettore arbitrario non è molto complicato. Se p — un vettore arbitrario, allora il prodotto 
ha la seguente forma:
La funzione che calcola il prodotto dell'Hessiana e di un vettore arbitrario viene passata come valore dell'argomento hessp della 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'Hessiana: 66Algoritmo del trust region dei gradienti coniugati di Newton
Una cattiva condizione della matrice di Hessiana 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 di definizione della matrice di Hessiana:
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 dell'hessiano: 19
[1. 1. 1. 1. 1.]Esempio con la funzione del prodotto dell'hessiano 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 dell'hessiano: 0
[1. 1. 1. 1. 1.]Metodi di tipo Krylov
Analogamente al metodo trust-ncg, i metodi di tipo Krylov sono ben adattati per risolvere problemi su larga scala, poiché utilizzano solo prodotti matriciali-vettoriali. La loro essenza consiste nel risolvere il problema in un'area fidata, limitata da uno spazio ridotto di Krylov. Per i problemi non definiti è meglio utilizzare questo metodo, poiché richiede un numero minore di iterazioni non lineari grazie a un numero inferiore di prodotti matriciali-vettoriali per ogni sottoproblema, rispetto al metodo trust-ncg. Inoltre, la soluzione del sottoproblema quadratico viene trovata in modo più preciso rispetto al metodo trust-ncg.
Esempio di definizione della matrice di Hessiana:
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'hessiano: 18
print(res.x)
[1. 1. 1. 1. 1.]
Esempio con la funzione del prodotto dell'hessiano 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'hessiano: 0
print(res.x)
[1. 1. 1. 1. 1.]
Algoritmo di soluzione approssimativa in un'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 dei gradienti coniugati implica la determinazione approssimativa della matrice inversa dell'hessiano. La soluzione viene trovata iterativamente, senza esplicitare la scomposizione dell'hessiano. Poiché è necessario solo definire una funzione per il prodotto dell'hessiano e un vettore arbitrario, questo algoritmo è particolarmente efficace nel lavorare con matrici sparse (diagonali a bande). Questo garantisce costi di memoria ridotti e un notevole risparmio di tempo.
Nei problemi di dimensioni medie, i costi di memorizzazione e la fattorizzazione dell' Hessiana non sono decisivi. Ciò significa che è possibile ottenere una soluzione in un numero minore di iterazioni, risolvendo quasi esattamente le sottoattività dell'area di fiducia. A tal fine, alcune equazioni non lineari vengono risolte iterativamente per ogni sottoattività quadratica. Questa soluzione richiede di solito 3 o 4 fattorizzazioni di Cholesky della matrice Hessiana. Di conseguenza, il metodo converge in un minor numero di iterazioni e richiede meno calcoli della funzione obiettivo rispetto ad altri metodi di area di fiducia implementati. Questo algoritmo prevede solo la determinazione della matrice Hessiana completa e non supporta la possibilità di utilizzare la funzione di prodotto dell'Hessiano 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.])Qui ci fermiamo. Nel prossimo articolo cercherò di raccontare le cose più interessanti sulla minimizzazione condizionale, sull'applicazione della minimizzazione nella soluzione di problemi di approssimazione, sulla minimizzazione di funzioni di una variabile, sui minimizzatori arbitrari e sulla ricerca delle radici di un sistema di equazioni utilizzando il pacchetto scipy.optimize.
Fonte:
Fonte: habr.com
