SciPy, optimeerimine

SciPy (hääldatakse kui sai pai) on rakenduslike matemaatiliste protseduuride paket, mis põhineb Python'i Numpy laiendusel. SciPy muudab interaktiivse Python'i seansi samasuguseks täielikuks andmete töötlemise ja keerukate süsteemide prototüüpimise keskkonnaks nagu MATLAB, IDL, Octave, R-Lab ja SciLab. Täna tahan lühidalt rääkida sellest, kuidas rakendada mõningaid tuntud optimeerimisalgoritme scipy.optimize pakendis. Ülevaatlikku ja ajakohast dokumentatsiooni funktsioonide rakendamise kohta saab alati käivitades käsku help() või kasutades Shift+Tab klahvikombinatsiooni.

Sissejuhatus

Kuna tahaksin end ja lugejaid säästa allikate otsimisest ja lugemisest, viidatakse meetodite kirjelduste osas peamiselt vikipeedia lehtedele. Üldiselt piisab sellest infost meetodite üldiseks mõistmiseks ja nende rakendustingimuste tundmiseks. Matemaatiliste meetodite olemuse mõistmiseks suundume linkide kaudu autoriteetsematesse väljaannetesse, mille leiate iga artikli lõpus või oma lemmikotsingumootorist.

Seega sisaldab modul scipy.optimize järgmisi protseduure:

1. Konditsionaalne ja tingimusteta mitme muutuja skalaarsed funktsioonid (minim) erinevate algoritmide (Nelder-Meade simplex, BFGS, Newtoni konjugaatgraafikud) abil. COBYLA ja SLSQP)
2. Globaalne optimeerimine (näiteks: basinhopping, diff_evolution)
3. Jääkide minimaliseerimine MNLK (least_squares) ja mittelineaarsete MNLK kõverate sobitamise algoritmid (curve_fit)
4. Ühe muutujaga skalaarsete funktsioonide minimaliseerimine (minim_scalar) ja juurte otsimine (root_scalar)
5. Mitmemõõtmelised võrrandite lahendajad (root) erinevate algoritmide abil (hübriidne Powell, Levenbergi-Marquardt või suuremassilised meetodid, nagu Newtoni-Krylov).

Käesolevas artiklis vaatleme ainult nimekirja esimest punkti.

Tingimusteta mitme muutuja skalaarsest funktsioonist

Paketi scipy.optimize funktsioon minim pakub üldist liidest mitme muutuja skalaarsed funktsioonid tingimusliku ja tingimusteta minimiseerimiseks. Selle töö demonstreerimiseks vajame sobivat mitme muutuja funktsiooni, mida me tahame erinevalt minimeerida.

Selleks sobib suurepäraselt Rosenbrocki funktsioon N muutuja jaoks, mille kuju on järgmine:

Kuigi Rosenbrocki funktsioon ja tema Jakobi ja Hesse maatriksid (esimene ja teine derivaat vastavalt) on juba määratletud scipy.optimize pakettis, määratlegem see 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 järgi.

Joonistamise kood

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

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

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

# Loome andmed joonise 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 pinda
surf = ax.plot_surface(X, Y, Z, cmap=cm.coolwarm)
plt.show()

Teades eelnevalt, et miinimum on 0, kui , vaatleme näiteid, kuidas leida Rosenbrocki funktsiooni miinimum erinevate scipy.optimize protseduuride abil.

Nelder-Meadi simplex meetod

Olgu meil algpunkt x0 5-mõõtmelises ruumis. Leiame selle lähedal asuva miinimumpunkti Rosenbrocki funktsiooni jaoks, kasutades Nelder-Meadi simplex'i (meetod on määratud parameetri method väärtusena):

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õpetati edukalt.
         Praegune funktsiooni väärtus: 0.000000
         Iteratsioonid: 339
         Funktsiooni hindamised: 571
[1. 1. 1. 1. 1.]

Simplex meetod on kõige lihtsam viis selgelt määratletud ja üsna sujuva funktsiooni miinimumini viimiseks. See ei vaja funktsiooni tuletiste arvutamist, piisab vaid selle väärtuste sisestamisest. Nelder-Meadi meetod on hea valik lihtsate optimeerimisülesannete jaoks. Kuid kuna see ei kasuta gradientide hinnanguid, võib miinimumi leidmiseks kuluda rohkem aega.

