PA081: Programování numerických výpočtů
2. Numerická stabilita (nejen) elementárních výpočtů
Aleš Křenek
jaro 2011
Kvadratické rovnice
PA081: Programování numerických
výpočtů
A. Křenek
► rovnice tvaru
► školní vzorec
ax2 + bx + c = 0
-b ± \Jb2 - 4ac 2a
Kvadratické rovnice
Elementárni funkce
Algebraické úpravy
Mocninné řady
Odmocnina
Algebra kvaternionu
2/33
Kvadratické rovnice
► rovnice tvaru
ax2 + bx + c = 0
► školní vzorec
-b ± \Jb2 - 4ac
x± =-~-
za
je-li b2 » ac, počítáme x+ z rozdílu dvou velmi blízkých čísel
► problém katastrofálního zrušení
► „»" neznamená až příliš velký rozdíl koeficientů
► typ float - přesnost na 7-8 dekadických číslic
► při a = 1 tedy stačí i? - 103c
PA081: Programování numerických
výpočtů
A. Křenek
Kvadratické rovnice
Elementární funkce
Algebraické úpravy
Mocninné řady
Odmocnina
Algebra kvaternionů
2/33
Kvadratické rovnice
► konkrétně pro b = 2000, c = 1
► „přesné" hodnoty
► VD = 1999.9989999997499998749999218749453...
► x+ = -.00050000012500006250003906252734377...
► x- = -1999.9994999998749999374999609374...
► pro f 1 oat
► VD = 1999.999
► x+ = -0.00048828125, chyba 2.5%
► X- = -1999.9995, chyba odpovídá přesnosti typu
Kvadratické rovnice
► numericky stabilní řešení pro x+
► vlastnosti kořenů rovnice:
ax2 + bx + c = (x - x+)(x - X-)   tedy  c = x+X-
a lze počítat x+ = c/x-\
► pro předchozí příklad dostáváme
► x+ = -0.00050000014, chyba odpovídá přesnosti typu
► analogicky pro b < 0 je nepřesné X-, vypočte se symetricky
Kvadratické rovnice
► numericky stabilní řešení pro x+
► vlastnosti kořenů rovnice:
ax2 + bx + c = (x - x+)(x - X-)   tedy  c = x+X-
a lze počítat x+ = c/x-\
► pro předchozí příklad dostáváme
► x+ = -0.00050000014, chyba odpovídá přesnosti typu
► analogicky pro b < 0 je nepřesné X-, vypočte se symetricky
► i v triviálním výpočtu může vzniknout problém
► řešení může být docela jednoduché
Elementární funkce
PA081: Programování numerických
výpočtů
A. Křenek
► goniometrické, hyperbolické, exponenciální, logaritmické, odmocniny, ...
Kvadratické rovnice
Elementární funkce
Algebraické úpravy
Mocninné řady
Odmocnina
Algebra kvaternionů
5/33
Elementární funkce
► goniometrické, hyperbolické, exponenciální, logaritmické, odmocniny, ...
► pohled matematika
► zaběhaný intuitivní aparát
► příjemné analytické vlastnosti (spojitost, derivace, ...)
► většina matematických modelů se bez nich neobejde
PA081: Programování numerických
výpočtů
A. Křenek
Kvadratické rovnice
Elementární funkce
Algebraické úpravy
Mocninné řady
Odmocnina
Algebra kvaternionů
5/33
Elementární funkce
► goniometrické, hyperbolické, exponenciální, logaritmické, odmocniny, ...
► pohled matematika
► zaběhaný intuitivní aparát
► příjemné analytické vlastnosti (spojitost, derivace, ...)
► většina matematických modelů se bez nich neobejde
pohled programátora
► samy o sobě numericky stabilní
► optimalizované implementace v knihovnách
► problematické chování v okrajových případech vede na numerickou nestabilitu
► netriviální výpočetní náročnost
PA081: Programování numerických
výpočtů
A. Křenek
Kvadratické rovnice
Elementární funkce
Algebraické úpravy
Mocninné řady
Odmocnina
Algebra kvaternionů
5/33
Elementární funkce
Co je špatně?
► celý aparát vede na iracionální čísla (tt, e, V2,...)
► ve většině případů zdroj nepřesnosti
► ex
► tendence k přetečení - e88 = 1.65 x 1038 tanx
► přetečení pro x -* j (a perioda) roste nade všechny meze
► sinx a cosx
► pro x -* y resp. x -* 0 vedou na špatně podmíněné rovnice
► malá změna vstupu vede k velké změně výstupu
► např. sinx = ř pro ř ->• 1
Kvadratické rovnice
Elementární funkce
Algebraické úpravy
Algebra kvaternionů
6/33
Co s tím?
► uvědomit si možné problematické chování
► potřebuji skutečně tyto funkce pro svůj výpočet?
► mnoho problémů má i jiné řešení
► provést transformace, kterými se vyhneme numerické nestabilitě
► podobné (a + b)2 = a2 + 2ab + b2
► podle kontextu zvolit vhodnou implementaci
Algebraické úpravy
► výraz VI + x - VI - x
pro malá x odčítání velmi blízkých čísel
► transformace
+ X
vr
X
(VI + x - VI - *)(Vl + x + VI - x)
vTTx + VT (1 +x) - (1 -
- X X)
VI + x + Vi
2x
VI + x + Vl - x ► sčítání velmi blízkých čísel - dostatečně přesné
Kvadratické rovnice
Elementární funkce
Algebraické úpravy
Algebra kvaternionů
8/33
Algebraické úpravy
► výraz ln \Jx + 1 - ln ^fx
► pro velká x rozdíl blízkých čísel
► transformace
ln     + 1 - ln V-*7 =
= ln(x + 1)^ - lnx^ = ^(ln(x + 1) - lnx)
= Trln-= -ln(l + -)
2x2 x
► možná ztráta přesnosti, ale i tak lepší
PA081: Programování numerických
výpočtů
A. Křenek
Kvadratické rovnice
Elementárni funkce
Algebraické úpravy
Mocninné řady
Odmocnina
Algebra kvaternionú
9/33
Mocninné řady
► Taylorův rozvoj základních funkcí
eX = l + Y\ + fr + --
sinx
x3 x5
Kvadratické rovnice
Elementární funkce
Algebraické úpravy
cosx
ln(l + x)
Vl + x 1
VI + x
x^ x^ + 4!
x2 x3
i +
X
2 x
2 + 2~A
x^ 2~A
3x2
3x3
2-4-6
3 ■ Sx3 2-4-6
Algebra kvaternionů
10/33
Mocninné řady
Demonstrační příklad
► vyraz
1 - cosx
se pro x ->• 0 blíží \
► pro malá x odečtení dvou blízkých čísel
► konkrétně pro x = 0.0005 ve float:
cos* = 0.99999988 1 - cosx = 1.1920928 x 10~7   (falešná přesnost)
1 ~C°SX = 0.47683709
xl
► správný výsledek je 0.49999999 (na 8 míst)
Mocninné řady
Demonstrační příklad
► náhrada cos x Taylorovým rozvojem
X2      X4 X6
cosx = 1---1-----h...
2!     4! 6!
► pro výpočet na 8 cifer stačí do řádu x4
i i i i x2 x4 1 ? 1 - COS X _ í - í + ~2--12  _ 1 X
x~2      ~ x~2 ~ 2 ~ 12
► výsledek výpočtu pro x = 0.0005 je 0.49999997
► podstatně lepší přesnost
řádově méně aritmetických operací
Mocninné řady
Nedostatky
► rychlá konvergence pro malá x, pomalá pro větší
► přímý důsledek Taylorovy věty
► vyplatí se použít vlastní vzorec pro malé x
► resp. blízko problematického bodu numerické stability
► v ostatních případech zůstat u knihovní funkce
► jak poznáme, co je „malé x"?
Mocninné řady
Kdy použít
► známe požadovanou přesnost
► bezpečné je použít přesnost daného typu
► pro f 1 oat a výsledky - 1 je to 10~7
► při méně přesných vstupních datech méně striktní
► nerovnost x < c musí zajistit, že největší zanedbaný člen řady je menší než požadovaná přesnost
PA081: Programování numerických
výpočtů
A. Křenek
Kvadratické rovnice
Elementární funkce
Algebraické úpravy
Mocninné řady
Odmocnina
Algebra kvaternionů
14/33
Mocninné řady
Kdy použít
►
známe požadovanou přesnost
► bezpečné je použít přesnost daného typu
► pro f 1 oat a výsledky ~ 1 je to 10~7
► při méně přesných vstupních datech méně striktní
nerovnost x < c musí zajistit, že největší zanedbaný člen řady je menší než požadovaná přesnost
konkrétně v předchozím příkladu x4
— < IQ'7   tedy  x < 0.09212 o!
zkontrolujeme, zdali dává smysl pro problematickou funkci
cosO.09212 = 0.99999871 je ještě v toleranci
Mocninné řady
Kdy použít
výsledný kód pro
1 - cosx
%2
Kvadratické rovnice
Elementární funkce
Algebraické úpravy
float x, y;
if (x < 0.09212)
y = 0.5 + x-x/12.0;
el se
y = (l-cos(x))/(x*x);
Algebra kvaternionů
15/33
Odmocnina
► ve většině používána k normalizaci vektorů
► silové působení, např. dva bodové náboje a a
► velikost síly
F = kcqaqb
kde r = ||a - b|| silový vektor
r
►
snadný výpočet
r2 = X(«í-kť)2 r už vyžaduje odmocninu
Odmocnina
Masivní použití
počítačová grafika - osvětlení plochy
► zpravidla interpolací počítáme normály k zakřivené ploše
► vektory mají správný směr ale nejsou jednotkové
► pro správný výpočet osvětlení potřebná normalizace
ke všem těmto výpočtům nepotřebujeme striktně -Jx J= se hodí mnohem více
► násobení je lepší než dělení
navíc zpravidla stačí velmi hrubá aproximace
► citlivý je hlavně směr vektoru, ten zůstává
přímo Taylorův rozvoj -yj= iterační metody (Newtonova atd.)
Kvadratické rovnice
Elementární íunkce
Algebraické úpravy
Algebra kvaternionů
17/33
Odmocnina
Fast inverse square root
► hra Quake III na konci 90. let
► převratně realistická grafika
► vděčila velmi rychlému výpočtu J=
► kód zřejmě pochází z grafických knihoven SGI
Odmocnina
Fast inverse square root
► hra Quake III na konci 90. let
► převratně realistická grafika
► vděčila velmi rychlému výpočtu J=
► kód zřejmě pochází z grafických knihoven SGI
Odmocnina
Fast inverse square root
► funguje na základě IEEE reprezentace čísla
x = (1 + mx) x 2ĚX
sign exponent (8 bits) significand (23 bits)
I I II
0	0	1	1	1	1	1	0	0	0	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
31 30 23 2°2 (bit index) Ó
► přitom mx = Mx/223 a ex = Ex - 127
► základ výpočtu
1
log2y = --log2x log2(l + my) + ey = -^log2(l + mx) -
Odmocnina
Fast inverse square root
► pro m e [0,1) si lze dovolit aproximaci
log2(l + m) = m + a
► dosazením a úpravami dostaneme
Eyx 223 +My=R- X 223 + Mx)
tj. červený řádek v kódu s empiricky stanovenou konstantou R
► více viz http ://en. wi ki pedi a .org/wi ki/Fast_ i nverse_square_root
Odmocnina
Fast inverse square root
poslední krok je jedna iterace Newtonovy metody
y'
► pro vstup x hledáme takové y aby y
► tj. hledáme kořen funkce f(y) - 1 Newtonovou metodou
v   ! - v     fiyn)     v -Í
X
X
y n +
n f'(yn)
Jn
což už odpovídá poslednímu výrazu v kódu
yn(l - xyp 2
Kvadratické rovnice
Elementární funkce
Algebraické úpravy
Algebra kvaternionů
21/33
Odmocnina
Fast inverse square root
► aproximace je překvapivě přesná
► v Quake III je zakomentována další Newtonova iterace
► nebyla potřeba
► cca. 4x rychlejší než dělení
► dnes překonána instrukcí rsqrtss
Odmocnina
Fast inverse square root
► aproximace je překvapivě přesná
► v Quake III je zakomentována další Newtonova iterace
► nebyla potřeba
► cca. 4x rychlejší než dělení
► dnes překonána instrukcí rsqrtss
► přesné pochopení konkrétního účelu výpočtu
► včetně zhodnocení „má smysl tvrdě optimalizovat"
► volba adekvátní metody
► přispěla k velkému komerčnímu úspěchu
► pro jiné účely by byla zcela nevhodná
Odmocnina
Lze se jí i vyhnout
► Lenard-Jonesuv potenciál
► nevazebná interakce dvou atomů
► významná pro modelování biochemických dějů
► vzorec pro energii (potenciál)
A     B   , J Elj = —jTj---g   kde r je vzdálenost atomu
Odmocnina
Lze se jí i vyhnout
► výpočet velikosti síly
_ dE _ _12A 6B dr       r13 r7
liché exponenty => bude potřeba odmocnina
Odmocnina
Lze se jí i vyhnout
►
výpočet velikosti síly
_ dE _ _12A 6B dr       r13 r7
► liché exponenty => bude potřeba odmocnina
► zajímá mě většinou silový vektor
ca-b    ,       , (6B VIA
F-= (a-b) —r
r
-14
►
tedy vystačím s r
2
Odmocnina
Lze se jí i vyhnout ► implementace
float n,r2,r4,r8,rl4,F[3]; i nt i ;
for (n=0,i=0; i<3; i++) { F[i] = a[i] - b[i]; n += F[i]*F[i];
}
r2 = 1.0/n;
r4 = r2*r2; r8 = r4*r4; r14 = r8*r4*r2;
n = 12*A*rl4 - 6*B*r8; for (i=0; i<3; i++) F [i]
••= n;
Kvadratické rovnice
Elementární funkce
Algebraické úpravy
Algebra kvaternionů
Parametrizace rotace
Aneb jak se vyhnout goniometrickým funkcím
► ve 2D jeden úhel <j>
► složení součtem, inverze -</> aplikace násobením maticí
^    . cos cj> -sinej) sinej)    cos cj)
goniometrickým funkcím se můžeme vyhnout
■ ^      21             * 1-t2 sm<p =--~      cosep =--~
► numericky příjemné, nevyjádříme cj> = ±j
Parametrizace rotace
► ve 3D eulerovské úhly
► značně nepřehledné
► větší množství singularit
► maticové vyjádření
► složení 3 rotací podél osy
/ 1      0 0 \
0   cos/? -sínfi \ 0   sin/?    cos/? /
aplikovatelná stejná substituce ► vede na komplikované polynomy
Parametrizace rotace
► přímo maticí 3x3
► numericky samo o sobě stabilní, pokrývá všechny případy
► příliš mnoho „parametrů navíc"
► musí být ortogonální
► snadno degeneruje
► nepřesností výpočtu ztrácí ortogonalitu korekce není jednoduchá
► zabere více paměti
Kvaterniony
► více ekvivalentních definic
nejjednodušší U = R.4 s kanonickou bází (1, i, j, k) a násobením
i2 = f = k2 = íjk = -1
► chápeme R.3 c Hl jako ryze imaginární kvaterniony
► potom pro \q\ = 1 je zobrazení x >-> qxq ortogonální lineární transformace
Kvaterniony
více ekvivalentních definic
► nejjednodušší U = R.4 s kanonickou bází (1, i, j, k) a násobením
i2 = f = k2 = íjk = -1
► chápeme R.3 c U jako ryze imaginární kvaterniony
► potom pro \q\ = 1 je zobrazení x >-> qxq ortogonální lineární transformace
► nakrytí 2:1 (q a -q reprezentují stejnou transformaci)
► skládání rotací je násobení
► inverze je q
► korekce degenerace normalizací
► snadná interpolace (SLERP)
Superpozice množiny bodů
Formulace problému
► množiny odpovídajících si bodů {Xí} a {yť} v R.3
► hledáme transformaci Q (rotaci + posunutí) tak, aby
XIIQxí-VíH2
bylo minimální
posunutí - stačí nezávisle posunout těžiště {Xí} do těžiště {yť}
► zbývá najít rotaci
Superpozice množiny bodů
Použití kvaternionů
► formulace pomocí kvaternionů
X WqJkq ~ YtW2 = XdlXíll2 + IIVíII2) - 2 ^Re^M)
► potřebujeme maximalizovat Y,^(YíH^-íH)
► lze vyjádřit jako kvadratickou formu U ->■ R.
► funkce q >-+ qTAq, kde A je symetrická matice 4x4
► prvky A lze vyjádřit jako lineární kombinace prvků matice
Superpozice množiny bodů
Nalezení maxima
► podle věty o hlavních osách existuje P tak, že A = PAPT kde P je ortogonální a A diagonální
► prvky A jsou vlastní hodnoty A
► řádky P jsou vlastní vektory A
► potom
q1Aq = Mql + \2O-2 + ^3^3 + ^4HÍ
kde qt jsou souřadnice q v bázi tvořené vlastními vektory A
protože báze je ortonormální, platí na základě \\q\\ = 1 S^2 = l
► předpokládejme např. Ai největší
► potom výraz je maximální pro q = (1,0,0,0)
Superpozice množiny bodů
Rekapitulace
► vyřešíme posun bodů do těžiště
► spočítáme matici koeficientů X X/YÍ
► z nich matici kvadratické formy
► najdeme nej větší vlastní vektor této matice
► na to existují rychlé metody
► pro 4 x 4 v podstatě triviální
PA081: Programování numerických
výpočtů
A. Křenek
Kvadratické rovnice
Elementární funkce
Algebraické úpravy
Mocninné řady
Odmocnina
Algebra kvaternionů
33/33
Superpozice množiny bodů
Rekapitulace
► vyřešíme posun bodů do těžiště
► spočítáme matici koeficientů X x{yf
► z nich matici kvadratické formy
► najdeme nej větší vlastní vektor této matice
► na to existují rychlé metody
► pro 4 x 4 v podstatě triviální
► práce s rotacemi je velmi frekventovaná
► grafika, fyzikální modely, ...
► kvaterniony jsou numericky velmi vhodná reprezentace
► transformace zdánlivě nesouvisejících problémů
► superpozice vs. hledání vlastních vektorů
► efektivní numerické řešení dané úlohy