Rozwiązywanie algebraicznych układów równań liniowych metodami iteracyjnymi

Rozwiązywanie algebraicznych układów równań liniowych metodami iteracyjnymi

Plan wykładu:

1. Przykłady macierzy rzadkich i formaty ich zapisu 2. Metody: Jacobiego, Gaussa-Seidla, nadrelaksacji 3. Zbieżność metod iteracyjnych

4. Metody: największego spadku, sprzężonego gradientu.

Literatura: Yousef Saad „Iterative methods for sparse linear systems”

DWT 87: WIeża DWT 234: Wieża z platformą DWT 607: Wirnik Wankela Matrix Market - DWT: Everstine's collection from the Harwell-Boeing Collection

Łopata turbiny 2802 x 2802, 6342 NZ

Cylinder z kołnierzem 2919 x 2919,

7593 NZ

3 Wybrane formaty zapisu macierzy rzadkich

- CSR (compressed sparse row) – trzy wektory: wartości, numery kolumn, początki wierszy (pierwsze nie-zero w wierszu)

- CSC (compressed sparse column) – trzy wektory: wartości, numery wierszy, początki kolumn (pierwsze nie-zero w kolumnie)

- CF (coordinate format) – trzy wektory dla: wartości, oraz numery kolumn i wierszy dla nie-zer

CSR dla macierzy symetrycznej – zapamiętujemy tylko macierz U

CSR dla macierzy niesymetrycznej

wartości = (1 -1 -3 5 4 6 4 7 -5) kolumna = (1 2 4 2 3 4 5 4 5) wiersz = (1 4 5 8 9 10)

wartości = (1 -1 -3 -2 5 4 6 4 -4 2 7 8 -5) kolumna = (1 2 4 1 2 3 4 5 1 3 4 2 5) wiersz = (1 4 6 9 12 14)

Mnożenie Ax=y w formacie CSR dla macierzy niesymetrycznej

a[ ] - elementy macierzowe, k[] - numery kolumn, w[]-indeksy z poczatkami wierszy

y=0;

for(i=1; i<n; i++){

l₁=w[i]; // początek indeksow dla i-tego wiersza l₂=w[i+1]-1; // koniec indeksów dla i-tego wiersza

for(l=l₁; l<=l₂; l++){

j=k[l]; //numer kolumny y[i]=y[i]+a[l]*x[j];

} }

X

a

x

= y

5 Mnożenie Ax=y w formacie CSR dla macierzy symetrycznej

A=U+D+L oraz L=U^T

(U+D+L)x = y =(D+U)x + (x

U)

Czyli

Jeśli zamienimy wskaźniki w drugiej sumie to element a_ij*x_i będzie dawać wkład do y_j , które zostanie wyznaczone później (j>i) – w ten sposób dochodząc do wiersza j-tego wartość drugiej sumy będzie znana.

a[ ] - wektor elementów macierzowych, k[] - numery kolumn, w[]-indeksy z początkami wierszy y=0;

for(i=1; i<=n; i++){

l₁=w[i]; // początek indeksow dla i-tego wiersza l₂=w[i+1]-1; // koniec indeksów dla i-tego wiersza

for(l=l₁; l<=l₂; l++){

j=k[l]; //numer kolumny y[i]=y[i]+a[l]*x[j];

if( i!=j ) y[j]=y[j]+a[l]*x[i]; // sumowanie tyko po elementach pozadiagonalnych // modyfikacja pochodzi od drugiej sumy

} }

X

a

x

= y

! X

j=i

a

x

+ X

j<i

x

a

Oznaczmy A jako sumę 3 macierzy

Metoda Jacobiego

Dla dowolnie wybranego przybliżenia rozwiązania x₀ chcemy tak przekształacać iteracyjnie wektor x^(k) aby doprowadzić do znikania składowych wektora reszt w k iteracjach

co można zapisać Szukamy rozwiązania układu n równań

liniowych

Dlaczego używamy metod iteracyjnych?

Przykład

N=50000 – liczba równań w układzie fl₂ = 8 bajtów/liczbę – podwójna precyzja a) Ograniczenia pamięci

P_d<N²fl₂= 20 GB (10GB) – zaalokowana pamięć w komputerze

