PA081: Programování numerických výpočtů 12. Automatické vyhodnocení derivací jaro 2012 Motivace ► kromě funkce dokážeme vyhodnotit i její derivace ► první, druhé, ... ► parciální ► mnohé numerické problémy jsou lépe řešitelné metodami s dostupnou derivací ► metody s derivací zpravidla rychleji konvergují ► Newton vs. metoda sečen Motivace PA081: Programování numerických výpočtů ► kromě funkce dokážeme vyhodnotit i její derivace ► parciální ► mnohé numerické problémy jsou lépe řešitelné metodami s dostupnou derivací ► metody s derivací zpravidla rychleji konvergují ► Newton vs. metoda sečen ► aproximace konečnou diferencí nestačí ► zpracovávaná funkce nebývá triviální ► analytické vyjádření derivace je náročné ► zdroj nebezpečných chyb ► komplikace při změnách ► první, druhé, ... Obrácené vyhodnocení ■O0.O Automatické derivování ► funkce je implementována jako program ► hledáme program, který vedle funkční hodnoty spočte i derivace ► http://www.autodiff.org Automatické derivování Obrácené vyhodnocení 3/16 Automatické derivování PA081: Programování numerických výpočtů ► funkce je implementována jako program ► hledáme program, který vedle funkční hodnoty spočte i I Automatické derivování derivace Obracené ► http://www.autodiff.org ► příklad s majákem bod dopadu paprsku na nábřeží vtancoř yi = —:-- yi = yyi y - tancoř 3/16 Automatické derivování ► vyjádřeno jako jednoduchý program a := tanout; yx := va/(y - a); yi:=yyv, *■ chceme spočítat i rychlost pohybu, tj. derivaci v čase ► výpočet současně s původní funkcí ř := 1 a := tanai ci2 := y - cl yi := va/a2 a\ := tor á\ := tor ä := ái/(cosai)2 áz '■= —d ýi := v(áa2 - aäi)la\ y2 ■= yyi y2 ■= yýi pouze mechanická aplikace základních pravidel Automatické derivování transformace zdrojového kódu (např. ADIC) přetížení operátorů a funkcí (ADOL-C) zjednodušený příklad Automatické derivování Obrácené vyhodnocení class adouble { double x; double dx; } adouble operátor * (adouble a, adouble b) { adouble r; r.x = a.x * b.x; r.dx = a.x * b.dx + a.dx * b.x; return r; } adouble vl, v2, v3; vl = v 2 * v3; Automatické derivování PA081: Programování numerických výpočtů Automatické derivování aktivní proměnná Obrácené vyhodnocení ► nezávislé proměnné, podle kterých chceme derivovat ► všechny další proměnné na nich během výpočtu závisející ► vyhrazený typ adoubl e ► přiřazení adoubl e — doubí e je chybné ► ztrácíme informaci o derivacích tam, kde je potřebná ► explicitní přístup metodou val ue () ► více nezávislých proměnných ► s každou aktivní proměnnou udržujeme vektor parciálních derivací ■O0.O Typický program #define NUMBER_DIRECTIONS 3 #define ADOLC_TAPELESS #include func(adouble[], const adouble[]); adouble x[3]; adouble y[4] ; for (i=0; i<3; for (j=0; i<3; x[i] .setADValueCj ,i==j ? 1.0 : 0.0); func(y.x); Automatické derivování Obrácené vyhodnocení Není to tak jednoduché ► nediferencovatelné funkce (fmod) ► nelze použít ► absolutní hodnota (f abs) ► v 0 pravá nebo levá derivace ±1 ► případně položíme uměle rovnu 0 minimum a maximum (fmi n, fmax) ► lze nahradit (a + b ± \a - b\)/2 ► nedefinované na celém R (sqrt, pow, tan, ► derivaci položíme uměle NaN, resp. Inf Obrácené vyhodnocení ► množství potřebných operací O(mk) pro /: R ► m počet nezávislých proměnných ► n počet závislých proměnných ► fc počet kroků výpočtu ► nepříjemné s rostoucím m *■ velká třída problémů velkým m ale malým n ► optimalizace funkce více proměnných (n = 1) Obrácené vyhodnocení PA081: Programování numerických výpočtů ► množství potřebných operací O(mk) pro /: Rm —■ Rn ► m počet nezávislých proměnných ► n počet závislých proměnných ► fc počet kroků výpočtu ► nepříjemné s rostoucím m *■ velká třída problémů velkým m ale malým n ► optimalizace funkce více proměnných (n = 1) ► počítali jsme citlivost pomocných proměnných na změny vstupu ► obrácený postup, citlivost výstupu (závislých proměnných) na změny pomocných ► výpočet od konce Automatické derivování Obrácené vyhodnocení Obrácené vyhodnocení PA081: Programování numerických výpočtů např. funkce y = sin(xi ► označujeme ä = By/Ba /x2)eX3 Automatické derivování Obrácené vyhodnocení x2 := xilx2 a2 '■= sinai a 3 := eX3 y := a2a3 y:= 1 ä2 := a3y ň-i := a2y x3 := a3ä3 ä\ := ä2 cosai := -ä\/x2 := ä2/x2 ~)(nk) operacích máme všechny By/Bxi Stopa výpočtu derivace se vyhodnocují od konce ► tj. až poté, co proběhl výpočet funkce ► proto „obrácené" (reverse) vyhodnocení výpočet funkce je nutné nějak zaznamenat ADOL-C: datová struktura (i soubor) páska (tape) záznam posloupnosti operací, nikoli konkrétních hodnot ► lze recyklovat pro jiný vstup ► případné rozdílnosti vyhodnocení (a>b ? c : d)jsou detekovány použije se pro vyhodnocení funkce, gradientu/Jakobiánu, Hessiánu, ... Automatické derivování Obrácené vyhodnocení Instrumentace pro zpětné vyhodnocení #include void func( double px[4], /* vstup */ double py[4], /* vystup */ ) { adouble x[4],y[4],aux; trace_on(0); /* paska 0 */ f oř (i=0; i<4; x[i] «= px[i]; aux = x[l] * x[4] + exp(x[2]); /* vlastni vypočet */ y[4] = si n(aux); f oř (i=0; i<4; y[i] py[i]; trace_off(); } Automatické derivování Obrácené vyhodnocení Drivery ► ADOL-C poskytuje celou řadu „driver" funkcí ► vstupem je vždy páska a konkrétní vstupní hodnoty ► driver provede výpočet na základě záznamu na pásce ► funkční hodnota funkce Rn —■ Km ► function(short tag, int m, int n, double x[n],double y[m]) ► gradient funkce W1 -* K ► gradient(short tag, int n, double x[n],double g[n]) Jakobián funkce Rn - Rm ► jacobian(short tag, int m, int n, double x[n],double J[m][n]) ► Hessián (matice druhých derivací) funkce Rn —■ R ► hessian(short tag, int m, int n, double x[n],double H[n][n]) Drivery ► optimalizované varianty pro Newtonovu metodu ► řešení různých typů diferenciálních rovnic ► speciální varianty pro řídké Jakobiány a Hessiány ► tenzory vyšších derivací ► implicitně definované funkce G(x,y) = 0 ► obecné - plná kontrola nad zpracováním pásky Rekapitulace obráceného vyhodnocení PA081: Programování numerických výpočtů Automatické derivování Obrácené ► aktivní proměnné deklarovat jako adoubl e vyhodnocení ► výpočet označit trace_on () ... trace_off () ► zavolat instrumentovanou funkci jednou ► s typickými hodnotami vstupů ► vyprodukuje záznam výpočtu - pásku ► opakovaně volat potřebné drivery ► případně na různé vstupy ► pro jinou posloupnost výpočtu zavolat instrumentovanou funkci s jinou páskou Další materiál ► rozsáhlá problematika, řada implementací ► včetně obsáhlé teorie ► http://www.autodiff.org ► přehled nástrojů, FAQ, ... https://projects.coi n-or.org/ADOL-C ► implementac a dokumentace ADOL-C ► včetně řady příkladů ► Griewank, Walter. Evaluating Derivatives. 2008