Zastosowanie teorii grafów w metodach numerycznych**

Maciej M. SYSŁO {Wrocław)

Stosowana teoria grafów*

Zastosowanie teorii grafów w metodach numerycznych**

Wstęp. Praca ta poświęcona jest zastosowaniom teorii graf ów w metodach numerycznych.

Nie jest to kompletne opracowanie tego tematu, dalsze szczegółowe rozważania można zna-

leźć w pracach podanych w bibliografii, przedstawione zagadnienia reprezentują jednak naj-

ważniejsze działy metod numerycznych, które ^stosują aparat teorii grafów.

Charakter zastosowań teorii graf ów wynika z charakteru samej teorii. Graf p~zyporządko­

wany ^jakiemuś obiektowi lub postępowaniu odzwierciedla bowiem jego ^strukturę topolo- giczną, a dodatkowe opisanie funkcji na wierzchołkach lub/i ^łukach tego grafu pozwala na otrzymanie rezultatów liczbowych ^zależnych przede wszystkim od struktury grafu, a co za tym idzie i od struktury rozpatrywanego obiektu lub postępowania. Tak ^więc w metodach numerycznych, teoria grafów znajduje zastosowanie przy organizacji obliczeń (rozdz. 2 i 3) oraz w metodach przekształcania struktury qbiektów ^obliczeń numerycznych (zwykle ma- cierzy) do postaci dogodnych dla ^obliczeń najbardziej ekonomicznych (rozdz. 3 i 4).

Numeracja wzorów jest oddzielna dla każdego rozdziału, rysunki natomiast numerowa- ne ^są kolejno w ^całej pracy.

Zamieszczona na ^końcu pracy bibliografia zawiera ^także niecytowane prace poświęcone

przedstawionym zagadnieniom.

1. Pojęcia wstępne. W tym rozdziale zdefiniujemy pojęcia teorii grafów używane w całej

pracy, inne natomiast ^(związane nie tylko z grafami) ^używane tylko w poszczególny_ćh roz-

działach, zdefiniowane ^są w pierwszych paragrafach tych rozdziałów.

Grafem zorientowanym (lub.digrafem) nazywamy parę uporządkowaną D = (X, R ),

gdzie X jest skończonym zbiorem elementów zwanych wierzchołkami a R ^{jest relacją} dwu-

argumentową RC X x X. Oznaczmy przez U = { [x, y]: x, y e X/\ xRy} zbiór par uporząd­

kowanych zwanych łukami. Graf można także określić jako parę uporządkowaną D = (X, U>.

Dalej używać będziemy tych ^oznaczeń zamiennie. ^Łuk postaci [x, x] nazywamy pętlą. Jeśli

R ^jest relacją symetryczną, to D .nazywamy grafem symetrycznym i zamiast dwóch łuków [x, y] i [y, x] wprowadzamy jedną parę nieuporządkowaną (x, y) zwaną krawędzią. Ciąg

*Pod tym ^tytułem, autor ma zamiar przedstawić kilka prac poświęconych zastosowaniom teorii graf ów w innych ^działach matematyki, w technice, ekonomii i innych dziedzinach nauki. Prace te w za- sadzie będą przeglądowe z obszerną "r>ibliografią. Jednocześnie zawierać będą pewne wyniki prac autora.

**Jest to wykład wygłoszony na IV Letniej Szkole Metod Numerycznych i Infwmatyki, która od- była się w ^Gdańsku w dniach 20.08-1.09.1973.

[69)

wierzchołków grafu [x ^I> x 2 , • •• , x _1] taki, że x/ix;+ ₁ (i= 1, 2, ... , I - 1) rtazywamy drogą;

a jeśli x ₁ = x _1, to ciąg ten nazywamy konturem. Dla grafów symetrycznych kontur ^nazywać

będziemy cyklem. Podgrafem grafu D nazywamy graf D' = _{-+ -+} <x', R'>. ^taki, _- ^{że X'} _- ^{C X oraz}

x, y e x' /\ xRy ^--+ xRy. Grafem częściowym grafu D nazywamy graf D = <X, U> taki, że

fJ <: U. Graf nazywamy silnie spójnym, jeśli dla każdej uporządkowanej pary wierzchołków istnieje w grafie łącząca je droga. Składową silnej spójności grafu D nazywamy maksymalny, silnie spójny podgraf D' grafu D. W przypadku grafów symetrycznych mówić będziemy o ^spójności i składowych spójności. Grafem obciążonym lub ^siecią nazywamy ^trójkę upo-

rządkowaną D1 =<X, U;f>, gdzie D =<X, U> jest grafem, af funkcją rzeczywistą opisaną na podzbiorze wierzchołków lub ^łuków.

Moc zbioru X oznaczać będziemy przez IXI. Dalsze pojęcia i fakty z teorii grafów zna-

leźć można w książkach [3], [9], [21] i [29].

2. Analiza wyrażeń arytmetycznych

