Pokročilé numerické metody I
Metoda nejmenších čtverců
Jiří Zelinka
Jiří Zelinka Metoda nejmenších čtverců 1 / 16
Metoda nejmenších čtverců
Systém lineárních rovnic
Ax = b
Reziduum:
rx = b − Ax
Definice
Vektor ˆx se nazývá řešením ve smyslu NČ, jestliže
∀x : ∥rx ∥ ≥ ∥rˆx ∥
Věta 1
Řešení ve smyslu NČ vždy existuje a je řešením systému
normálních rovnic
A∗
Ax = A∗
b.
Důkaz: Pomocí ortogonálního rozkladu vektoru b = b1 + b2,
kde b1 ∈ R(A).
Jiří Zelinka Metoda nejmenších čtverců 2 / 16
Důsledek
Systému normálních rovnic je vždy řešitelný.
Věta 2
Každé řešení systému normálních rovnic je řešením ve smyslu
NČ.
Věta 3
Řešení ve smyslu NČ je jediné právě tehdy, když sloupce A
jsou LNZ.
Jiří Zelinka Metoda nejmenších čtverců 3 / 16
Numerický přístup k MNČ
Aproximace dat daným typem funkce
x = (x0, . . . , xN) – uzly
f = (f0, . . . , fN) – funkční hodnoty v uzlech neznáme funkce f
φ0, . . . , φM – zadané bázové funkce (např. φk = xk
)
Data (xi , fi ) aproximujeme funkcí φ =
M
k=0
ckφk
Hledáme minimum výrazu
σ2
(c0, . . . , cM) =
N
j=0
(fj − φ(xj ))2
Vektor c = (c0, . . . , cM)T
je řešením normálního systému
AT
Ax = AT
f pro matici A =



φ0(x0) · · · φM(x0)
...
...
φ0(xN) · · · φM(xN)


.
Jiří Zelinka Metoda nejmenších čtverců 4 / 16
Statistický přístup k MNČ
Lineární regresní model:
Y = β0X0 + . . . βMXM + ε
Y , ε – náhodné veličiny, βj – neznámé parametry.
Pro N pozorování:
Y0 = β0X00 + . . . βMXM0 + ε0
...
YN = β0X0N + . . . βMXMN + εN
t.j.
Y = Xβ + ε
X – matice plánu, β – neznámé parametry
Minimalizujeme výraz (rozptyl odhadu β)
(Y − Xβ)T
(Y − Xβ)
Řešení je odhad parametru β (BLUE)
ˆβ = (XT
X)−1
XT
Y
Jiří Zelinka Metoda nejmenších čtverců 5 / 16
Mooreova-Penroseova pseudoinverze
Definice
Buď A matice typu m × n Matice X se nazývá její
pseudoinverzní maticí, jestliže platí následující vlastnosti
(Penroseovy axiomy)
(1) AXA = A
(2) XAX = X
(3) (AX)∗
= AX (AX je hermitovská)
(4) (XA)∗
= XA (XA je hermitovská)
Věta 4
Pseudoinevrzní matice vždy existuje a je určena jednoznačně.
Označení
Pseudoinverzní matici k matici A značíme A+
.
Jiří Zelinka Metoda nejmenších čtverců 6 / 16
Lemma (skeletní rozklad)
Nechť A je typu m × n, r(A) = r. Pak existují matice B typu
m × r a C typu r × n, r(B) = r(C) = r, že platí A = BC.
Důkaz: Jako B vybereme nezávislé sloupce matice A, C pak
dopočítáme.
Důkaz vety o existenci a jednoznačnosti
Jednoznačnost dostaneme z Penroseových axiomů.
Pro nulovou matici A je A+
= A∗
.
Pokud má A nezávislé sloupce, pak A+
= (A∗
A)−1
A∗
,
pro matici s nezávislými řádky je A+
= A∗
(AA∗
)−1
.
V obecném případě uděláme skeletní rozklad A = BC a
dokážeme, že A+
= C+
B+
.
Jiří Zelinka Metoda nejmenších čtverců 7 / 16
Věta 5 Vlastnosti pseudoinverzní matice
1 pro regulární A je A+
= A−1
2 (A+
)+
= A
3 (A∗
)+
= (A+
)∗
4 (A∗
A)+
= A+
A+∗
, (AA∗
)+
= A+∗
A+
5 A+
= (A∗
A)+
A∗
= A∗
(AA∗
)+
6 h(A) = h(A+
) = h(A∗
) = h(A+
A) = h(A∗
A) = . . .
7 R(A+
) = R(A∗
)
8 N(A+
) = N(A∗
)
9 obecně neplatí (AB)+
= B+
A+
10 c ̸= 0 ⇒ (cA)+
= 1
c
A+
11 A+
= (QA∗
A)+
QA∗
pro regulární Q
Jiří Zelinka Metoda nejmenších čtverců 8 / 16
M–P pseudoinverze a MNČ
Lemma
Buď P = AA+
, Q = A+
A. Pro libovolné vektory x, y platí
∥Ax + (I − P)y∥2
= ∥Ax∥2
+ ∥(I − P)y∥2
∥A+
y + (I − Q)x∥2
= ∥A+
y∥2
+ ∥(I − Q)x∥2
Důkaz: dokážeme kolmost.
Věta 6
Buď A typu m × n, b typu n × 1, ˆx = A+
b. Pak pro každé
x ̸= ˆx platí: ∥b − Ax∥ ≥ ∥b − Aˆx∥, a v případě, že
∥b − Ax∥ = ∥b − Aˆx∥, pak ∥x∥ > ∥ˆx∥.
Interpretace:
ˆx = A+
b je řešení ve smyslu NČ s minimální normou a je
jediné s touto vlastností.
Jiří Zelinka Metoda nejmenších čtverců 9 / 16
Důkaz
Ax − b = Ax − AA+b + AA+b − b = A(x − A+b) + (I − AA+)(−b)
Podle předchozího lemmatu
∥Ax − b∥2
= ∥A(x − ˆx)∥2
+ ∥Aˆx − b∥2
≥ ∥Aˆx − b∥2
,
přičemž rovnost nastane pouze pro
∥A(x − ˆx)∥2
= 0, tj. Ax = Aˆx.
Pokud Ax = Aˆx, pak A+
Ax = A+
Aˆx = A+
AA+
b = A+
b = ˆx.
Potom x = ˆx + (x − ˆx) = A+
b + (I − A+
A)x.
Podle druhé rovnosti předchozího lemmatu platí
∥x∥2
= ∥ˆx∥2
+ ∥x − ˆx∥2
, takže ∥x∥2
≥ ∥ˆx∥2
a rovnost
nastane pouze pro x = ˆx.
Jiří Zelinka Metoda nejmenších čtverců 10 / 16
M–P pseudoinverze a systém lineáních rovnic
Věta 7 (Penrose 1955)
Systém rovnice AXB = C je řešitelný právě tehdy, když platí
AA+
CB+
B = C, V tom případě je obecné řešení systému ve
tvaru
X = A+
CB+
+ Y − A+
AYBB+
pro libovolnou matici Y .
Důsledek
Systém Ax = b je řešitelný právě tehdy, když platí AA+
b = b.
Obecné řešení tohoto systému je pak
x = A+
b + (I − A+
A)y
pro libovolný vektor y.
Jiří Zelinka Metoda nejmenších čtverců 11 / 16
Grevilleův algoritmus
Výpočet Mooreovy–Penroseovy pseudoinverze po
řádcích (sloupcích)
A: matice, Ak: prvních k řádků, ak: k-tý řádek, známe A+
k−1
Položme:
dk = akA+
k−1
ck = ak − dkAk−1
bk =