Ale jeśli układ jest np. pięcioprzekątniowy to do zapisu macierzy A (w postaci wektorowej)

potrzebujemy tylko P_i<5Nfl₂ =2MB pamięci

b) większa wydajność dla macierzy rzadkich (liczba elementów macierzy różnych od 0 jest rzędu N) w stosunku do metod bezpośrednich Macierze takie często pojawiają się w

obliczeniach naukowych i inżynierskich (FEM, PDE)

Ax = b; A 2 R

; x; b 2 R

x = [»

; »

; : : : ; »

] b = [¯

; ¯

; : : : ; ¯

]

(b ¡ Ax

)

= 0

¯

¡ X

j

a

»

= 0 A = L + D + U

D

U

L

7 Składowe wektora reszt znikają w kolejnych

iteracjach, więc możemy zapisać

oraz dla całego wektora

W metodzie Jacobiego obliczamy kolejno wszystkie składowe nowego przybliżenia wektora rozwiązań.

Metoda Gaussa-Seidla

Różni się od metody Jacobiego tym, że obliczone już składniki

wykorzystywane są w obliczeniach składników j+1,j+2,...,n.

»

= 1 a

0

B @¯

¡ X

j j6=i

a

»

1 C A a

»

= ¯

¡

X

j j6=i

a

»

; i = 1; 2; : : : ; n

»

; i = 1; 2; : : : ; j

»

=

= 1 a

0

@¡

i¡1

X

j=1

a

»

¡

X

a

»

+ ¯

1 A

x

= ¡D

(L + U )x

+ D

b

b ¡ Lx

¡ Dx

¡ U x

= 0 x

= ¡D

Lx

¡ D

U x

+ D

b

¯

¡

i¡1

X

j=1

a

»

¡a

»

¡

X

a

»

= 0

8 Metody Jacobiego i GS można zapisać ogólnie

w postaci

M x

= N x

+ b = (M ¡ A)x

+ b A = M ¡ N

metoda Jacobiego:

metoda Gaussa-Seidela:

Metoda relaksacji

Metoda nadrelaksacji (SOR) (Successive Over Relaxation)

M = D

»

= !»

+ (1 ¡ !)»

! 2 (1; 2)

Macierze iterujące i ich przekształcenia (preconditioning)

Ogólny schemat iteracyjny

przy podziale macierzy A

definiujemy iterację do ustalonego punktu w jako

Z porównania obu zapisów dostajemy

A = M ¡ N

x

= M

N x

+ M

b

G = M

N = M

(M ¡ A) = I ¡ M

A f = M

b

x

= Gx

+ f G

(A) = I ¡ D

A

G

(A) = I ¡ (D + L)

A M = D + L

(D + !L)x

= [ ¡!U + (1 ¡ !)D]x

+ !b A = L + D + U

!A = !D + !L + !U

!A = (D + !L) + (!U ¡ (1 ¡ !)D)

b ¡ Lx

¡ Dx

¡ U x

= 0

9 Zbieżność metod iteracyjnych

Dla macierzy

definiujemy liczbę

którą nazywamy promieniem spektralnym macierzy.

Dla dowolnej macierzy kwadratowej zgodnej z normą wektorów prawdziwa jest nierówność

Lemat

Tw. Dla każdego wektora elementy ciągu

dążą do zera wtedy i tylko wtedy gdy

Dowód Proces iteracyjny

możemy potraktować także jako problem rozwiązania układu

co dla G=I-M^-1A daje układ równań

Układ ten ma identyczne rozwiązanie jak układ pierwotny. Co nam to daje?

Z przepisu iteracyjnego

wynika, że musimy w każdej iteracji obliczyć

Ponieważ M-1 nie znamy, więc chcemy

niewielkim kosztem rozwiązać układ równań

Dla metod Jacobiego, Gaussa-Seidla i SOR macierz ta ma postać:

x

= Gx

+ f

(I ¡ G)x = f

M

Ax = M

b

A 2 R

j¸

j · jjAjj; ¸

2 Z

x 2 R

½(A) < 1 Ax; A

x; : : : ; A

x; : : :

^

">0

_

jjAjjp

jjAjj

· ½(A) + "

M

= D

M

= D + L M

= 1

! (D + !L)

½(A) = max