2.0. Drzewem (tree) nazywamy graf symetryczny, spójny i nie zawierający cykli, oznacz- my go przez T = (X, V>. ^Jeśli w definicji tej pominiemy spójność, to otrzymamy określenie

lasu (forest). Wyróżnijmy dowolny wierzchołek drzewa, tzw. korzeń (root). W dalszej części tego rozdziału przez drzewo rozumieć będziemy drzewo z ustalonym korzeniem. Za stopień wierzchołka. danego drzewa przyjmujemy liczbę krawędzi incydentnych z tym ^wierzchoł­

kiem pomniejszoną o jeden. Wierzchołek nazywamy końcowym, jeśli jego stopień jest równy zero, pozostałe natomiast nazywamy pośrednimi lub_rozgałęzieniami. Jeśli stopień każdego wierzchołka pośredniego jest równy dwa, to drzewo nazywamy binarnym. Zauważmy, że

liczba wierzchołków pośrednich drzewa binarnego jest o jeden mniejsza od liczby jego wierz-

chołków końcowych.

Oznaczmy przez I; ^długość drogi (liczbę łuków) z korzenia do wierzchołka i (pośrednie- go lub końcowego). Dla danego drzewa z ustalonym korzeniem liczby te określone są jedno- znacznie. Niech każdemu z N wierzchołków końcowych przyporządkowana będzie liczba dodatnia wi' tzw. waga wierzchołka końcowego. Długością drzewa binarnego T nazywamy

sumę

(1) ITI= L ^{N liwi"}

j= 1

Jednym z szerokich zastosowań drzew jest analiza wyrażeń arytmetycznych (patrz [ 12], str. 362-405). Każdemu bowiem wyrażeniu arytmetycznemu przyporząd~ować można drzewo, którego wierzchołkami końcowymi są argumenty, a pośrednimi - działania arytme- tyczne. Dla przykładu, na rysunku 1 pokazane jest drzewo odpowiadające wyrażeniu aryt- metycznemu (a* b )tc/(-d). Zauważmy, że jeśli wyrażenie zawiera tylko działania dwuargu- mentowe, to odpowiada mu drzewo binarne.

2.1. W tym paragrafie zakładać-będziemy, że (a) błąd względny zaokrąglenia liczby

zmiennopozycyjnej x = m · 1 oc o mantysie m z t cyframi dziesiętnymi po przecinku nie

przekracza 5 . 1 o-t oraz (b) zmiennopozycyjny wynik dodawania liczb zmiennopozycyjnych

ma błąd względny zaokrąglenia nie większy niż 5.10-t.

Rys. 1

P r o b 1 e m: Określić, jak należy obliczać wartość sumy N

W= L ^W;,

j= 1

gdzie w;> ^{O (j} = 1, ... , N), aby oszacowanie błędu zaokrągleń a było najmniejsze. Rozwią­

zanie tego problemu zostało podane w pracy [20].

W rozważaniach naszych ograniczymy się do takich sposobów (algorytmów) obliczania sumy W, które polegają11a wykonaniu tylko n - 1 ^dodawań, który<;h argiimentami mogą być

liczby w 1 , ••• , w ^N lub ich sumy częściowe.

Przy obliczaniu wartości sumy W ^metodą „od lewej do prawej" błąd zaokrągleń nie prze- kracza wielkości

a ₁ = 5 · 10- t c~: (N - N j + 1) W; - w _{1) ,} j= 1

która przyjmuje wartość minimalną dla liczb w 1 , w 2 , ••• , w ^N uporządkowanych od najmniej- szej do największej. Tej metodzie obliczania wartości sumy W odpowiada drzewo binarne T _{1 ,}

. N

przedstawione na rysunku 2, którego długość wynosi IT ₁ I= L ^{(N - j + 1) w; - w} 1 , a więc

a ₁ = ^{5 · 10- t I T} ₁ I. i= 1

Łatwo zauważyć, że każdemu sposobowi obliczania wartości sumy W odpowiada drzewo binarne To N wierzchołkach końcowych obciążonych wagami w _{1 ,} w 2 , ••• , w ^N i N - 1

wierzchołkach pośrednich odpowiadających operacjom dod;iwania, i na odwrót, oraz, że błąd zaokrągleń sumy W zgodnie z drzewem T nie przekracza a = ^{5 · 1} o- ^{t I Tl.}

Tak więc minimalizacja oszacowania błędu zaokrągleń a jest równoważna minimalizacji sumy (1) po wszystkich d·opuszczalnych ciągach liczb naturalnych I= (/ 1 , / 2 , „„ /N),.tj. od-

powiadających drzewom binarnym. Łatwo można wykazać, że ciąg I liczb naturalnych jest

dopuszczalny wtedy i tylko wtedy, jeśli spełnia następujące warunki: ·

Rys. 2

(i) Niech m = max Ir ^Ciąg I zawiera parzystą liczbę elementów równych m, które wy-

stępują w odcinkach parzystej długości.

(ii) Ciąg zredukowany !', utworzony przez zastąpienie (od lewej do prawej) kolejnych dwóch elementów ciągu I równych m przez m - 1, spełnia warunek (i).

(iii) Powtarzając proces (ii) dla kolejno otrzymywanych ciągów zredukowanych speł­

niających (i) otrzymujemy ostatecznie ciąg zredukowany 1,1.

Ponumerujmy wierzchołki pośrednie drzewa binarnego w sposób jednoznaczny. Niech bi będzie wartością sumy częściowej obliczaną w j-tym wierzchołku pośrednim (j = 1, ... , N - 1 ). Łatwo wykazać, że

(2) ^N-1 _L ~ b. =~I.w _{I L I I} ^N ..

j= 1 j= 1

Rozwiązaniem postawionego problemu jest drzewo binarne o najmniejszej długości,

tzw. drzewo Huffmana ([ 11 ], [ 12]). Odpowiada mu następująca kolejność obliczania sumy W (algorytm Huffmana):

O) { _v ^{o o} _{1 , v 2 ,} „ „ v ^o Nf ^ł +- { w _{1 ,} ^„ „ w N , } g ^{d . ~1e} o < _v ^o _{1 ....,,} ^~ v ^{o ..,....} 2 -:::::: „ . ....,, ^~ v ^o N' k +-o.

1) Wykonać N - 1 razy następujące czynności:

{V~+ 1 ' V~+ 1 , ... ' vt ^{~ ~} - ^{1} +-l} V~ + V~, V~, ^{••. ,} vt _ ^k t

d . ^{k+ 1 ~ k+ 1 ~ ~ k+ 1} g zie V l '"""" V 2 .._,,,, ••• .._,,,, V N -k- 1 . k+-k+l.

Ostatecznie W = ~ - l.

Można łatwo wykazać (metodą indukcji względem k), ^że algorytm Huffmana prowadzi do konstrukcji drzewa binarnego o minimalnej długości. Na podstawie (2) wystarczy bowiem

