SciPy, optimeerimine

SciPy (hääldatakse nagu sai pai) on rakenduslike matemaatiliste protseduuride paket, mis põhineb Python'i Numpy laiendusel. SciPy muudab interaktiivse Python seansi täielikuks andmete töötlemise ja keerukate süsteemide prototüüpimise keskkonnaks, võrreldav MATLAB'i, IDL'i, Octave'i, R-Labi ja SciLabiga. Täna tahan ma lühidalt rääkida, kuidas rakendada mõningaid tuntud optimeerimisalgoritme scipy.optimize paketis. Podrobnejat ja ajakohast teavet funktsioonide kasutamise kohta võib alati saada käsu help() või Shift+Tab abil.

Sissejuhatus

Kuna soovin ennast ja lugejat säästa esmaste allikate otsimisest ja lugemisest, on linkide viidatud meetodite kirjeldustele peamiselt Wikipedia. Üldiselt on see teave piisav, et mõista meetodeid laiemalt ja nende rakendustingimusi. Matemaatiliste meetodite südamiku mõistmiseks järgime linke autoriteetsematele väljaannetele, mida saab leida iga artikli lõpus või oma lemmikotsingumootorist.

Seega sisaldab moodul scipy.optimize järgmiste protseduuride teostust:

1. Mitme muutuja skalarfunktsioonide tingimuslik ja tingimusteta minimeerimine (minim) erinevate algoritmidega (Nelder-Mead simplex, BFGS, Newtoni konjugaatgradendid, COBYLA ja SLSQP)
2. Globaalne optimeerimine (näiteks: basinhopping, diff_evolution)
3. Jääkide minimeerimine MNLK (least_squares) ja mittelineaarse MNLK (curve_fit) kõverate sobitamise algoritmid
4. Skalarse funktsiooni minimeerimine ühe muutuja korral (minim_scalar) ja juurte leidmine (root_scalar)
5. Mitme mõõtmelistele võrrandite süsteemide lahendamine (root) erinevate algoritmide abil (hübriidne Powell, Levenbergi-Marquardt või suurte meetodite, nagu Newtoni-Krylov).

Selles artiklis käsitleme vaid loendi esimest punkti.

Tingimusteta skalarse funktsiooni minimeerimine mitme muutuja korral

Scipy.optimize paketi funktsioon minim pakub üldist liidest mitme muutuja skalarfunktsioonide tingimusliku ja tingimusteta minimeerimise probleemide lahendamiseks. Selle töö demonstreerimiseks vajame sobivat mitme muutuja funktsiooni, mida me eri viisil minimeerime.

Selleks sobib suurepäraselt Rosenbrocki funktsioon, mis on määratud N muutuja jaoks.

Kuigi Rosenbrocki funktsioon ja selle Jakobian ja Hesseni maatriksid (esimene ja teine järeldus vastavalt) on juba määratletud paketis scipy.optimize, määratleme selle ise.

import numpy as np

def rosen(x):
    """Rosenbrocki funktsioon"""
    return np.sum(100.0*(x[1:]-x[:-1]**2.0)**2.0 + (1-x[:-1])**2.0, axis=0)

Selguse huvides joonistame 3D-s Rosenbrocki funktsiooni väärtused kahe muutuja osas.

Kood joonistamiseks

from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter

# Seadistame 3D diagrammi
fig = plt.figure(figsize=[15, 10])
ax = fig.gca(projection='3d')

# Määrame vaatenurga
ax.view_init(45, 30)

# Loome andmed diagrammi jaoks
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]))

# Joonistame pinna
surf = ax.plot_surface(X, Y, Z, cmap=cm.coolwarm)
plt.show()

Teades ette, et miinimum on 0, kui , vaatleme Rosenbrocki funktsiooni miinimumväärtuse kindlaksmääramise näiteid erinevate scipy.optimize protseduuride abil.

Nelder-Meadi simplex meetod