Powelli meetod

Teine optimeerimise algoritm, kus arvutatakse ainult funktsioonide väärtusi, on Powelli meetod.Seda kasutamiseks peab meetodi väärtuseks olema seatud '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õpetati edukalt.
         Praegune funktsiooni väärtus: 0.000000
         Iteratsioonid: 19
         Funktsiooni hindamised: 1622
[1. 1. 1. 1. 1.]

Broyden-Fletcher-Goldfarb-Shanno (BFGS) algoritm

Kiirema lahenduse saavutamiseks kasutage protseduuri BFGS kasutab sihtfunktsiooni gradienti. Gradient saab olla määratud funktsioonina või arvutatuna esimese järgu erinevustega. Igal juhul nõuab meetod BFGS tavaliselt vähem funktsioonikutsungeid kui simplex-meetod.

Leiame Rosenbrocki funktsiooni derivaadi analüütiliselt:

See avaldis on kehtiv kõigi muutujate derivaatide jaoks, välja arvatud esimene ja viimane, mis määratakse järgmiselt:

Vaadakem Python'i 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ääratakse minim-funktsiooni parameetri jac väärtusena, nagu on näidatud allpool.

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]

Konjugaatide gradientide algoritm (Newtoni)

Algoritm Newtoni konjugaatide gradientide on muudetud Newtoni meetod.
Newtoni meetod põhineb funktsiooni ligikaudsel esitlemisel kohalikus piirkonnas teise astme polünoomiga:

kus on teise järgu derivaatide maatriks (Hessi maatriks, hessiaan).
Kui hessiaan on positiivselt määratud, saab selle funktsiooni kohalikku miinimumi leida, seades ruutfunktsiooni nullgradienti nulliks. Tulemuseks on avaldis:

Hessi tagasi arvutatakse konjugaatide gradientide meetodi abil. Allpool on näide selle meetodi kasutamisest Rosenbrocki funktsiooni minimeerimiseks. Newtone-CG meetodi kasutamiseks tuleb määrata funktsioon, mis arvutab hessiaani.
Rosenbrocki funktsiooni hessiaan analüütiliselt võrdub:

kus ja , määravad maatriksi .

Ülejäänud mitte null elemendid maatriksis on:

Näiteks viie mõõtmelises ruumis N = 5, Rosenbrocki funktsiooni Hessi maatriks on lintkujuline:

Kood, mis arvutab selle hessiaani koos koodiga Rosenbrocki funktsiooni minimeerimiseks konjugaatide gradientide meetodi (Newtoni) abil:

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õpetati edukalt.
         Praegune funktsiooni väärtus: 0.000000
         Iteratsioonid: 24
         Funktsiooni hindamised: 33
         Gradientide hindamised: 56
         Hessianide hindamised: 24
[1.         1.         1.         0.99999999 0.99999999]

Näide Hessi maatriksi ja juhusliku vektori määramisest

Reaalsetes ülesannetes võib Hessi maatriksi täieliku arvutamise ja salvestamise protsess nõuda märkimisväärseid ressursse, nagu aega ja mälu. Samas pole tegelikult vaja määrata Hessi maatriksit, kuna minimaalsete protseduuride jaoks on vajalik ainult vektor, mis võrdub Hessi maatriksi ja mõne muu juhusliku vektori korrutisega. Seetõttu on arvutuslikult palju mugavam kohe määratleda funktsioon, mis tagastab Hessi maatriksi korrutise juhusliku vektoriga.

Vaatleme funktsiooni hess, mis võtab minimaalse vektori esimesena argumendina ning juhusliku vektori teise argumendina (koos teiste minimiseeritava funktsiooni argumentidega). Meie puhul pole Rosenbrocki funktsiooni Hessi maatriksi korrutise arvutamine juhusliku vektoriga väga keeruline. Kui p — juhuslik vektor, siis korrutis on kujul:

Hessi maatriksi ja juhusliku vektori korrutuse arvutamiseks edastatakse funktsioon argumendina hessp minimiseerimise funktsioonile:

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.
         Praegune funktsiooni väärtus: 0.000000
         Iteratsioonid: 24
         Funktsiooni hindamised: 33
         Gradientide hindamised: 56
         Hessianide hindamised: 66