i=1;2;:::;n

j¸

j

x

= M

N x

+ M

b

M

N x

= M

y

= z

M z

= y

" = 1 ¡ ½(A) 2

jjAjj

· 1 + ½(A)

2 < 1

jjA

x jj

· jjAjj

jjxjj

! 0

A

x ! 0

Tw. Ciąg wektorów

którego elementy wyznaczamy według wzoru

jest zbieżny do jedynego punktu granicznego wtedy i tylko wtedy gdy

Dowód

x

; x

; : : : ; x

; : : :

x

= Gx

+ f; i = 0; 1; : : :

x

= Gx

+ f = G(Gx

+ f ) + f = : : :

= G

x

+ (G

f + G

f + : : : + f )

i

lim

G

x

! 0

jjf + G

f + : : : + G

f + : : : jj

·

· X

i=0

jjfjj

jjGjj

= jjfjj

1 ¡ jjGjj

Zbieżność w metodzie SOR

Jeśli macierz układu jest symetryczna,

dodatniookreślona i nieosobliwa to procedura iteracyjna jest zawsze zbieżna dla

G

= (D + !L)

[ ¡!U + (1 ¡ !)D]

det(G

) = det ¡

(D + !L)

¢

£ det (¡!U + (1 ¡ !)D)

det ( ¡!U + (1 ¡ !)D) = det((1 ¡ !)D)

= (1 ¡ !)

det(D) det(G

) = (1 ¡ !)

det(G

) = ¸

¸

: : : ¸

j1 ¡ !j · max

i=1;:::;n

¸

= ½(G

) < 1 0 < ! < 2

0 < ! < 2 det ¡

(D + !L)

¢

= 1

det(D + !L) = 1 det(D)

½(G) < 1

11 Minimalizacja formy kwadratowej

Jeśli Ax=b i r=b-Ax to możemy utworzyć formę kwadratową postaci

Która jest dodatniookreślona i przyjmuje wartość minimalną dla dokładnego rozwiązania x.

W dalszych rozważaniach zakładamy że macierz A jest symetryczna i dodatniookreślona, wówczas możemy użyć formy kwadratowej postaci

która ma minimum w x, ponieważ

Proces poszukiwania rozwiązania dokładnego przebiega iteracyjnie, tj. szukamy ciągu przybliżeń

gdzie:

Od sposobu wyznaczania α_i i v_i zależy zbieżność i szybkość metody.

R = r

r = (b ¡ Ax)

(b ¡ Ax)

Q(x + ¢x) ¡ Q(x) = 1

2 ¢x

A¢x > 0

x

; x

; x

; : : : x

= x

+ ®

v

x₁

x₂

x₆ x₅ x₄ x₃

Związek gradientu Q z kierunkiem poszukiwania przybliżonego rozwiązania. Prosta interpretacja

geometryczna w 2D – powierzchnie o stałej wartości Q mają kształt elipsy (hiperelipsy w przestrzeni o większej liczbie wymiarów).

Q = 1

2 x

Ax ¡ x

b

Metoda największego spadku

Przybliżone rozwiązanie w i+1 iteracji ma postać

Jako v_i wybieramy kierunek gradientu Q

W celu znalezienia współczynnika _i obliczamy Q(x_i+1)

i różniczkujemy je po parametrze

wariacyjnym w celu znalezienia minimum

x

= x

+ ®

v

rQ = Ax

¡ b = ¡r

v

= ¡r

Q(x

¡ ®

r

) = ¡ 1

2 x

r ¡ 1 2 x

b + 1

2 ®

r

Ar

+ ®

r

r

@Q

@®

= r

r

+ ®

r

Ar

@Q

@®

= 0 ! ®

= ¡ r

r

r

Ar

x

= x

+ r

r

r

Ar

r

Kolejne przybliżenie w metodzie największego spadku opisuje wyrażenie

dla którego zachodzi warunek

Metoda może być jednak wolnozbieżna w przypadku gdy hiperelipsoida ma wydłużony kształt co

odpowiada złemu uwarunkowaniu układu.

Q(x

) < Q(x

)

13 Metoda sprzężonego gradientu

Założenia:

- x_d jest rozwiązaniem dokładnym - ciąg wektorów

