Minimalizace součtu čtverců- úvod

Optimalizuje teoretický model tak, aby co nejvíce odpovídal naměřeným datům.

=> Minimalizuje odchylku modelu od naměřených hodnot.

Využití:

Všude, kde máme co do činění s analýzou nějakého přírodního nebo technického systému.

Minimalizace součtu čtverců- úvod II

Naměřená data si můžeme představit jako dvojice:

(ti, yi), i = 1, ..., m

kde:

ti Rk bod mìøení (napøíklad èas nebo místo mìøení

nebo obojí)

yi hodnota, namìøená v ti

Minimalizace součtu čtverců- úvod III

Dále pak máme nìjaký matematický model M:

Rk+n -> R, který je závislý na n volných parametrech x1, x2, ..., xn a pro který požadujeme, aby:

M(ti, x) yi

kde:

x = (x1, ..., xn)

i = 1, ..., m

(m je tedy počet naměřených bodů, se kterými

budeme pracovat)

Minimalizace součtu čtverců- úvod IV

V úlohách tohoto typu tedy pro m-prvkovou množinu naměřených bodů (ti, yi) hledáme parametry x1,..., xn modelu M tak, aby daný model co možná nejlépe popisoval tuto množinu.

=> Minimalizujeme odchylku modelu od naměřených dat.

Minimalizace součtu čtverců- příklad

Ohmùv zákonData: ((Ui), Ii) kde Ui je napìtí na svorkách

rezistoru a Ii je proud, který prochází rezistorem

Model: Obecně: M(ti, x) pro data (ti, yi)

Konkrétně:

R

U)R(),U(M i

i

Parametry modelu: x = (R), kde R je odpor rezistoru.

Minimalizace součtu čtverců- příklad II

Radioaktivní rozpadData: ((ti), Ni) kde ti je èas od poèátku

mìøení a Ni je poèet atomù v èase ti

Model: Obecně: M(ti, x) pro data (ti, yi)

Konkrétně:Parametry modelu: x = (N0, T), kde N0 je poèet

atomù v èase 0 a T je poloèas rozpadu.

2ln.T

t

00i

i

e.N)T,N(),t(M

Minimalizace součtu čtverců- příklad III

Potenciální energie molekulyData: ((s1,..., sn), Epot) kde s1, ..., sn jsou souřadnice

jednotlivých atomů molekuly a Epot je potenciální energie molekuly

Model: Obecně: M(ti, x) pro data (ti, yi)

Konkrétně:

Parametry modelu: je jich velmi mnoho – ideální vazebné vzdálenosti, vazebné úhly a torzní úhly, konstanty úměrnosti, atd ...

elvdwtoranbnn1 EEEEE(...)),s,...,s(M

Minimalizace součtu čtverců- obecně

Chceme minimalizovat odchylku modelu od naměřených dat =>

Chceme tedy, aby hodnoty rozdílù

ri(x) = M(ti, x) - yi

byly v absolutní hodnotì co nejmenší.

To se dá interpretovat jako minimalizace normy vektoru:

r(x) = (r1(x), ..., rm(x))T

Minimalizace součtu čtverců- obecně II

Nejèastěji se používá euklidovská (L2) norma, pro kterou dostáváme minimalizovanou funkci ve tvaru:

Namísto L2 normy je také možno použít normu L1 (souèet absolutních hodnot ri) nebo L (maximum z absolutních hodnot ri). Tyto normy mají svoje opodstatnìní: napøíklad L1 norma lépe eliminuje body mìøení, které „uletìly“, tj. jsou výraznì mimo prùbìh zadaný ostatními body, èasto v dùsledku chyby pøi mìøení. V pøednáškách budeme dále pracovat pouze s Euklidovskou normou.

m

1i

2i

T )(r)()()(f xxrxrx

Minimalizace součtu čtverců- obecně III

Na funkce typu:

je možno přímo aplikovat minimalizační metody, probrané v předchozích kapitolách. Při výpočtech ale můžeme ušetřit mnoho času i paměti tím, že využijeme speciální vlastnosti tohoto problému :-).

m

1i

2i