Usalduspiirkonna (trust region) paarisgradientide (Newtoni) algoritm

Hessi maatriksi halb määramatus ja vale otsingu suunad võivad põhjustada seda, et paarisgradientide Newtoni algoritm võib olla ebatõhus. Sellistes olukordades eelistatakse usalduspiirkonna meetodit (trust-region) paarisgradientide Newtoni meetodit.

Näide Hessi maatriksi määramisest:

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
         Hessiani hindamised: 19
[1. 1. 1. 1. 1.]

Näide hessiani ja juhusliku vektori korrutamise 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
         Hessiani hindamised: 0
[1. 1. 1. 1. 1.]

Krylov'i tüüpi meetodid

Sarnaselt trust-ncg meetodile sobivad Krylov'i tüüpi meetodid hästi suurte probleemide lahendamiseks, kuna need kasutavad ainult maatriksi-vektori korrutamisi. Nende tuum on usaldusväärse piirkonna lahendamine, mida piirab Krylov'i lõigatud alamruum. Määramata probleemide korral on parem kasutada seda meetodit, kuna see vajab vähem mittekonkreetseid iteratsioone madalama maatriksi-vektori korrutamiste arvu tõttu ühe alaprobleemi puhul võrreldes trust-ncg meetodiga. Lisaks leitakse kvadrilise alaprobleemi lahendus täpsemalt kui trust-ncg meetodiga.
Näide Hessi maatriksi määramisest:

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

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

print(res.x)

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

Näide hessiani ja juhusliku vektori korrutamise 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.
         Praegune funktsiooni väärtus: 0.000000
         Iteratsioonid: 19
         Funktsiooni hindamised: 20
         Gradientide hindamised: 20
         Hessiani hindamised: 0

print(res.x)

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

Usaldusväärse piirkonna lähenemisviisi algoritm

Kõik meetodid (Newton-CG, trust-ncg ja trust-krylov) sobivad hästi suurte probleemide lahendamiseks (tuhandete muutujatega). See tuleneb asjaolust, et nende aluseks olev konjugeeritud gradientide algoritm eeldab Hessiani pöördmatrise ligikaudset leidmist. Lahendus leitakse iteratiivselt, ilma Hessiani selges vormis lagundamiseta. Kuna on vajalik määrata ainult funktsioon Hessiani ja juhusliku vektori korrutamiseks, on see algoritm eriti hea haruldaste (vööndi diagonaalsete) maatriskatega töötamiseks. See tagab madalad mälukulud ja olulise ajakokkuhoiu.

Keskmise suurusega ülesannetes salvestamise ja Hessiani faktorisatsiooni kulud ei ole määravad. See tähendab, et lahenduse võib saada vähemate iteratsioonide arvuga, lahendades usalduspiirkonna osauuringud peaaegu täpselt. Selleks lahendatakse mõned mittejoonised iteratiivselt iga ruutülesande jaoks. Selline lahendus nõuab tavaliselt 3-4 Cholesky dekompositsiooni Hessiani maatriksi jaoks. Seetõttu konverdib meetod vähemate iteratsioonidega ja nõuab vähem eesmärgi funktsiooni arvutusi kui teised rakendatud usalduspiirkonna meetodid. See algoritm eeldab ainult täieliku Hessiani maatriksi määratlemist ning ei toeta võimalust kasutada Hessiani funktsiooni ja mis tahes vektori produktsiooni.

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õpetatud eduka 
         Praegune funktsiooni väärtus: 0.000000
         Iteratsioonide arv: 13
         Funktsiooni arvutused: 14
         Gradientide arvutused: 13
         Hessiani arvutused: 14

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

Sellele paneme tõenäoliselt punkti. Järgmises artiklis püüan rääkida kõige huvitavamast tingimuslikust minimeerimisest, minimeerimise rakendamisest approximatsiooni ülesannete lahendamisel, ühe muutuja funktsiooni minimeerimisest, mis tahes minimeerijatest ja võrrandisüsteemide juurte leidmisest scipy.optimize paketi abil.

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

Allikas: habr.com