Jak zrychlit Powell Method (hledání globálního optima funkce)?

[GloGlo]

Rozumíte někdo mechanismu toho, jak funguje Powell Method - hledání minima/maxima funkce bez použití gradientů? Používám scipy.optimize.minimize ve snaze najít globální maximum mé funkce MyFunction s mnoha proměnnými (500+). (MyFunction je bohužel step function a její derivace je pro jakoukoli hodnotu zero gradient, takže nejde použít metody využívající derivace).

Když to hledání maxima spustím s nějakými náhodnými výchozími hodnotami, třeba to začne na úrovni, kdy moje funkce dá hodnotu -7000000 - tak to většinou udělá celkem rychle progres a najde hodnotu funkce o dost vyšší než ta výchozí hodnota - třeba skončí na -1000000 nebo +500000. Měl jsem zato, že to je tedy maximální hodnota funkce, kterou to dokáže najít a tím jsem měl věc za uzavřenou.

Jenže sem si všimnul, že když si uložím ty koeficienty z toho jakoby nejúspěšnějšího pokusu a spustím hledání znovu s těmito koeficienty jako výchozími hodnotami - tak to v některých případech vrátí výsledek, který je malinko lepší.

Takže jsem si udělal loop, který vypadá takhle:

Kód: [Vybrat]

start = np.random.sample(532)
res = minimize(MyFunction, start, method='Powell')

while True:
   res = minimize(MyFunction, start, method='Powell')
   start = res.xTakhle se to postupně propracovává k lepšímu a lepšímu výsledku, příklad:

Kód: [Vybrat]

   
      2 1210270.0
      1 1210780.0
      1 1213950.0
    141 1214800.0
      1 1216230.0
      7 1219820.0
      2 1220800.0
      1 1222210.0
      7 1222500.0
      1 1228050.0
      4 1229650.0
     43 1231510.0
     95 1232520.0
     22 1234020.0
     38 1234320.0
      3 1243910.0
     60 1245490.0
     20 1245700.0
     14 1245840.0
      1 1246830.0
      1 1247830.0
    117 1247940.0
     15 1250220.0
     38 1255260.0
     11 1256940.0
     44 1257950.0
     11 1258060.0
    362 1265900.0
     94 1267220.0
      8 1269460.0
    428 1270720.0
    604 1271960.0
     94 1272000.0
      5 1272840.0
    141 1275320.0
      5 1275520.0
    200 1276930.0Ten druhý sloupeček je aktuální hodnota globálního maxima MyFunction, je vidět, že to hledání funguje. Problém je v tom prvním čísle - to je kolikrát skript spustil minimize() s koeficienty odpovídajícími jedné hodnotě funkce než to našlo vyšší.

Problém je v tom, že jedno to spuštění funkce minimize udělá cca 3-5000 evaluací MyFunction a to na AMD EPYC 7B13 (CPU MHz: 3049.996, BogoMIPS: 6099.99) trvá cca 90-120 sekund. Takže když to bylo aktuálně na hodnotě 1271960, tak trvalo asi 16 hodin než to malinko popostoupilo a našlo vyšší hodnotu.

Máte někdo nápad, jak tohle zrychlit? (Pozn: Už jsem to hodně zrychlil použitím JIT - jedna evaluace MyFunction je teď podle mého odhadu (neměřil jsem) 100x rychlejší než bez JIT - výše uvedeného hodnoty už jsou s JIT)

[Reklama]

[a6b]

nikdo ti nezaruci globalni minimum!

ja mel jen cca 5 parametru a derivaci bych mohl pocitat, ale nechtel jsem slozite, takze jsem pouzival nelder-mead simplex, co derivaci nepotrebuje, akorat musis spocitat funkcni hodnotu v pocet_parametru+1 bodech pro kazdy krok.

[a6b]

prvni video o powellove metode a bylo to jasne. vzdycky je to problem nalezeni vhodneho startovaciho bodu a smerovych vektoru kam dal pokracovat, nikdy nevis ze jsi spadl jen do lokalniho minima.

vylezt z lokalniho minima muze pomoct metoda simulovaneho zihani, ale jistota na globalni minimum neexistuje.

zrychleni jedine snad paralelnim prohledavanim prostoru parametru, tj. paralelne to apoustet s ruznymi startovacimi body a vektory.

[GloGlo]

nikdo ti nezaruci globalni minimum!