stanowi bazę w n-wymiarowej przestrzeni euklidesowej

Różnicę rozwiązania dokładnego i przybliżonego możemy zapisać w postaci kombinacji liniowej elementów bazy

Jeśli elementy bazy są ortogonalne to można łatwo wyznaczyć współczynniki kombinacji liniowej

Ale powyższy wzór wymaga modyfikacji ponieważ nie znamy wektora x_d, wiemy jednak, że Ax_d=b więc

v

; v

; v

; : : :

Żądamy więc, aby wektory bazy spełniały

warunek A-ortogonalności (wektory A-sprzężone)

Dla macierzy dodatniookreślonej zachodzi warunek

Jak skonstruować bazę A-ortogonalną?

Jeśli dysponujemy zwykłą bazą wektorów

to możemy ją poddać procesowi ortogonalizacji Grama-Schmidta

Jak utworzyć ciąg wektorów u_i?

v

Av

= 0 $ i 6= j

v

Av

6= 0

u

; u

; u

; : : :

v

= u

v

= u

+ X

¯

v

¯

= ¡ v

Au

v

Av

x

¡ x

=

X

®

v

®

= v

(x

¡ x

)

v

v

; j = 1; 2; : : :

®

= v

A(x

¡ x

)

v

Av

= v

r

v

Av

W metodzie CG bazę stanowią wektory reszt (kierunki gradientów), które dzięki A-

ortogonalizacji są sprzężone.

Kolejne przybliżenia w podstawowej metodzie CG wyznaczamy zgodnie z poniższym schematem:

Dzięki A-ortogonalności w każdej iteracji

wystarczy wyznaczyć tylko jeden współczynnik  (reszta współczynników znika).

W podstawowej metodzie CG w każdej iteracji należy wykonać dwa mnożenia macierz-wektor

i to te dwie operacje determinują nakład obliczeń.

Algorytm metody CG można przedstawić w alternatywnej postaci, gdzie wymagamy tylko jednego mnożenia macierz-wektor:

Maksymalna liczba iteracji w metodzie CG wynosi n+1 – więc jest metodą skończoną. Zazwyczaj do uzyskania akceptowalnego rozwiązania wystarcza wykonanie znacznie mniejszej liczby iteracji.

Av

Ar

v

= r

= b ¡ Ax

¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡

®

= r

r

v

Av

x

= x

+ ®

v

r

= r

¡ ®

Av

¯

= ¡ r

r

r

r

v

= r

+ ¯

v

v

= r

= b ¡ Ax

¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡

®

= v

r

v

Av

x

= x

+ ®

v

r

= r

¡ ®

Av

¯

= ¡ v

Ar

v

Av

v

= r

+ ¯

v

15 W metodzie CG bazę stanowią wektory reszt

(kierunki gradientów), które dzięki A- ortogonalizacji są sprzężone.

Kolejne przybliżenia w podstawowej metodzie CG wyznaczamy zgodnie z poniższym schematem:

Dzięki A-ortogonalności w każdej iteracji

wystarczy wyznaczyć tylko jeden współczynnik  (reszta współczynników znika).

W podstawowej metodzie CG w każdej iteracji należy wykonać dwa mnożenia macierz-wektor

i to te dwie operacje determinują nakład obliczeń.

Algorytm metody CG można przedstawić w alternatywnej postaci, gdzie wymagamy tylko jednego mnożenia macierz-wektor:

Maksymalna liczba iteracji w metodzie CG wynosi n+1 – więc jest metodą skończoną. Zazwyczaj do uzyskania akceptowalnego rozwiązania wystarcza wykonanie znacznie mniejszej liczby iteracji.

Av

Ar

v

= r

= b ¡ Ax

¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡

®

= r

r

v

Av

x

= x

+ ®

v

r

= r

¡ ®

Av

¯

= ¡ r

r

r

r

v

= r

+ ¯

v

v

= r

= b ¡ Ax

¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡

®

= v

r

v

Av

x

= x

+ ®

v

r

= r

¡ ®

Av

¯

= ¡ v

Ar

v

Av

v

= r

+ ¯

v

16 Macierz A_nxn

a

= 1

1 + ji ¡ jj ; ji ¡ jj · 5 _ a

= 0; ji ¡ jj > 5