wykazać, że sumak wag b. _I (łub innymi słowy, las k-tej sumy częściowej), otrzymanych po k krokach algorytmu Hu ff mana, jest minimalna w porównaniu z ^s~mą k wag bi otrzymanych innym, dowolnym sposobem wykonania k dodawań (patrz [ 11 ], [ f2]).

W ogólnym przypadku, może istnieć wiele drzew binarnych minimalizujących sumę (l ).

Jeśli w ^każdym kroku algorytmu Huffmana elementy następnego ciągu spełniać będą nie- równości v~ ~ ... < vf ^{~ v~ + v~} < vf ⁺ 1, to utworzone drzewo Huffmana minimalizuje jednocześnie m_ax li oraz L N ^{li pośród} wszystkich drzew o minimalnej długości.

l ^~J~N i= I

Przedstawiony ^wyżej problem można zmodyfikować, żądając wyznaczania tzw. minimal- nego drzewa alfabetycznego, tj. takiego które minimalizuje ^sumę (1) przy ustalonej kolejno-

ści (od lewej do prawej) wag wierzchołków końcowych.

Drzewo alfabetyczne o minimalnej ^długości odpowiada sposobowi obliczania sumy W z optymalnie (w sensie oszacowania błędu zaokrągleń) rozstawionymi nawiasami. Algorytm wyznaczania takiego drzewa podany został w pracy [ 11] a jego realizacja w języku ALGOL-60 w pracy [28].

PRZYKŁAD. Na rysunkach 3-5 pokazane są drzewa binarne odpowiadające omawia- · nym wyżej sposobom obliczania sumy W= w ₁ + w ₂ + w ₃ + w 4 = 3 + 7 + 2 + 4. {W wierz-

chołkach pośrednich umieszczone są wartości sum częściowych).

Sposób obliczania W „ od lewej do prawej "

W=lfw ₁ + w ₂ )+w ₃ J+w,_

L ^Wj li= 38·

Rys. 3

Podobnemu problemowi poświęcona jest praca [ 17]. Udowodniono w niej, że dla zmien-

nopozycyjnych liczb binarnych o różnych cechach ^kolejność dodawania ^według algorytmu

Huffmana daje maksymalną dokładność sumy W. (Przez błąd sumy rozumie się obcięcie naj-

mniej znaczących bitów sumy.)

Drzewo o minimalnej ^długo$ci W= {(w ₃ + w ₁ J+w ₄ }+ w ₂

~WJ /J = 30 Rys. 4

Minimalne drzewo alfabetyczne W= (w ₁ + w ₂ J+ (w ₃ + w ₄ J

Ewi ₁₁ = 32 Rys. 5

3. Obliczanie ^wartości wyznacznika i rozwiązywania układu równań liniowych

3.0. D. Konig jest autorem pierwszych prac (patrz [ 13]) poświęconych zastosowaniu gra- fów w algebrze liniowej, a mianowicie do obliczania ^wartości wyznaczników. Następnie wpro- wadzone przez Masona grafy ^sygnałów (signal-flow graph) a przez Coatesa grafy ^przepływu (flow graph) zostały zastosowane przez nich do rozwiązywania równań liniowych (patrz [3], str. 140-320).

W tym rozdziale przedstawimy metodę obliczania wartości wyznacznika macier~y A i

metodę rozwiązywania układu równań liniowych Ax = b

wykorzystujące grafy Coatesa. Metody te ^znajdują zastosowanie w przypadku rozrzedzonych macierzy A lub macierzy A o elementach alfanumerycznych.

I-faktorem grafu D =<X, U>, oznaczmy go przez h, nazywamy regularny graf ^częściowy

stopnia 1 grafu D, tzn. dla każdego wiarzchołka x mamy

l{y: [x, y] € U}I = l{y: [y, x] ^€ U}I = I,

lub bardziej intuicyjnie, I-faktor jest zbiorem rozłącznych konturów grafu D zawierających

wszystkie wierzchołki grafu. 1-faktorowym połączeniem wierzchołków i zj w D, oznaczmy je przez Hii' nazywamy graf częściowy grafu D, który zawiera (i) ^drogę Pij z i do j, (ii) zbiór

rozłącznych konturów zawierających wszystkie wierzchołki grafu z wyjątkiem należących

do Pif

Dla ustalonej macierzy kwadratowej A = (aii) stopnia n grafem przepływu D(A) =

=<X, U; AT) nazywamy obciążony graf zorientowany taki, że IXI =n oraz każdy wierz-

chołekj €X odpowiadaj-temu wierszowi ij-tej kolumnie macierzy A, i ^jeśli aii :#=O, to [i, j] € U oraz ^{wagą łuku} [i, j] jest aii' natomiast jeśli aii = O, to [i, j) 4 ^U.

I jeszcze jedno oznaczenie. Dla podgrafu F grafu D przyjmujemy f(F) = . TJ ^{a;;• lub}

f (F) = O ^jeśli F jest pod grafem pustym. [ 1,1] e F 3.1. Z definicji, wartość wyznacznika macierzy kwadratowej A stopnia n jest równa (1),

gdzie sumowanie rozciągnięte jest na wszystkie permutacje j = (j _{1 ,} j 2 , ••• , j n) ciągu liczb (I, 2, ... , n) a e. _{I 1 , •• „/} . n = ^{1 jeśli j jest} permutacją parzystą i -1 ^jeśli - nieparzystą. Niech ..

a _{1 . / 1} a 2 . • •• a :#= O dla pewnej permutacji j = ^(j 1 , •.• , j n). Wiadomo z teorii grup symetry-

/2 ni

cznych, że każda n permutacja może być jednoznacznie przedstawiona w postaci rodziny cykli

rozłącznych. Tak więc graf częściowy odpowiadający temu składnikowi utworzony z łuków

1/ _{1 ,} ^1), 1/ _{2 ,} ^2], .„, Un' ^{n] jest} sumą rozłącznych konturów grafuD(A) pokrywających wszystkie wierzchołki, a ^więc jest I-faktorem. Podobnie w drugą stronę. Istnieje więc odpo- wiedniość wzajemnie jednoznaczna między różnymi od zera składnikami sumy w (1) i I-fak- torami grafu D(A).