ja mel jen cca 5 parametru a derivaci bych mohl pocitat, ale nechtel jsem slozite, takze jsem pouzival nelder-mead simplex, co derivaci nepotrebuje, akorat musis spocitat funkcni hodnotu v pocet_parametru+1 bodech pro kazdy krok.

Derivaci nemusíš počítat ručně, jsou knihovny, který ti ji spočítají - je to hrozně jednoduchý. (Zkoušel sem JAX, funkce grad)

[a6b]

kdyz derivace nejdou, tak bud ten tvuj powell nebo co jsem mel ja simplexovou metodu.
je to podobne.

ale zkus to zkombinovat se simulated annealing, to umoznuje vyhnuti se lokalnimu minimu. pekny priklad je traveling salesman a simulated annealing.

vzdycky jde o to nalezt vhodnou heuristiku jak generovat dalsi body a jak je zhodnotit.

[Reklama]

[Quba]

Bez znalosti toho co optimalizuješ se těžko radí.

Zaměřil bych se na snížení počtu proměnných. Nevím jakou úlohu řešíš, ale optimalizace 500 proměnných bude obecně velmi pomalá. Nelze některé z těch proměnných zafixovat na nějaké rozumné hodnotě a až po nalezení blízkosti extrému tyto proměnné uvolnit a dooptimalizovat?

Jsou opravdu všechny parametry na sobě nezávislé? Nedají se některé vyloučit a dopočítat z ostatních?

Nezkoušel jsi nějaký genetický algoritmus? (Nastavení je zdlouhavé, ale mohlo by to fungovat.  Velmi záleží na typu problému.)

Zajištění globálního extrému stejně nemáš zaručeno. Jediná možnost bude náhodně generovat startovací vektor hodnot a zkoušet optimalizaci.

[a6b]

https://www.google.com/url?q=https://citeseerx.ist.psu.edu/viewdoc/download%3Fdoi%3D10.1.1.722.2134%26rep%3Drep1%26type%3Dpdf&sa=U&ved=2ahUKEwisqbmtrer5AhWRtKQKHai0C3wQFnoECAIQAg&usg=AOvVaw11XOdSnSbbcLxnOJfvkIv8

[a6b]

quba radi dobre, snizit pocet parametru.
nebo pokud maji nektere parametry nejake omezeni (obor hodnot), treba 0-1, to taky muze pomoci.

ja bych masivne paralelizoval, spustit na mnoha strojich s ruznymi starty. nebyl by vypocet v C/C++ rychlejsi nez v pythonu?

je open source vedecka knihovna GSL, tam jsou ruzne metody, nebo Linpack, Lapack.

[GloGlo]

ja bych masivne paralelizoval, spustit na mnoha strojich s ruznymi starty. nebyl by vypocet v C/C++ rychlejsi nez v pythonu?

je open source vedecka knihovna GSL, tam jsou ruzne metody, nebo Linpack, Lapack.

Zkoušel sem všechno, co mě napadlo. Různý starty, různý solvery, různý parametry. S tou metodou, kterou aktuálně používám, se to vždycky někde na dlouho zasekne, než to pak pokračuje - ale to kde se to zasekne je asi dost náhodný, nejde o to, že čím by bylo maximum vyšší, tím by se to sekalo víc.

Jinak ten JIT kompiluje Python podobně jako C/C++, jak sem psal, tak mi to zrychlilo evaluaci funkce odhadem asi 100x.

S tou paralelizací to je zajímavý, jen zatím nevím, jak to udělat.

[Jakub Štech]

Já na tyhle diskrétní/step function problémy rád používám genetiku, stovky až tisíce parametrů jsou spíš ještě menší genotypy. Je normální nechat vyvíjet hotovou neuronovou síť s tisíci neuronů a statisíci synapsemi (např. systém NEAT). Stačí definovat fitness funkci (což už MyFunction přímo je), funkci pro odvození fenotypu z genotypu (což v případě vektoru argumentů pro MyFunction je no-op), a pak už jen pár funkcí pro vytvoření nové generace z předchozí, jednak výběrem nejlepších jedinců, jednak pohlavní kombinací více jedinců a jednak náhodnými mutacemi, kterými to úspěšně překonává lokální maxima.

S tou paralelizací to je zajímavý, jen zatím nevím, jak to udělat.

Spustit pro každé CPU jádro jednu instanci a sbírat výsledky (třeba do souborů) by mělo stačit, ne?

[a6b]

ppripadne to zkusit na grafice, masivni paralelismus tam jde.

[GloGlo]

 

