Metody rozwiązywania układów równań
liniowych
n
n
nn
n
n
n
n
n
n
b
x
a
x
a
x
a
b
x
a
x
a
x
a
b
x
a
x
a
x
a
2
2
1
1
2
2
2
22
1
21
1
1
2
12
1
11
b
Ax
n
nn
n
n
n
n
b
b
b
a
a
a
a
a
a
a
a
a
2
1
2
1
2
22
21
1
12
11
b
A
0
det
A
Metody skończone:
• Metoda Gaussa
• Metoda Gaussa-Jordana
• Metody Choleskiego
• Metoda Householdera
• Metoda sprzężonych gradientów
Metody iteracyjne dla dużych układów równań:
• Metoda Jacobiego
• Metoda Gaussa-Seidla
n
n
nn
n
n
n
n
c
x
r
c
x
r
x
r
c
x
r
x
r
x
r
2
2
2
22
1
1
2
12
1
11
Metoda eliminacji Gaussa z wyborem elementu głównego w
kolumnie
Układ równań sprowadzamy do postaci trójkątnej
Układ z macierzą trójkątną można następnie łatwo rozwiązać
zaczynając od obliczenia wartości x
n
z n-tego równania,
następnie wstawić x
n
do równania n-1 i wyliczyć z niego x
n-1
,
następnie wstawić x
n
oraz x
n-1
do równania n-2 i wyliczyć x
n-2
aż
do dotarcia do równania pierwszego i wyznaczenia x
1
.
1. Wybieramy równanie i takie, że |a
i1
| jest największym elementem w
pierwszej kolumnie po czym przestawiamy i-te równanie na początek i
eliminujemy x
1
z równań od 2 do n.
)
1
(
)
1
(
2
2
)
1
(
2
)
1
(
2
)
1
(
2
2
)
1
(
22
)
0
(
1
)
0
(
1
2
)
0
(
12
1
)
0
(
11
n
n
n
n
n
n
n
n
b
x
a
x
a
b
x
a
x
a
b
x
a
x
a
x
a
1
1
,
1
)
0
(
11
)
0
(
1
)
0
(
1
)
0
(
)
1
(
)
0
(
11
)
0
(
1
)
0
(
1
)
0
(
)
1
(
i
a
a
b
b
b
k
i
a
a
a
a
a
i
i
i
i
k
ik
ik
2. Procedurę powtarzamy z macierzą A
(1)
o rozmiarach (n-1)x(n-1) i
wektorem b
(1)
o rozmiarze n-1, eliminując z nich drugą zmienną i
otrzymując macierz A
(2)
o rozmiarach (n-2)x(n-2) i wektor b
(2)
o
rozmiarze n-2. W ten sam sposób postępujemy z kolejnymi macierzami
A
(2)
, A
(3)
,..., A
(n-1)
oraz wektorami b
(2)
, b
(3)
,..., b
(n-1)
.
j
i
a
a
b
b
b
j
k
j
i
a
a
a
a
a
j
jj
j
ij
j
j
j
i
j
i
j
jj
j
ij
j
jk
j
ik
j
ik
)
1
(
)
1
(
)
1
(
)
1
(
)
(
)
1
(
)
1
(
)
1
(
)
1
(
)
(
,
Dla j-tego kroku
Po zakończeniu operacji otrzymujemy układ równań z macierzą
trójkątną
)
1
(
)
1
(
22
)
0
(
11
)
1
(
)
1
(
)
2
(
3
)
0
(
3
3
)
2
(
33
)
1
(
2
)
0
(
2
3
)
1
(
23
2
)
1
(
22
)
0
(
1
)
0
(
1
3
)
0
(
13
2
)
0
(
12
1
)
0
(
11
)
1
(
det
n
nn
p
n
n
n
n
nn
n
n
n
n
n
n
a
a
a
b
x
a
b
x
a
x
a
b
x
a
x
a
x
a
b
x
a
x
a
x
a
x
a
A
p jest liczbą przestawień wierszy macierzy A podczas
sprowadzania układu równań do postaci trójkątnej.
1
,...,
2
,
1
1
)
1
(
)
1
(
)
1
(
)
1
(
)
1
(
)
1
(
n
n
j
x
a
a
a
b
x
a
b
x
k
n
j
k
j
jj
j
jk
j
jj
j
j
j
n
nn
n
n
n
3. Z otrzymanego układu równań z macierzą trójkątną
wyznaczamy po kolei x
n
, x
n-1
,..., x
1
.
Wysiłek obliczeniowy (liczba mnożeń i dzieleń) w metodzie
eliminacji Gaussa:
Faktoryzacja macierzy A: n(n
2
-1)/3 operacji
Przekształcenie wektora b: n(n-1)/2 operacji
Obliczenie x: n(n+1)/2 operacji.
Razem: n
3
/3+n
2
-n/3≈n
3
/3 operacji.
Metody typu Choleskiego dla macierzy symetrycznych silnie
nieosobliwych
T
LDL
A
L
T
L
D
L
T
LL
A
klasyczna metoda Choleskiego
tylko dla macierzy dodatnio
określonych.
Postępowanie przy rozwiązywaniu układów równań metodą
faktoryzacji Choleskiego.
1. Wyznaczenie faktorów L i D. Układ przyjmuje postać
LDL
T
x=b
2. Obliczenie pomocniczego wektora w.
w=L
-1
b przez rozwiązanie układu równań Lw=b.
Ponieważ L jest macierzą trójkątną dolną układ ten rozwiązuje się
wyliczając kolejno w
1
, w
2
,…, w
n
podobnie jak w koncowym etapie
eliminacji Gaussa.
3. Obliczenie z=D
-1
w (D jest macierzą diagonalną więc po prostu
dzielimy w
i
przez d
ii
. Ten etap nie występuje w klasycznej metodzie
Choleskiego.
4. Obliczenie x poprzez rozwiązanie układu równań z macierzą
trójkątną górną
L
T
x=z
Ten etap jest identyczny z ostatnim etapem metody eliminacji
Gaussa.
Metoda wymaga ok. n
3
/6 operacji (2 razy mniej niż metoda eliminacji
Gaussa). Uwaga: klasyczna metoda Choleskiego wymaga ponadto
n pierwiastkowań.
Klasyczna faktoryzacja Choleskiego (A=LL
T
)
1
1
1
1
2
11
1
1
11
11
1
i
k
jk
ik
ij
ii
ji
i
k
jk
ii
ii
j
j
l
l
a
l
l
l
a
l
l
a
l
a
l
Faktoryzacja “bezpierwiastkowa”
1
0
1
0
0
1
1
,
1
21
2
1
n
n
n
n
l
l
l
d
d
d
L
D
i
n
k
ik
jk
k
ji
ji
j
j
n
k
nk
k
nn
n
i
k
ki
k
ii
i
d
l
l
d
a
l
d
a
l
l
d
a
d
l
d
a
d
a
d
1
1
1
1
1
1
1
2
1
1
2
11
1
/
Regresja nieliniowa
)
,
,
,
;
,
,
(
2
1
2
1
2
1
2
2
22
21
1
1
12
11
m
m
n
nm
n
n
m
m
p
p
p
x
x
x
f
y
y
x
x
x
y
x
x
x
y
x
x
x
f jest funkcją nieliniową względem parametrów p
1
, p
2
,
…, p
m
.
Przykład problemu nieliniowego linearyzowalnego: kinetyka reakcji
pierwszego rzędu
o
o
o
k
C
k
p
C
y
t
x
kt
C
t
A
kt
C
t
A
B
A
ln
ln
ln
ln
exp
3
2
1
3
2
1
2
2
1
2
1
1
3
2
1
/
exp
exp
1
3
2
1
k
k
k
C
C
y
t
x
t
k
k
k
k
t
k
k
k
k
k
C
t
C
C
B
k
k
k
B
A
C
A
o
o
k
k
k
p
Przykład problemu nieliniowego nielinearyzowalngo:
kinetyka reakcji pierwszego rzędu z produktem
przejściowym
2
2
2
2
1
T
1
2
2
1
2
1
2
1
T
2
2
1
2
1
1
0
0
0
0
1
0
0
0
1
)
;
(
)
;
(
,
,
,
;
,...,
,
1
)
;
(
)
;
(
,
,
,
;
,...,
,
n
n
i
m
im
i
i
i
i
n
i
m
im
i
i
i
p
p
p
x
x
x
f
y
p
p
p
x
x
x
f
y
W
p
X
f
y
W
p
X
f
y
p
X
f
y
p
X
f
y
regresja
nieważona
regresja ważona
Metoda Newtona-Gaussa
Rozwijamy funkcję f dla każdego punktu pomiarowego
w otoczeniu arbitralnego przybliżenia parametrów p
0
w
szereg Taylora
m
j
j
ij
i
m
j
j
j
p
p
p
x
x
x
i
m
in
i
i
m
in
i
i
p
J
f
p
p
p
f
p
p
p
x
x
x
f
p
p
p
x
x
x
f
m
in
i
i
1
)
0
(
)
0
(
)
0
(
1
0
,
,
,
,
,
,
0
0
2
0
1
2
1
2
1
2
1
0
0
2
0
1
2
1
,
,
,
;
,
,
,
,
,
,
;
,
,
,
p
)
0
(
)
0
(
)
0
(
1
)
0
(
)
0
(
)
0
(
)
0
(
p
J
r
p
m
j
j
ij
i
i
i
p
J
f
y
r
W ten sposób dostajemy funkcję liniową względem
przyrostów parametrów p
(0)
. Macierz J nazywa się macierzą
Jacobiego zagadnienia a wektor r wektorem reziduów.
n
i
m
j
j
ij
i
n
i
m
j
j
ij
p
J
r
w
p
J
r
i
i
1
)
0
(
)
0
(
)
0
(
T
)
0
(
)
0
(
)
0
(
2
1
)
0
(
)
0
(
)
0
(
)
0
(
1
)
0
(
)
0
(
)
0
(
T
)
0
(
)
0
(
)
0
(
2
1
)
0
(
)
0
(
)
0
(
)
0
(
p
J
r
W
p
J
r
p
J
r
p
J
r
regresja
nieważona
regresja ważona
(0)
(0)T
(0)
(0)
(0)T
(0)
(0)T
(0)
(0)
(0)T
Wr
J
p
WJ
J
r
J
p
J
J
Algorytm Newtona-Gaussa
1. Wybieramy przybliżenie parametrów p
(0)
.
2. Obliczamy wektor reziduów oraz sumę kwadratów odchyleń
(0)
.
3. Obliczamy macierz Jacobiego a następnie macierz i wektor
wyrazów wolnych układu równań.
4. Obliczamy wektor przyrostów parametrów p
(0)
.
5. Obliczamy p
(1)
=p
(0)
+p
(0)
a następnie nową sumę kwadratów
odchyleń.
6. Jeżeli przyrosty parametrów są odpowiednio małe, zmiana
sumy kwadratów odchyleń jest odpowiednio mała lub
przekroczono dopuszczalną liczbę iteracji, kończymy procedurę.
Jeżeli nie, przechodzimy do punktu 3 wstawiając p
(1)
za p
(0)
i
(1)
za
(0)
.
Zbieżność metody Gaussa jest na ogół rzędu pierwszego. Tylko
jeżeli rezidua dla optymalnego rozwiązania są zerowe zbieżność
jest kwadratowa. Szybkość zbieżności zależy bardzo silnie od
wielkości końcowych reziduów.
Metoda Newtona-Gaussa może być niestabilna szczególnie jeżeli
rezidua odpowiadające początkowym przybliżeniom parametrów
są duże i wskaźnik uwarunkowania macierzy J
T
J jest duży. Można
ją jednak poprawić stosując skracanie przyrostów parametrów.
Szukamy takiego , że (p
(0)
+p
(0)
)<(p
(0)
). Teoretycznie takie
zawsze istnieje, ponieważ kierunek metody Newtona-Gaussa jest
kierunkiem poprawy.
0
2
2
2
)
0
(
2
/
1
)
0
(
1
0
0
0
)
0
(
)
0
(
r
J
J
J
r
J
J
J
J
r
p
J
p
J
r
p
J
r
p
J
r
p
p
T
T
T
T
(0)
(0)
(0)
(0)
T
T
T
d
d
d
d
Wyznaczanie :
Metoda Hartleya: wyznaczamy tak, żeby zminimalizować jako
funkcję a (minimalizacja kierunkowa) Metoda ta zwykle nie jest
stosowana w takiej wersji ze względu na czasochłonność; zwykle
kończy się jak tylko
(1)
<
(0)
.
Proste skracanie kroku: Zakładamy, że =2
-k
gdzie k jest liczbą
całkowitą dobraną tak, że
F(1)<F(0).
Jeżeli
macierz J
T
J jest źle uwarunkowana to zarówno metoda Hartleya jak i
metoda skracania kroku działają bardzo wolno. Metodą sprawdzoną dla takich
przypadków jest metoda Levenberga-Marquardta.
Metoda Levenberga-Marquardta (metoda Marquardta)
Zamiast układu równań normalnych z metody Newtona-Gaussa
rozwiązujemy układ następujący:
r
J
p
D
J
J
T
T
Macierz D jest dodatnią macierzą diagonalną; zwykle D=I lub
D=diag{(J
T
J)
11
,…,(J
T
J)
mm
}.
dobiera się tak aby uzyskać zmniejszenie sumy kwadratów
odchyleń w stosunku do poprzedniej iteracji. Dowód, że metoda jest
metodą kierunku poprawy przebiega jak w przypadku metody
Newtona-Gaussa ze skracaniem kroku.
W pierwotniej wersji Levenberga a było dobierane tak aby uzyskać
minimum w funkcji (kosztowne).
Metoda Marquardta:
1. Na początku przyjmujemy
0
.
2. W każdej kolejnej iteracji przyjmujemy najpierw
(p)
=
(p-1)
/V.
3. Obliczamy
(p)
().
4. Jeżeli
(p)
()<
(p-1)
przechodzimy do następnej iteracji z nową
wartością .
5. Jeżeli
(p)
()>=
(p-1)
i nie przekroczono maksymalnej
dopuszczalnej wartości liczby krokow modyfikacji wstawiamy
:=*V i przechodzimy do punktu 3. Jeżeli przekroczono
maksymalną liczbę kroków kończymy optymalizację.
Metoda Marquardta jest jedną z “najodporniejszych” metod
minimalizacji sum kwadratów. Niepowodzenie optymalizacji jest
na ogół sygnałem błędu w programie (najczęściej błędnemu
obliczeniu pochodnych).
n
i
i
i
r
f
y
m
n
1
2
2
1
p
Macierz wariancji-kowariancji parametrów:
Wariancja resztowa:
1
2
T
*
*
J
J
p
p
p
p
T
r
E
1
2
T
*
*
WJ
J
p
p
p
p
T
r
E
Odchylenia standardowe:
ii
r
p
i
1
WJ
J
T
Metoda Newtona w rozwiązywaniu problemu regresji
nieliniowej
n
i
m
i
i
i
m
n
i
m
i
i
m
p
p
p
f
y
p
p
p
p
p
p
f
y
p
p
p
1
2
2
1
2
2
1
1
2
2
1
2
1
,
,
,
;
1
,
,
,
,
,
,
;
,
,
,
x
x
Warunek konieczny minimum
0
p
Wr
p
J
p
0
p
r
p
J
p
p
p
x
*
*
2
*
*
*
2
*
,
,
2
,
1
0
*
*
2
,
,
,
;
2
1
*
*
2
*
1
1
,
,
,
,
,
,
*
*
2
*
1
*
*
2
*
1
T
T
n
k
k
ki
m
i
k
n
k
p
p
p
i
k
p
p
p
i
m
i
r
J
p
p
p
f
y
p
f
p
m
m
regresja nieważona
regresja ważona
regresja nieważona
regresja ważona
Otrzymujemy zatem układ m równań nieliniowych z m
niewiadomymi, który można rozwiązać metodą Newtona.
r
J
Δp
H
r
J
Δp
r
:
f
J
J
T
pp
T
pp
T
m
i
r
p
f
p
r
p
p
f
p
f
p
f
m
i
p
p
p
p
p
p
n
k
k
i
k
m
j
j
n
k
n
k
k
j
i
k
j
k
i
k
m
j
j
j
p
p
p
j
i
p
p
p
i
p
p
p
i
m
m
m
,
,
2
,
1
,
,
2
,
1
0
1
1
1
1
2
1
*
,
,
,
2
,
,
,
,
,
,
2
1
2
1
*
*
2
*
1
H – hesjan (Hessian) funkcji . Macierz J
T
J: dodatnio określona
część hesjanu (występuje w metodzie Newtona-Gaussa.
Algorytm Newtona dla zagadnienia regresji nieliniowej:
1. Wybrać początkowe przybliżenia parametrów p(0).
2. Dla k-tej iteracji obliczyć rezidua, funkcję minimalizowaną, macierz
Jacobiego oraz drugie pochodne wielkości mierzonych względem
parametrów.
3. Rozwiązać układ równań H
(k)
p
(k)
=J
T(k)
r
(k)
4. Jeżeli ||p
(k)
||<e zakończyć proces iteracyjny, w przeciwnym
przypadku wstawić p
(k+1)
=p
(k)
+p
(k)
i przejść do punktu 2.
Metoda Newtona jest zbieżna kwadratowo, jeżeli hesjan w minimum
jest nieosobliwy.
Metoda Newtona różni się od metody Newtona-Gaussa tym, że
linearyzuje się tutaj gradient funkcji minimalizowanej a nie od razu
teoretyczną zależność pomiędzy zmiennymi objaśniającymi i
niezależnymi. Metoda Newtona-Gaussa „gubi” drugie pochodne
tych wielkości względem parametrów.