Oletame, et algne punkt x0 asub 5-mõõtmelises ruumis. Otsime läheima miinimumpunkti Rosenbrocki funktsiooni jaoks, kasutades algoritmi. Nelder-Meadi simplex (algoritm on määratud parameetri method väärtuseks):

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)

Optimeerimine lõppes edukalt.
         Praegune funktsiooni väärtus: 0.000000
         Iteratsioonid: 339
         Funktsiooni hindamisi: 571
[1. 1. 1. 1. 1.]

Simpelksimeetod on kõige lihtsam viis minimaalsete selgelt määratletud ja üsna sujuvate funktsioonide saavutamiseks. See ei nõua funktsiooni tuletiste arvutamist, piisab ainult selle väärtuste määramisest. Nelder-Meadi meetod on hea valik lihtsate optimeerimisülesannete jaoks. Kuna see ei kasuta gradientide hinnanguid, võib miinimumi leidmiseks kuluda rohkem aega.

Pauli meetod

Teine optimeerimisalgoritm, kus arvutatakse ainult funktsiooni väärtused, on Pauli meetod. Selle kasutamiseks tuleb seada method = 'powell' funktsioonis 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)

Optimeerimine lõppes edukalt.
         Praegune funktsiooni väärtus: 0.000000
         Iteratsioonid: 19
         Funktsiooni hindamisi: 1622
[1. 1. 1. 1. 1.]

Broyden-Fletcher-Goldfarb-Shanno (BFGS) algoritm

Kiirema lahenduse saavutamiseks protseduur BFGS kasutab sihtfunktsiooni gradienti. Gradient võib olla määratletud funktsiooni abil või arvutatud esmakordsete vahede abil. Igal juhul nõuab BFGS meetod tavaliselt vähem funktsioonikutsungeid kui simplex-meetod.

Leidke Rosenbrocki funktsiooni tuletis analüütiliselt:

See avaldis kehtib kõigi muutujate tuletiste kohta, välja arvatud esimene ja viimane, mis määratletakse järgmiselt:

Vaatame Pythonis funktsiooni, mis arvutab selle gradienti:

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

Gradienti arvutamise funktsioon määratletakse minim funktsiooni parameetri jac väärtusena, nagu allpool näidatud.

res = minimize(rosen, x0, method='BFGS', jac=rosen_der, options={'disp': True})
print(res.x)

Optimeerimine lõpetati edukalt.
         Praegune funktsiooni väärtus: 0.000000
         Iteratsioonid: 25
         Funktsiooni hindamised: 30
         Gradientide hindamised: 30
[1.00000004 1.0000001  1.00000021 1.00000044 1.00000092]

Kaheksanda astme gradientide (Newtoni) algoritm

Algoritm Newtoni konjugaatide gradientide on muudetud Newtoni meetod.
Newtoni meetod põhineb funktsiooni lähialas teise astme polünoomi abil lähendamisel:

kus on teise tuletise maatriks (Hesse maatriks, hessián).
Kui hessián on positiivselt määratud, saab selle funktsiooni lokaalse miinimumi leida, seades ruutfunktsiooni nulligradiendi nulliks. Tulemusena saadakse väljend:

Tagumine hessián arvutatakse konjugeeritud gradientide meetodi abil. Allpool on toodud näide selle meetodi kasutamisest Rosenbrocki funktsiooni minimeerimiseks. Newton-CG meetodi kasutamiseks peab olema määratud funktsioon, mis arvutab hessiáni.
Rosenbrocki funktsiooni hessián analüütiliselt on:

kus ja , määravad maatriksi .

Ülejäänud mitte nullid maatriksi elemendid on:

Näiteks viie mõõtmelises ruumis N = 5, Hesse maatriks Rosenbrocki funktsiooni jaoks on ribakujuline:

Kood, mis arvutab selle hessiáni koos koodiga Rosenbrocki funktsiooni minimeerimiseks konjugeeritud gradientide meetodi abil (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)