prvni video o powellove metode a bylo to jasne. vzdycky je to problem nalezeni vhodneho startovaciho bodu a smerovych vektoru kam dal pokracovat, nikdy nevis ze jsi spadl jen do lokalniho minima.

vylezt z lokalniho minima muze pomoct metoda simulovaneho zihani, ale jistota na globalni minimum neexistuje.

zrychleni jedine snad paralelnim prohledavanim prostoru parametru, tj. paralelne to apoustet s ruznymi startovacimi body a vektory.

Tohle sem našel na netu. Asi bych to měl nějak zkombinovat s tou scipy.optimize.minimize? Netuším jak.  Chvíli sem si s tím hrál a zkoušel upravovat ty parametry - ale že by to dávalo nějaký zajímavý výsledky samo o sobě, to teda fakt vůbec 

Kód: [Vybrat]

def simAnneal(utility_func, initial_params, data, numSteps=100000,
  noise_magnitude=0.00001, cooling_rate=0.999):
 
optimal_params = initial_params
params = initial_params.copy()  # lists are mutable, so .copy()
best_utility = utility = utility_func(initial_params, data)
temperature = 1.0

for i in range(numSteps):
temperature *= cooling_rate
# consider using numpy/scipy for params and noise
new_params = [param+gauss(0, noise_magnitude)
  for param in params]
# explicitly passing multiple parameters
new_utility = utility_func(new_params, data)

if (new_utility > best_utility
or random.random()*temperature > new_utility / best_utility):
params, utility = new_params, new_utility
if new_utility > best_utility:
optimal_params, best_utility = params, utility
# print (best_utility)
return optimal_params


[Quba]

Já dělám optimalizace fyzikálních úloh, kde jsou parametry vždy omezeny, aby měl výsledek fyzikálně smysl. Jde to i ve tvé úloze?
Mají ty optimalizační parametry nějaký reálný základ, tedy obor hodnot, které mohou nabývat? Nebo jsou hodnoty zcela libovolné.

Asi by bylo vhodné v nalezeném minimu zobrazit 2D řezy té optimalizační funkce. Vybrat dvě proměnné a vykreslit graf v okolí minima. Možná se ukáže, že funkce je "velmi divoká" a v tom případě bude optimalizace dosti obtížná. Asi by stálo za to vygenerovat ty řezy pro všechny proměnné a podívat se které parametry tu "divokost" způsobují a zvážit případně omezení hodnot, které mohou nabývat.

[a6b]

jestli ma 500+ parametru tak si na tech rezech vykouka ocicka :-)

tolik parametru by mohly mit nejake obchodni transakce?

nebo muze hledat nejlepsi konformaci molekuly a 500 parametru jsou uhly mezi atomy v molekule a pocita nejjnizsi konformacni energii.

simulovane zihani jen k jakekoliv optimalizacni metode pridava upraveny test na mensi nalezenou hodnotu, simulovane zihani v zavislosti na teplote jen umoznuje prijmout i vysledek, ktery ma vyssi hodnotu objektivni funkce a umoznuje tak vylezt na kopec z lokalniho minima.

[GloGlo]

Já dělám optimalizace fyzikálních úloh, kde jsou parametry vždy omezeny, aby měl výsledek fyzikálně smysl. Jde to i ve tvé úloze?
Mají ty optimalizační parametry nějaký reálný základ, tedy obor hodnot, které mohou nabývat? Nebo jsou hodnoty zcela libovolné.

V podstatě mám nějakejch 15 základních faktorů, kdy každej má určitou základní "sílu", jenže každej z těch faktorů má několik vlastností, který ovlivňujou celkovou sílu toho faktoru - no a jednak nevím, jak velkou váhu/důležitost ty vlastnosti mají a pak hlavně mají x možných hodnot a každá z těch hodnot zase hraje jinou roli, pravděpodobně v tom není nějaká zákonitost nebo jí minimálně neznám. V podstatě je ale teda pravda, že mám základní faktory a zbytek jsou jen koeficienty, kdy se to může pohybovat v určitém rozmezí, třeba 0.0-1.0 - tak sem podle toho zkusil nastavit ty bounds a je pravda, že se to zrychlilo - jedna iterace teď netrvá kolem 90 sekund, ale kolem 60. Nevýhoda tohoto postupu ale je, že to resetovalo současnej stav a začalo optimalizovat od začátku Tím sem jakoby zahodil cca tejden výpočtů.. Takže to spíš spouštím paralelně a sem zvědavej za jak dlouho ten novej thread dohoní ten starej