T )(r)()()(f xxrxrx

Minimalizace součtu čtverců- obecně IV

Gradient funkce f se dá vyjádøit jako:

kde J(x) Rnm je Jakobiho matice funkce f(x)

(i-tý øádek matice J(x) je gradient ri v bodì x).

Hessián funkce f se pak dá zapsat jako:

)(.)(2)(r).(r2)( Tm

1iii xrxJxxxg

G x x x x x( ) ( ). ( ) ( ). ( ) 2 2

1

2

1

r r r ri iT

i

m

i ii

m

Minimalizace součtu čtverců- obecně V

Pokud je model M v dobré shodì s daty, pak v minimu x* má funkce f velmi malou (kladnou) hodnotu. Pak tedy ri(x) jsou malá èísla a v okolí x* se proto dá zanedbat druhá suma ve vztahu:

Hessián funkce f pak může být aproximován jako:

G x x x x x( ) ( ). ( ) ( ). ( ) 2 2

1

2

1

r r r ri iT

i

m

i ii

m

)()(2)(r).(r2)( Tm

1i

Tii xJxJxxxG

Minimalizace součtu čtverců- Gauss-Newtonovy metody

Metody, které kombinují tuto aproximaci s Newtonovskými metodami se nazývají Gauss-Newtonovy metody.

Klasický Newtonovský vztah pro výpočet směru přesunu (sk) z bodu xk do bodu xk+1:

je tedy přeformulován na vztah:

)k()k()k( . gsG

J J s J r( ) ( ) ( ) ( ) ( ). . .k T k k k T k

Minimalizace součtu čtverců- Levenberg-Marquardtovy metody

Použití dané aproximace v metodě s omezeným krokem dává Levenberg-Marquardtovu metodu pro součet čtverců. Původně byla Levenberg-Marquardtova metoda vyvinuta právě pro tuto aplikaci. Rovnice:

kterou využívá Levenberg-Marquardtova metoda se v tomto speciálním případě přepisuje do tvaru:

 

G I g( ) ( ) ( )k k k

J J I J r( ) ( ) ( ) ( ) ( ). .k T k k k T k

Minimalizace součtu čtverců- další metody

Situace se komplikuje v případě, kdy neexistuje dostatečně dobré minimum funkce f, tj. když r(x) je nezanedbatelně velké pro všechna x.

Pak je možné například kombinovat metodu BFGS a Gauss-Newtonovu metodu tímto způsobem:

Podle velikosti relativního poklesu funkce f v k-té iteraci se hessián pro další iteraci aproximuje buď pomocí BFGS nebo pomocí Gauss-Newtonovy metody.

Minimalizace součtu čtverců- další metody II

Pøesnìji øeèeno, pokud

f(k) - f(k+1) .f(k)

kde (0, 1) je vhodì zvolený parametr, pak se pro odhad použije Gauss-Newtonova metoda, v opaèném pøípadì pak BFGS. Pro hodnotu doporuèuje napø. Fletcher volit = 0.2.

Lineární úloha nejmenších čtverců- úvod

V tomto případě je model lineární vzhledem k aproximovaným parametrùm:

M(ti, x) = 1(ti).x1 + ... + n(ti).xn

Pro odchylku modelu od reálného výsledku měření platí:

=> ri(x) = M(ti, x) - yi = 1(ti).x1 + ... + n(ti).xn - yi

Funkce, kterou budeme v rámci metody minimalizovat, má tedy tvar:

f r t x t x yii

m

i n i n ii

m

( ) ( ) ( ). ... ( ).x x 2

11 1

2

1

Lineární úloha nejmenších čtverců- úvod II

Budeme tedy minimalizovat funkci:

V minimu musí pro všechny parametry x1, ..., xn modelu platit:

Po odderivování tedy platí:

f

x j

0

m

1i

2inin1i1 yx).t(...x).t()(f x

f

xt t x t x y

jj i i n i n i

i

m

2 01 1

1

. ( ). ( ). ... ( ).

Lineární úloha nejmenších čtverců- úvod III

Rovnici:

budeme dále upravovat:

Soustavu rovnice v tomto tvaru můžeme zapsat pomocí matice:

A.x = b

m

1iinin1i1ij 0yx).t(...x).t().t(.2

)t(.x).t(...x).t()t(.y ij

m

1i

m

1inin1i1iji

m

1iijinn

m

1iiji11 )t().t(.x...)t().t(.x

Lineární úloha nejmenších čtverců- úvod IV

Soustavu rovnic:

lze zapsat ve tvaru A.x = b následovně:

kde k, j {1, …, n}

Můžeme tedy obejít náročný proces minimalizace a získat minimum přímo řešením této soustavy.

m

1iijinn

m

1iiji11

m

1iiji )t().t(.x...)t().t(.x)t(.y

m

1iijikkj )t().t(a

m

1iijik )t(.yb

Lineární úloha nejmenších čtverců- příklad

Chceme řešit tento problém:

Mějme objekt, pohybující se v čase t rychlostí v.

Naměřili jsme, že objekt se v čase t1 = 1s pohyboval rychlostí v1 = 1 m/s

t2 = 2s pohyboval rychlostí v2 = 6 m/s

t3 = 3s pohyboval rychlostí v3 = 2 m/s

Pohyb tohoto objektu považujeme za rovnoměrně zrychlený a můžeme ho tedy modelovat pomocí vztahu:

v = a.t + v0, kde:

a zrychlení objektu

v0 počáteční rychlost objektu

Úkolem je určit parametry a a v0.

Lineární úloha nejmenších čtverců- příklad II

Obecný vztah pro model:

M(ti, x) = 1(ti).x1 + ... + n(ti).xn

můžeme tedy v našem případě přepsat do tvaru:

M(ti, (a, v0)) = 1(ti).a + 2(ti).v0 = a.t + v0

=> n = 2

Pro 1 a 2 tedy platí:

1(ti) = ti

2(ti) = 1

Lineární úloha nejmenších čtverců- příklad III

Měřením jsme získali 3 body: (1, 1); (2, 6) a (3, 2) (=> m = 3)

Problém lze obecně zapsat pomocí soustavy n rovnic A.x = b

kde k, j {1, …, n}

m

1iijikkj )t().t(a

m

1iijik )t(.yb

V našem případě tedy bude platit: a11 = t1. t1 + t2. t2 + t3. t3 = 14

a12 = t1 + t2 + t3 = 6

a21 = t1 + t2 + t3 = 6

a22 = 1 + 1 + 1 = 3

b1 = y1. t1 + y2. t2 + y3. t3 = 19

b2 = y1 + y2 + y3 = 9

x1 = a

x2 = v0

Lineární úloha nejmenších čtverců- příklad IV

Budeme tedy řešit soustavu rovnic:

Řešení: a = 0.5 m/s2

v0 = 2 m/s

Součet čtverců odchylek:

9

19

36

614

3

1i

2i0i

m

1i

2ii yvt.ay),t(MS x

5,135,1)3(5,1

)223.5,0()622.5,0()121.5,0(222

222

Lineární úloha nejmenších čtverců- příklad V

Grafické znázornění výsledku:

y = 0,5x + 2

0

1

2

3

4

5

6

7

0 0,5 1 1,5 2 2,5 3 3,5

t [s]

v [

m/s

]

Lineární úloha nejmenších čtverců- úlohy se dvěma parametry

S tímto typem úloh se setkáváme v praxi velmi často, jedná se o úlohy typu:

Máte zadáno m bodů (ti, yi), proložte těmito body přímku.

= nalezněte koeficienty k a q v rovnici y = k.t + q.

V tomto případě lze obecnou soustavu rovnic A.x = b

kde k, j {1, …, n}

Pøepsat do tvaru:

m

1iijikkj )t().t(a

m

1iijik )t(.yb

m

1ii

m

1iii

m

1ii

m

1ii

m

1i

2i

y

t.y

q

k.

mt

tt

Lineární úloha nejmenších čtverců- úlohy se dvěma parametry II

Řešením této soustavy:

Pak získáme pro k a q vztahy:

m

1ii

m

1iii

m

1ii

m

1ii

m

1i

2i

y

t.y

q

k.

mt

tt

m

t.tt

m

y.tt.y

k m

1ii

m

1iim

1i

2i

m

1ii

m

1iim

1iii

m

t.kyq

m

1ii

m

1ii

Lineární úloha nejmenších čtverců- úlohy se dvěma parametry III

Příklad:

Zadání stejné jako u předchozího příkladu o „objektu, pohybujícím se v čase t rychlostí v“.

Máme tedy body: (1, 1); (2, 6) a (3, 2) a chceme jimi proložit přímku: v = a.t + v0

Konkrétní výpočet:

19t.v i

3

1ii

6t3

1ii

9v3

1ii

14t3

1i

2i

5,0

36.6

14

39.6

19k

3

6.5,09q

Kvadratická úloha nejmenších čtverců- úloha se třemi parametry

Naměřenými body tedy chceme proložit rovnici:

y = a.t2 + b.t + c

Analogicky jako v lineárním případě lze i v kvadratickém případě tuto speciální úlohu zapsat pomocí soustavy rovnic A.x = b, a to následovně:

t t t

t t t

t t cm

a

b

c

t y

t y

y

i i i

i i i

i i

i

i

i

m

i

m

i

m

i

m

i

m

i

m

i

m

i

m

ii

m

ii

m

ii

m

4

1

3

1

2

1

3

1

2

1 1

2

1 1

2

1

1

1

.

Kvadratická úloha nejmenších čtverců- úloha se třemi parametry II

Metodou nejmenších čtverců najděte polynom 2.stupně, který je nejblíže bodům:

[1,1], [2,3], [4,6].

Řešíme tedy soustavu rovnic:

t t t

t t t

t t m

a

b

c

t y

t y

y

i i i

i i i

i i

i

i

i

m

i

m

i

m

i

m

i

m

i

m

i

m

i

m

ii

m

ii

m

ii

m

4

1

3

1

2

1

3

1

2

1 1

2

1 1

2

1

1

1

.

Kvadratická úloha nejmenších čtverců- úloha se třemi parametry III

=> soustava: 273a + 73b + 21c = 109

73a + 21b + 7c = 31

21a + 7b + 3c = 10

Výsledek je:

a = -1/6,

b = 5/2,

c = -4/3.

y t t

1

6

5

2

4

32

Cvičení- příklad 1

Metodou nejmenších čtverců najděte přímku, která je nejblíže bodům: (1,1); (2,6) a (3,2)

Použijte „přímý přístup“ - vytvořit funkci, odderivovat, derivaci položit rovnu nule, dořešit :-).