Zauważmy dalej, że znak permutacji

u.„„,jn)=(il,„„in )(k.; ... ,kn ) ... (11' .. „ln )

1 2 • Li

jest równy

(2)

Ostatecznie otrzymujemy

detA = (-l)n ~(-1{hf(h),

h

gdzie h jest I -faktorem graf u D (A), a L h oznacza liczbę konturów w h.

Łatwo można pokazać także, że

(3)

gdzie qh oznacza liczbę parzystej długości konturów w I-faktorze h.

PRZYKŁAD. Niech

(4) [

au A= O

o

Graf ^przepływu dla powyższej macierzy pokazany jest na rysunku 6 a jego I-faktory na ry- sunku 7.

Rys. 6

Rys. 7 Z definicji mamy

natomiast na podstawie (2) otrzymujemy

detA = (-I) ³ [(-1) ³ /(h) + (--1) ² /(g)].

Przejdźmy teraz do wyprowadzenia wzoru na wartość Llii' dopełnienia algebraicznego (i, j)-tego elementu macierzy A. W tym celu zauważmy, że L1ij = A ii, gdzie A ij jest macierzą

kwadratową stopnia n utworzoną z macierzy A przez zastąpienie i-tej kolumny i-tą kolumną

jednostkową e;- Macierzy Aii odpowiada graf D(Aii), który można otrzymać z grafu D(A)

przez usunięcie wszystkich łuków prowadzących z wierzchołkaj (włącznie z pętlą) i doda-

nia łuku U. ^{i] z wagą I. Stosując} do macierzy A ii i grafu D(A ii) wzór (2) otrzymujemy

(5) _{LJij = (-l)n} 2 _{H .. (-I)} ^L _q ^-1 _f(Hii),

l}

gdzie Hij jest I-faktorowym połączeniem z i do jw grafie D(A), a Lq liczbą konturów w Hif Z podobnego rozważania otrzymujemy

(6) _LJii ^/1 _=(-Ir f ^{)' (-1) L h ,-1 f(/1'),}

gdzie h' jest I-faktorem w podgrafie otrzymanym z grafu D(A) przez usunięcie i-tego wierz- chołka wraz z incydentnymi z nim łukami, aLh~ jest liczbą konturów w h ^{1 •}

PRZYKŁAD (c.d). Na rysunku 8 pokazany jest graf D(A 32 ) dla macierzy A określonej przez (4), a na rysunku 9 jedyne jego I-faktorowe połączenie z 3 do 2 w grafie D(A). Stąd natychmiast otrzymujemy

,132 = (-I)3[(-l) ^{1 -} 1/(H32)1 = ^-a11a23·

a,,

Rys. 8

Rys. 9

3.2. Wzory wprowadzone w poprzednim paragrafie zastosujemy teraz do wyznaczenia

rozwiązania układu równań liniowych.

(7)

Wiadomo, ^że dla nieosobliwej macierzy A na podstawie wzorów Cramera mamy n

xk = ^~ bi',:1ik/detA i= 1

(k= 1,2, ... ,n).

Oznaczmy przez Au macierz powstałą z macierzy A przez dodanie (n + I )-szej kolum-.

ny postaci -b i (n + I )-szego wiersza złożonego z n + I zer.

Zauważmy, że

L n ^{b/lik = (-I)n-k} (-l)detA~,

i= 1

gdzie macierz A~ powstaje z macierzy A u prze~ skreślenie ostatniego wiersza i k-tej kolumny.

Ponieważ (-I)n+ ¹ +kdetA~ = L1n+ l ,k(Au), więc na podstawie (5) otrzymujemy

Hn+l,k

gdzie Hn+ l ,k jest 1-faktorowym połączeniem z n + 1 dok w grafie D(Au), a LH jest ^liczbą konturów w Hn+ ₁ ' k' ^.,:.

Ostatecznie, na podstawie (2), (7) i (8) otrzymujemy xk = L (-1{HJ(Hn+l,k){~)-1{hf(h)

H n+ l ,k h .

(9) _(k = ^{1, 2, ...} ,n).

PRZYKŁAD (c.d). Na rysunku 10 pokazany jest graf D(Au) dla macierzy określonej przez ( 4 ), a na rysunku 11 pokazane są jego 1-faktorowe połączenia H _{4 ,2.}

Stąd

On

Rys. 10

X2 = [(-1) ² (-b2)ou033 + (-1 ) ¹ (-b3)0230u )/(o u 023032 - o u 022033) =

= (b3023 - b2033)/(023032 - 022033).

W książce Chena [3), str. 120 znaleźć można wiele innych wzorów na obliczanie warto-

ści det A ,L'.1ii i x k, które różnią się od wyżej przedstawionych tym, że sumowania rozciągnię-

te ^są na inne rodzaje pod graf ów związanych z grafem D (A).

{(;) a,,\._)-"

Rys. 11

3.3. Jak już wyżej wspomnieliśmy, metody obliczania ^wartości wyznacznika i ^rozwiązy­

wania układu równań liniowych oparte ria metodach teorii grafów ^znajdują zastosowanie przy obliczeniach z ^dużymi macierzami rozrzedzonymi lub w przypadku ^obliczeń z macierza- mi o elementach alfanumerycznych. Z takimi macierzami spotykamy ^się np. przy oblicze- niach związanych z obwodami elektrycznymi; elementami macierzy A ^są wtedy funkcje wy- mierne zmiennej zespolonej s i bezpośrednie stosowanie ^którejś z tradycyjnych metod nume- rycznych bardzo komplikuje obliczenia lub czyni je wręcz niemożliwymi. Użyteczność przed- stawionych metod ^zależy od efektywności metod generowania dla danego grafu wszystkich podgrafów, na które rozciągnięte jest sumowanie. Dla wzorów wprowadzonych w tym roz- dziale należy podać metody generowania I-faktorów i I-faktorowych ^połączeń zadanego grafu.

Po pierwsze, moce zbiorów I-faktorów i 1-faktorowych połączeń mogą być wyznaczo- ne przy użyciu (+)-wyznacznika (permanent), który dla macierzy kwadratowej A stopnia n definiujemy w następujący sposób

per A · =La _i 11. . .. _l a . , _n/n

gdzie sumowanie rozciągnięte jest na wszystkie permutacje j = (j 1 , .•• , j n) ^ciągu liczb (I, 2, ... , n).

Niech C(D) = (cii) będzie macierzą kwadratową stopnia n o elementach cii = k, gdzie k jest liczbą łuków z i do jw grafie D. Łatwo wykazać, że liczba I-faktorów grafu D jest równa ^wartości per C(D), a liczba 1-faktorowych ^połączeń z i do jw grafie D jest równa per C(D'ji, gdŻie macierz C(D)ii jest macierzą otrzymaną z macierzy C(D) przez skreślenie j-tego wiersza i i-tej kolumny.

W pracy [ I8] podana jest procedura w ^języku ALGOL-60 obliczania.wartości (+)-wy- znacznika.

· Po drugie, niech Z oznacza zbiór I-faktorów a Zii - zbiór 1-faktorowych ^połączeń z i do jw grafie D oraz Ui i 0 - odpowiednio, zbiory ^łuków grafu D wychodzących z wierz-

chołka j i dochodzących do j.

Wtedy

n n

Z=(P uk)n(P u;),

k=l, /=l

n n Z;i = ( P uk)n ( ^J> uj).

k = ^{l /=I}

k*-i /-:/:.i

gdzie J> oznacza iloczyn kartezjański.

n

Dla dowodu pierwszej ^równości wystarczy zauważyć, że elementami puk ^{są grafy}

• k =I

częściowe takie, ^że z każdego wierzchołka. wychodzi dokład~ie jeden ^łuk, a elementami /> n l1t ^{są grafy częściowe takie, że do} każdego wierzchołka dochodzi dokładnie jeden łuk, /= ^ł

a ^więc elementami zbioru Z ^są 1-faktory. Podobnie można wykazać prawdziwość drugiej

równości.

4. Redukcja macierzy rozrzedzonych

4.0. W tym rozdziale zajmować się będziemy metodami sprowadzania zadanej macie- rzy rozrzedzonej A do postaci, dla której ^możliwa jest redukcja kosztu (czasu i/lub ^pamięci)

ob.liczeń koniecznych dla otrzymania rozwiązania postawionego problemu.

Na wstępie podamy pewną liczbę nowych ^pojęć i faktów z teorii grafów. Graf zoriento- wany D = ^(X, U> nazywamy bez konturowym, ^jeśli nie zawiera konturu. Łatwo można wy-

kazać, że graf jest bezkonturowy wtedy i tylko wtedy, gdy istnieje tzw. uporządkowana nu- meracja jego wierzchołków, tj. jeśli [i, jJ ^{EU, to i <j.} Oznaczmy przez D* = ( Y, V> graf zorientowany będący kondensacją grafu D, tj. zbiór wierzchołków Y grafu D* odpowiada składowym silnej spójności grafu Di dla y ₁ :Yi ^{E Y zachodzi}

[y,,y2)E v=3x1 EY1 3x2 EY2 [X1,X2JEU.

Można wykazać także, że graf jest bezkonturowy wtedy i tylko wtedy, gdy D* = D.

Określmy dla danego graf u D jeszcze jeden graf pochodny. Tranzytywnym domknięciem

grafu D nazywamy graf D ^{= <X,} U>,<Ua którego zbiór łuków zdefiniowany jest w następują­

cy sposób: dla x 1 , x 2 EX mamy [x 1 , x 2 ] E U wtedy i tylko wtedy, gdy w grafie D istnieje droga z x ₁ do x _{2 .} Graf D =<X, U> nazywamy grafem dwudzielnym lub grafem Koniga, ^jeśli istnieje podział zbioru wierzchołków na dwa niepuste i rozłączne podzbiory X 1 , X 2 takie,

że jeśli [x _{1 ,} x 2] EU, tox ¹ E X ₁ ix ₂ e X ₂ lubx 1 e X ₂ ix ₂ e X _{1 •}

Niech P będzie macierzą permutacyjną stopnia n. Wiadomo, że p-i =PT, oraz macierz PAPT powstaje z macierzy A przez dokonanie tej samej permutacji na jej wierszach i kolum- nach. Macierz A= PAPT nazywać będziemy permutacją macierzy A.

~acierz A można rozłożyć, jeśli istnieje jej permutacja postaci

gdzie A ii (i = I, 2, ... , k) są macierzami kwadratowymi. Natomiast macierz A jest redukowal- na, jeśli istnieje jej permutacja postaci

(2) _A - = lA11 ^A12 ~lkl

A22 . .

O Akk

Rozkład lub redukcja macierzy A jest bardzo przydatna w wielu zagadnieniach metod numerycznych. Pozwala bowiem na ^redukcję rozmiaru zagadnienia.

I tak, dla macierzy A dającej się rozłożyć lub zredukować mamy:

k • (i) detA = JI. ^detA;;,

i= I k

(ii) det (A - Al) = Jl ^{det (A;; - AJ), a więc zbiór} wartości własnych macierzy A jest i= I

sumą zbiorów wartości własnych pod macierzy A 11 , ••• , A kk .

(iii) A- 1 ₌ pT(Ar1 P. Na przykład, w przypadku macierzy redukowalnej dla k = ²

otrzymujemy

A -I A-1]

- 11 ^{·A 12 • 22 •} A-1 22

(iv) Dla macierzy A dającej się rozłożyć układ równań liniowych Ax = b rozpada się

na niezależne podukłady. Wystarczy bowiem rozwiązać układ Ay =Pb, a stąd x = pTy.

Dla zagadnień (i), (iii) oraz (iv) możliwość rozkładu na zagadnienia o mniejszych roz- miarach ma miejsce także w przypadku doprowadzenia macierzy A do uogólnionej postaci rozłożonej lub uogólnionej postaci zredukowanej A = ^PAQ. ·

Dla przykładu, dla macierzy

A = [~ ~ ~ ~i _{o· o} 3 2

I O 2 O

nie istnieje postać zredukowana A = PAPT, gdyż graf D(A) jest silnie spójny, zawiera bo- wiem kontur dłogości n: (I, 2, 3, 4, 1). Istnieje natomiast dla macierzy A uogólniona postać

zredukowana A= PAQ, gdzie

[ I O P= O O O I

o o O ¹ Ol ^o

o o ' o ¹

1 o o o

[

O O O Il

Q = .o o ¹ o .'

O·l O O

skąd A =r~1+ _L~J o o 1 5 .

O O 2 I

Metoda otrzymywania macierzy permutacyjnych Pi Q dla ^rozkładu macierzy A podana

została w pracy [22], a dla redukcji - w pracach [5], [6] i [19]. Wszystkie te metody wyzna-

czają macierze P i Q korzystając z informacji o strukturze macierzy A uzyskanych z różnego

typu graf ów przyporządkowanych tej macierzy.

Procedury w języku ALGOL-60 realizujące metody wyznaczania macierzy P lub Pi Q

zostaną opublikowane w czasopiśmie Zastosowania Matematyki.

4.1. W tym paragrafie przedstawimy ^metodę wyznaczania macierzy permutacyjnej P, takiej że macierz A jest postaci (2) (patrz [7], [8]).

Przyporządkujemy macierzy A graf zorientowany D(A) =<X, U> taki, ^że IXI =n i wierzchołek i-ty tego grafu odpowiada i-temu wierszowi i i-tej kolumnie macierzy A oraz [i, i] e ^U~ aii =F O ^A i =F j. Oznaczmy przez M(A) lub krócej przez M macierz incydencji grafu D(A), tj. macierz o elementach mii = 1 jeśli [i, j] e U i mii = ^{O jeśli [i, j] łf U.}

Algorytm wyznaczania postaci zredukowanej macierzy opiera ^się na fakcie, ^że maksy- malne nieredukowalne podmacierze macierzy A odpowiadają wzajemnie jednoznacznie

składowym silnej spójności grafu D(A) (patrz [7], [8]).

Algorytm

K rok 1. Dla danej macierzy A ^utworzyć macierz incydencji M grafu D(A).

Kr ok 2. Utworzyć macierz incydencji M tranzytywnego domknięciaD(A) grafu D(A).

Bardzo ^efektywną ze względu na czas ^obliczeń i obszar zajmowanej ^pamięci maszyny cyfrowej metodą wyznaczania macierzy Mjest metoda Warshalla (patrz [14], procedura paths). Niech 1jM = (r _1k), gdzie _r1k ~ _m1k U mu n ^{mik (/, k =} 1, 2, ... , n), wtedy

M = TnTn-l ... T V ₁ M.

Łatwo można sprawdzić, że liczba działań wynosi co najwyżej 2n ³ oraz, że macierzM

może być utworzona na miejscu macierzy M (patrz [14]).

Kr ok 3. Wyznaczyć składowe silnej ^spójności grafu D(A). Oznaczmy je p·rzez Dl, ... ,Dk. Niech s =MX _MT= (s;j), gdzie sij = m;;. mji (i, i= 1, 2, ... ,n). Widzimy, że

wierzchołki i, i ^należą do tej samej ^składowej silnej spójności wtedy i tylko wtedy gdy sii = 1.

Tak więc z wierzchołkiemj do tej samej składowej silńej spójności należą wierzchołki

{/: sil= 1}.

Jeśli k = 1, to macierz A nie jest redukowalna (koniec algorytmu).

Kr ok 4. Utworzyć kondensację D* = ( ^{Y, V>} grafuD(A), tj. Y = ^{D 1 ,D _{2 , •••} ,Dn}

oraz [Di' Di] e V wtedy i tylko wtedy, gdy ^istnieją g e Di oraz h e Di takie, że [g, h] e U.

Jeśli V= </J, to dla macierzy A istnieje jej permutacja postaci (1).

Kr ok 5. Wyznaczyć macierz permutacyjną Q stopnia k taką, że macierz QM(D*)QT jest macierzą górną trójką tnę. Istnienie m$ierzy Q wynika z bezkonturowości grafu D*, a wyznaczenie jej jest równoznaczne z W}znacz;niem uporządkowanej numeracji składo­

wych silnej spójności grafu D(A) (patrz [14), funkcja niestandardowa circuit).

Niech D. , D. , ... , D. będzie tym uporządkowaniem.

'1 '2 l k

Kr ok 6. ^Wyznaczyć macierz permutacyjną P. W tym celu należy ponumerować kolej- nymi liczbami 1, 2, ... , n wszystkie wierzchołki kolejnych składowych spójności D; , D. ₁₂ , .„,D . . _{1 k} Niech (j _{1 , {,} i _{2 , •••} ,jn) będzie nową numeracją wierzchołków (1, 2, „., n). 1

Wtedy macierz permutacyjna P jest postaci

P=[e. ,e. , ... ,e. ],

11 Iz 'n gdzie ei jest i-tym wektorem jednostkowym.

PRZYKŁAD. Niech macierz A będzie postaci 1 o ^{1 1} o o ^{l 1} o o o

A= o o 1 o o

1 1 1 1 1 1 o ^{1 1} o o o o

Odpowiadający jej graf D(A) pokazany jest na rys. 12, a macierz incydencji M grafu D(A)

różni się od macierzy A tylko elementami głównej przekątnej (m;; = O dla i = 1, 2, ... , n).

Rys. 12

Graf D(A) składa się z trzech składowych silnej spójności D ₁ = ^{2ł,D 2 ⁼ {1, 4, Sfi D ₃ = {3, 6}, dla których graf kondensacji D* pokazany jest na rysunku 13.

Rys. 13

Uporządkowanie składowych silnej ^spójności jest następujące D _{2 ,} D _{1 ,} D _{3 ,} skąd otrzymuje- my nową numerację wierzchołków

stara numeracja 1 2 3 4 5 6 nowa numeracja 1 4 5 2 3 6, i macierz permutacyjną P:

1 o o o o o o o o ¹ o o

P= o o o o ^I o o ^I o o o o o o ^I o o o o o o o o ^I

Ostatecznie

To, czy dana macierz jest redukowalna ^zależy od liczby elementów zerowych i od ich rozmieszczenia. Dla· przykładu, może być nieredukowalna macierz o n ² -n zerach, której od- powiada graf ^będący konturem o długości n, ^{a może być} redukowalna macierz o n - I ze- rach (które w macierzy A ^są n - I pierwszymi elementami ostatniego wiersza lub ostatniej kolumny). W praktyce ^są nawet ^kłopoty z generowaniem losowych macierzy, którym odpo-

wiadają grafy nie silnie spójne. Prawdziwy jest bowiem następujący fakt. Niech dla dowolne- go p e (O, I] element macierzy kwadratowej A stopnia n będzie różny od zera z prawdopo- dobieństwem p. Wtedy prawdopodobieństwo tego, że graf D(A) macierzy A jest silnie spój·

ny ^dąży do 1 przy n dążącym do ^{00 •}

4.2. W poprzednim paragrafie omówione ^zostały metody pozwalające w przypadkach zagadnień (i)-(iv) na ich ^rozkład na niezależne zagadnienia tego samego typu o mniejszych . wymiarach. Jeśli rozkład zagadnienia nie jest ^możliwy, to należy zbadać, czy nie istnieje uogólniona permutacja A = ^PAQ macierzy A, dla której obliczenia ulegają uproszczeniu.

Załóżmy, że pamiętanie macierzy rozrzedzonej A (lub jej uogólnionej permutacji), dla której wykonany zostanie odpowiedni algorytm, polega na pamiętaniu jej elementów ^róż­

nych od zera (packed form). Taki sposób pamiętania zapewniam. in. (a) pamiętanie całej macierzy A w pamięci operacyjnej lub jeśli nie jest to możliwe, to (b) szybką transmisję między pamięcią operacyjną a zewnętrznymi, (c) zmniejszenie czasu obliczeń dzięki wyko·

nywaniu tylko koniecznych operacji arytmetycznych, tj. takich których argumentami nie są zera. Można więc zażądąć takiej uogólnionej permutacji A macierzy A, dla której zasto- sowanie odpowiedniego algorytmu wprowadza minimalną liczbę elementów ^różnych od zera na miejsce zer. Dla przykładu rozpatrzmy pierwszy etap eliminacji Gaussa (GE).

Wiadomo na przykład, że jeśli macierz A jest postaci

A11 A12

A22

A= o

Akl Ak2

gdzie podmacierze Au (i = 1, 2, ... , k) są nieosobliwe, to GE nie wprowadza nowych blo- ków niezerowych.

Wiadomo także, że jeśli macierz A jest macierzą wstęgową (w skrócie BF), tzn. jeśli aii =I= O, to (i, j) EI= l ^{(i, j):} I i-jl ~A.} a dla (i, j) t ^{I mamy aii} = O, oraz jako elementy głów­

ne wybierane są elementy diagonalne, to GE może utworzyć nowe elementy różne od zera tylko dla (i, j) ^{E /.}

O innych postaciach macierzy A upraszczających obliczenia (nie tylko GE) traktują pra- ce [24] i [25).

Dalszą część tego paragrafu poświęcimy na omówienie sprowadzania macierzy A do ma- cierzy A postaci BF.

Jeśli macierz A jest symetryczna i dodatnio określona, to celowym jest zachowanie syme- trii w macierzy A, ^{a więc} w tym przypadku należy założyć, że Q =PT. W pracach [1] i [4]

podane zostały efektywne metody przybliżone wyznaczania macierzy permutacji P. Wyko-

rzystują one informacje o strukturze macierzy A dające się otrzymać z grafu postaci <X, U>

przyporządkowanego tej macierzy w następujący _{. .} sposób _. [i,. i] ^E U= a .. _ij =I= O.

Autorowi nie jest znana żadna praca poświęcona sprowadzaniu dowolnej macierzy A (nie koniecznie symetrycznej) do postaci wstęgowej o minimalnej szerokości A.

Przedstawimy teraz pewien problem teorii grafów (patrz [27]), którego rozwiązanie (je- szcze nie znane!) może być pomocne w tym przypadku.

Niech A będzie macierzą prostokątną m x n, a K(A) grafem dwudzielnym przyporząd­

kow~nym tej macierzy, tj. K(A) = (XU Y, U), gdzie IXI = m, I YI =n oraz dla i EX, j E Y mamy [i, j) E U= aii =I= O. Załóżmy, że wierzchołki ze zbiorów X i Y ponumerowane są od- powiednio kolejnymi liczbami 1, 2, ... , m i 1, 2, ... , n. ^{Widać, że} zmiana kolejności wierszy (kolumn) macierzy A odpowiada zmianie numeracji wierzchołków X (odpowiednio, Y).

Niech k = ^(k 1 , k 2 , .•• ,km) i I= (/ 1 , / 2 , ••• ,In) ^będą odpowiednio numeracjami wierzchoł- ków X i Y. Określmy szerokość wstęgi macierzy A przez

v =max lk; - Iii'

gdzie max jest brane po wszystkich numerach k;, li wierzchołków połączonych. Przyjmijmy,

że graf K (A) narysowany jest na płaszczyźnie w następujący sposób: Elementy zbiorów

X i Y ^leżą na dwóch prostych równoległych w kolejności zgodnej z numeracją, natomiast

elementy zbioru U są odcinkami prostej przechodzącej przez odpowiednie punkty (wierz-

chołki) na prostych równoległych. Dla grafu K(A) z numeracją wierzchołków k, I oznacz-

my przez c(k, I) liczbę par łuków przecinających się nie w wierzchołkach. Interesować

nas będzie wyznaczenie numeracji wierzchołków k* i I* takich, że

c(k*, /*) = minc(k, [).

k,l

Zauważmy, że każdej parze łuków grafu K(A) przecinających się nie w wierzchołkach

odpowiada w macierzy A podmacierz postaci ..

(3)

gdzie a;·* O i agh *O, a gwiazdką oznaczono elementy o nieistotnej ^wartości. Tak więc mi-

/

nimalizacja c(k, [)może być zastąpiona przez minimalizację liczby podmacierzy (3) w ma- cierzy po wszystkich macierzach otrzymanych z A przeż permutację wierszy i permutację

kolumn.

Nie jest znany algorytm wyznaczania numeracji k* i I*. Intuicyjnie wydaje ^{się że} przy szukaniu numeracji k* i /* należy dążyć do przemieszczania niezerowych elementów ma- cierzy A jak najbliżej przekątnej.

Otwartym problemem jest podanie ewentualnego związku między numeracjami k* i I*

orai numeracjami k' i / ¹ minimalizującymi szerokość wstęgi. Można podać przykład, że mi- nimalna szerokość wstęgi macierzy nie zapewnia minimalnej liczby podmacierzy typu (3) w tej macierzy. Autor przypuszcza, ^że prawdziwa jest ^zależność odwrotna.

Bibliografia

[ l J I. Ar a ny, W. F. Sm y t h and L. S z od a, An improved method for reducing the bandwidth of sparse symmetric matrices, lnformation Processing 71, Proceedings I FIP Conf erence 1971, vol. 2, str. 1246-1250.

[ 2] W: K. Che n, On directed graph solution of linear a/gebraic equations, SIAM Rcview 9(1967), str. 692-707.

[ 3] - Applied Graph Theory, Amsterdam 1971.

(4 l E. Cu t hi l l, Severa/ strategies forreducing the bandwidth of matrices, w [ 16 ), str. 157-166.

[ 5] A. L. D u l m age and N. S. M e n de l s o h n, On the inversion of sparse matrices, Math.

Comp. 16(1962), str. 494-496.

[ 6) - - Graphs and matrices, w Graph Theory and Theoretica/ Physics (red. Harary), 1967, str.•

167-227.

(7) F. Ha r ary, A graph theoretic method for the complete reduction of a matrix with a view to- ward finding its eigenvalues, J. Math. Physics 39(1959), str. 104-111.

[ 8] - A graph theoretic approach to matrix inversion by partitioning, Num . Math. 4 (1962), str.

128-135.

[ 9] - R. N o r m a n and D. C a r t w r i g h t, Structural Mode/s: An Introduction to the Theory of Directed Graph, New York 1965.

[ 10] - Sparse Matrices and graph theory, w [ 15), str. 139-150.

[ 11] T. C. Hu and A. C. Tu ck er, Optima/ computer search trees and variable-length alphabeti- cal codes, SIAM J. Appł~ Math. 21 (1971), str. 514-532.

[ 12) D. E. Knut h, The Art ofComputer Programming, vol. l, New York 1968.

[ 13 I ^D. Ko ^{n i g,} Theorie der endlichen und unendlichen Graphen. New York 1950.

[ 14 J J. Kuch arc z y k i M. M. S y ^sł o, Realizacja algorytmów optymalizacji w ^języku ALGOL-60, Warszawa 1974.

[ 15 I J. K. R e i d (red.), Large Sparse Sets of Linear Equations. London 1970.

[ 16] D. J. Rose and R. A. W il I o u g h by (red.), Sparse Matrices and Their Applications, New

York 1972.

[ 17 J A. S c h u r m a n n, On the minimum error in addition processes of positive jloating-poit num- bers, Zastosow. Matem. 13(1973), str. 351-364.

[18) B. S hi r i v er, P. J. E ber Ie i n and R. D. Dix o n, Algorithm 361. Permanent Junction of a square matrix I and 1/, Comm. ACM 12(1969), str. 634.

[19] D. V.Stew ar d,Partitioningand tPnringsystemsofequationa, SIAM J. Numer. Anal. 2(1965), str. 345..,..365.

[20) M. M. S y sł o, Remarks on addition processes of positive floating-point numbers, Zastosow.

Matem. 14(1974).

[21] L. Szam ^koło w i cz, Teoria grafów skończonych, Wrocław 1971.

[22] R. P. Te w ars o n, Row-column permutation ofsparse matrices, Computer J. 10(1967), str.

300-305.

[23) - The product form of inverses ofsparse matrices and graph theory, SIAM Review 9(1967), str. 91-99.

L24] - Computations with sparse matrices, SIAM Review 12(1970), str. 527-543.

( 25] - Sorting and ordering Sparse linear systems, w [ 15], str. 151-16 7.

[26] - SparseMatrices, New York 1973.

[ 27) M. E. W a t ki n s, A sper:ial crossing number for bipartite graphs: a research problem. Ann. N.

York Acad. Sc. 175 (1970), str. 405-410.

[28] M. Yo he, Algorithm 428. An Algol procedure for the Hu-Tucker minimum redundancy alpha- betic coding met!iod, Comm. ACM 15 (1972), str. 360.

[29) A. A. Z y ko w, Teoria grafów skończonych (ros.), Nowosybirsk 1969.