Optimeerimine lõpetatud edukalt.
         Praegune funktsiooni väärtus: 0.000000
         Iteratsioonide arv: 24
         Funktsiooni hindamised: 33
         Gradientide hindamised: 56
         Hessiani hindamised: 24
[1.         1.         1.         0.99999999 0.99999999]

Näide hessiaani ja juhusliku vektori funktsiooni määratlemisest

Tegelike probleemide korral võib kogu Hessiani maatriksi arvutamine ja salvestamine nõuda märkimisväärseid aja- ja mäluresursse. Samas pole Hessiani maatriksi määratlemine tegelikult vajalik, kuna minimaalsuse protseduuriks on vajalik ainult vektor, mis on võrdsed Hessiani tulemusega, mis on korrutatud mingi teise juhusliku vektoriga. Seega on arvutuslikult palju eelistatavam kohe määratleda funktsioon, mis tagastab Hessiani ja juhusliku vektori korrutise tulemuse.

Analüüsime funktsiooni hess, mis võtab optimeerimise vektori esimeseks argumendiks ja suvalise vektori teiseks argumendiks (koos teiste minimaalsete funktsioonide argumentidega). Meie puhul ei ole Rosenbrocki funktsiooni hessi vektori suvalise vektoriga korrutamine väga keeruline. Kui p — suvaline vektor, siis korrutamine on järgmine:

Funktsioon, mis arvutab hessi ja suvalise vektori korrutise, edastatakse argumendi hessp väärtusena funktsioonile 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})

Optimeerimine lõpetati edukalt.
         Hetke funktsiooniväärtus: 0.000000
         Iteratsioonid: 24
         Funktsiooniväärdustuste arv: 33
         Gradientide hindamised: 56
         Hessi hindamised: 66

Usalduspiiri (trust region) lineaarsete gradientide (Newtoni) algoritm

Halva Hesse'i maatriksi halb tingimus ja vale otsingu suunad võivad põhjustada, et Newtoni seotud gradientide algoritm osutub ebaefektiivseks. Sellistel juhtudel eelistatakse usalduspiiri meetodit (trust-region) Newtoni seotud gradientide.

Hesse'i maatriksi määratlemise näide:

res = minimize(rosen, x0, method='trust-ncg',
               jac=rosen_der, hess=rosen_hess,
               options={'gtol': 1e-8, 'disp': True})
print(res.x)

Optimeerimine lõpetati edukalt.
         Praegune funktsiooni väärtus: 0.000000
         Iteratsioonid: 20
         Funktsiooni hindamised: 21
         Gradientide hindamised: 20
         Hesse'i hindamised: 19
[1. 1. 1. 1. 1.]

Näide Hesse'i maatriksi ja suvalise vektori korrutise funktsioonist:

res = minimize(rosen, x0, method='trust-ncg', 
                jac=rosen_der, hessp=rosen_hess_p, 
                options={'gtol': 1e-8, 'disp': True})
print(res.x)

Optimeerimine lõpetati edukalt.
         Praegune funktsiooni väärtus: 0.000000
         Iteratsioonid: 20
         Funktsiooni hindamised: 21
         Gradientide hindamised: 20
         Hesse'i hindamised: 0
[1. 1. 1. 1. 1.]

Krylov'i tüüpmeetodid

Trust-ncg meetodile sarnased Krylov tüüpi meetodid sobivad suuremate ülesannete lahendamiseks hästi, kuna need kasutavad ainult maatriks-vektorite korrutusi. Nende olemus seisneb ülesande lahendamises usaldusväärses valdkonnas, mis on piiratud kärbitud Krylov'i alamruumiga. Määramatute probleemide puhul on parem kasutada seda meetodit, kuna see nõuab vähem mittelineaarseid iteratsioone, tänu väiksemale maatriks-vektorite korrutuste arvule iga alamülesande jaoks võrreldes trust-ncg meetodiga. Lisaks on ruutude alamülesande lahendus täpsem kui trust-ncg meetodil.
Hesse'i maatriksi määratlemise näide:

res = minimize(rosen, x0, method='trust-krylov',
               jac=rosen_der, hess=rosen_hess,
               options={'gtol': 1e-8, 'disp': True})

Optimeerimine lõpetati edukalt.
         Hetke funktsiooni väärtus: 0.000000
         Iteratsioonid: 19
         Funktsiooni hindamised: 20
         Gradientide hindamised: 20
         Hessiaani hindamised: 18

print(res.x)

    [1. 1. 1. 1. 1.]

Näide Hesse'i maatriksi ja suvalise vektori korrutise funktsioonist:

res = minimize(rosen, x0, method='trust-krylov',
               jac=rosen_der, hessp=rosen_hess_p,
               options={'gtol': 1e-8, 'disp': True})

Optimeerimine lõpetati edukalt.
         Hetke funktsiooni väärtus: 0.000000
         Iteratsioonid: 19
         Funktsiooni hindamised: 20
         Gradientide hindamised: 20
         Hessiaani hindamised: 0

print(res.x)

    [1. 1. 1. 1. 1.]

Usaldusväärses valdkonnas ligikaudse lahendamise algoritm

Kõik meetodid (Newton-CG, trust-ncg ja trust-krylov) sobivad suurepäraselt ulatuslike probleemide lahendamiseks (tuhandete muutujatega). See on seotud sellega, et aluseks olev gradientide paaritusalgoritm eeldab Hessiani pöördmatriisi ligikaudset leidmist. Lahendust otsitakse iteratiivselt, ilma Hessiani avaldamiseta. Kuna on vajalik määrata ainult funktsioon, mis sisaldab Hessiani ja suvalise vektori korrutamist, on see algoritm eriti efektiivne hõredate (lint-diagonaalsete) maatriksitega töötamisel. See tagab madalad mälukulud ja märkimisväärse ajakasutuse kokkuhoiu.

Keskmise suurusega ülesannetes ei ole Hessiandi salvestamise ja faktoriseerimise kulud kriitilise tähtsusega. See tähendab, et lahendus saab kätte vähemate iteratsioonidega, lahendades usalduspiirkonna alateemasid peaaegu täpselt. Selleks lahendatakse mõned mitte-lineaarsed võrrandid iteratiivselt iga ruutfunktsiooni alateema jaoks. Selline lahendus nõuab tavaliselt 3 või 4 Hessiandi Cholecki lagundamist. Tulemusena konvergib meetod vähemate iteratsioonidega ja vajab vähem eesmärgifunktsiooni arvutusi kui teised rakendatud usalduspiirkonna meetodid. See algoritm eeldab vaid Hessiandi täieliku maatriksi määratlemist ja ei toeta Hessiandi ja suvalise vektori tootefunktsiooni kasutamise võimalust.

Rosenbrocki funktsiooni minimeerimise näide:

res = minimize(rosen, x0, method='trust-exact',
               jac=rosen_der, hess=rosen_hess,
               options={'gtol': 1e-8, 'disp': True})
res.x

Optimeerimine lõpetati edukalt.
         Praegune funktsiooni väärtus: 0.000000
         Iteratsioonide arv: 13
         Funktsiooni hindamised: 14
         Gradientide hindamised: 13
         Hessiandi hindamised: 14

array([1., 1., 1., 1., 1.])

Käesolevaga lõpetame. Järgmistes artiklites püüan rääkida huvitavamatest teemadest, sealhulgas tingimuslikust minimeerimisest, minimeerimise rakendamisest approximatsiooni ülesannete lahendamisel, ühe muutuja funktsiooni minimeerimisest, erinevatest minimeerijatest ning võrrandisüsteemi juurte otsimisest paketiga scipy.optimize.

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

Allikas: habr.com