Řešení: Hledáme přímku ve tvaru y = kx + q

Minimalizovaná funkce:

f k q t k q yi ii

n

( , ) ( . )

1

2

f k q k q k q k q( , ) ( ) ( ) ( ) 1 2 6 3 22 2 2

f k q k kq k q q( , ) 14 12 38 3 18 412 2

Cvičení- příklad 1 II

Derivace:

Řešíme tedy a

Minimum je v bodě k = 0.5 a q = 2

f

kk q 28 12 38

f

qk q 12 6 18

0 28 12 38 k q 0 12 6 18 k q

Cvičení- příklad 2

Metodou nejmenších čtverců najděte přímku, která je nejblíže bodům: (1, 2); (2, 2.3) a (3, 3)

Vytvořte soustavu rovnic A.x = b a vyřešte ji, obecně:

kde k, j {1, …, n}

m

1iijikkj )t().t(a

m

1iijik )t(.yb

Soustava:a11 = 14 a12 = 6 a21 = 6 a22 = 3

b1 = 15,6 b2 = 7,3

Řešení: k = 0,5 q = 1,43

Cvičení- příklad 3

Metodou nejmenších čtverců najděte přímku, která je nejblíže bodům: (-2,0); (0,1) a (2,3)

Využijte vztahy:

m

t.tt

m

y.tt.y

k m

1ii

m

1iim

1i

2i

m

1ii

m

1iim

1iii

m

t.kyq

m

1ii

m

1ii

y ti ii

1

3

6 t ii

1

3

0 yii

1

3

4 t ii

2

1

3

8

k

60 43

80 03

0 75

.

., q

4 0

31 33

34 .

,

Cvičení- příklad 4

Metodou nejmenších čtverců najděte polynom 2.stupně, který je nejblíže bodům:

[1,2], [2,0], [3,3], [4,4].

Řešení: Opět hledáme polynom ve tvaru

Tedy řešíme soustavu rovnic:

y ax bx c 2

t t t

t t t

t t m

a

b

c

t y

t y

y

i i i

i i i

i i

i

i

i

m

i

m

i

m

i

m

i

m

i

m

i

m

i

m

ii

m

ii

m

ii

m

4

1

3

1

2

1

3

1

2

1 1

2

1 1

2

1

1

1

.

Cvičení- příklad 4 II

=> soustava: 354a + 100b + 30c = 93

100a + 30b + 10c = 27

30a + 10b + 4c = 9

Výsledek je:

a = 3/4,

b = -57/20,

c = 15/4.

y t t 3

4

57

20

15

42