c+
k ck ̸= o
1
1+dk d∗
k
A+
k−1d∗
k ck = o
A+
k = A+
k−1 − bkdk|bk
Jiří Zelinka Metoda nejmenších čtverců 12 / 16
Sloupcová verze
A: matice, Ak: prvních k sloupců, ak: k-tý sloupec, známe
A+
k−1
Položme:
dk = A+
k−1ak
ck = ak − Ak−1dk
bk =



c+
k ck ̸= o
1
1+dk d∗
k
d∗
k A+
k−1 ck = o
A+
k =


A+
k−1 − dkbk
− − − − −−
bk


Využití: MNČ – přidání pozorování nebo vysvětlující proměnné
Jiří Zelinka Metoda nejmenších čtverců 13 / 16
Nelineární MNČ
Aproximace dat nelineární funkcí závislé na parametrech, které
se snažíme určit
Příklad:
Michaelisův-Mentenové model kinetiky enzymů
f (x, β) =
β1x
β2 + x
x0, . . . , xN – uzly
y0, . . . , yN – funkční hodnoty, yj ≈ f (xj , β)
Rezidua: rj (β) = yj − f (xj , β) Hledáme minimum výrazu
S(β) =
N
j=0
rj (β)2
Jiří Zelinka Metoda nejmenších čtverců 14 / 16
Nutná podmínka pro minimum:
0 =
∂S
∂βi
= 2
N
j=0
rj (β)
∂rj (β)
∂βi
= −2
N
j=0
Jji rj (β) ∀i,
pro
∂rj (β)
∂βi
= −
∂f (xj ,β)
∂βi
= −Jji ,
J = (Jji ) je Jacobiova matice, J = J(β).
Máme systém nelineárních rovnic – řešíme iterační metodou:
β0
- počáteční iterace, βk
+ ∆βk
= βk+1
≈ β
Taylorův rozvoj v parametrech β:
f (xj , β) ≈ f (xj , βk
) +
i
∂f (xj ,βk
)
∂βi
(βi − βk
i ) ≈
≈ f (xj , βk
) +
i
Jji ∆βk
i .
Jiří Zelinka Metoda nejmenších čtverců 15 / 16
Označení: ∆yk
j = yj − f (xj , βk
)
rj (β) = yj −f (xj , β) = (yj −f (xj , βk
))+(f (xj , βk
)−f (xj , β)) ≈
≈ ∆yk
j −
i
Jji ∆βk
i = ∆yk
j −
l
Jjl ∆βk
l
Dosazením do nutné podmínky pro minimum dostaneme
−2
j
Jji ∆yk
j −
l
Jjl ∆βk
l = 0, ∀i
Po úpravě dostaneme systém normálních rovnic
JT
J∆βk
= JT
∆yk
, pro J = J(βk
).
Gaussův-Newtonův algoritmus: volíme počáteční iteraci β0
,
z předchozího systému normálních rovnic určíme přírůstek
∆β0
, z něj iteraci β1
atd.
Jiří Zelinka Metoda nejmenších čtverců 16 / 16