Dunnbesetzte Matrizen: Algorithmen und¨ Datenstrukturen · 1 Motivation, Beispiele und...

Dunnbesetzte Matrizen: Algorithmen undDatenstrukturen

1

Inhaltsverzeichnis

1 Motivation, Beispiele und Definitionen 41.1 Lagrange-Polynome . . . . . . . . . . . . . . . . . . . . . . . . 41.2 Was ist eine dunnbesetzte Matrix? . . . . . . . . . . . . . . . 6

1.2.1 Wie kommt man auf ein dunnbesetztes Problem? . . . 81.3 Wiederholung . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

1.3.1 Vektoren . . . . . . . . . . . . . . . . . . . . . . . . . . 81.3.2 Matrizen . . . . . . . . . . . . . . . . . . . . . . . . . . 91.3.3 Lineare Gleichungssysteme . . . . . . . . . . . . . . . . 13

1.4 Graphen und dunnbesetzte Matrizen . . . . . . . . . . . . . . 16

2 Dunnbesetzte Vektoren und Matrizen 242.1 Dunnbesetzte Vektoren . . . . . . . . . . . . . . . . . . . . . . 24

2.1.1 Das SDOT-Problem . . . . . . . . . . . . . . . . . . . 252.1.2 Das SAXPY-Problem . . . . . . . . . . . . . . . . . . . 26

2.2 Speicherung von dunnbesetzten Matrizen . . . . . . . . . . . . 272.2.1 Koordinatenform . . . . . . . . . . . . . . . . . . . . . 282.2.2 Compressed Sparse Row Format (CSR) . . . . . . . . . 282.2.3 CSR mit extrahierter Diagonale . . . . . . . . . . . . . 292.2.4 Diagonalweise Speicherung . . . . . . . . . . . . . . . . 302.2.5 Jagged Diagonalformat . . . . . . . . . . . . . . . . . . 312.2.6 Zusammenfassung . . . . . . . . . . . . . . . . . . . . . 32

2.3 Elementare Matrixoperationen . . . . . . . . . . . . . . . . . . 332.3.1 Das Matrix/Vektor-Produkt . . . . . . . . . . . . . . . 332.3.2 Das Matrix/Matrix-Produkt . . . . . . . . . . . . . . . 35

2.4 Algorithmen . . . . . . . . . . . . . . . . . . . . . . . . . . . . 362.4.1 Matrix/Vektor-Produkt . . . . . . . . . . . . . . . . . 362.4.2 Matrix/Vektor-Produkt mit vertauschtem i und j . . . 372.4.3 Matrix/Matrix-Produkt . . . . . . . . . . . . . . . . . 372.4.4 Matrix/Matrix-Produkt mit vertauschtem k und j . . . 37

3 Gauss-Elimination 383.1 kij-Standardform . . . . . . . . . . . . . . . . . . . . . . . . . 383.2 ikj-Form . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 393.3 jki-Form . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 403.4 Crout-Form . . . . . . . . . . . . . . . . . . . . . . . . . . . . 403.5 Verfahren fur symmetrisch positiv definite Matrizen . . . . . . 41

3.5.1 Borderingform . . . . . . . . . . . . . . . . . . . . . . . 423.5.2 Innere Produktform . . . . . . . . . . . . . . . . . . . . 42

2

4 Band- und Envelope-Methoden fur symmetrisch positiv de-finite Matrizen 444.1 Band . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 444.2 Envelope . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 454.3 Der Cuthill-McKee-Algorithmus . . . . . . . . . . . . . . . . . 484.4 Der Reverse Cuthill-McKee-Algorithmus . . . . . . . . . . . . 50

5 Gauss-Elimination und der Eliminationsgraph 525.1 Quotient-Tree-Methoden . . . . . . . . . . . . . . . . . . . . . 565.2 Dissection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 585.3 Lokale Pivotstrategien . . . . . . . . . . . . . . . . . . . . . . 60

5.3.1 Hauptkriterien fur lokale Pivot-Strategien . . . . . . . 615.3.2 Die algebraische Struktur der Inversen . . . . . . . . . 62

6 Iterative Verfahren zur Losung linearer Gleichungssysteme 656.1 Das GMRES-Verfahren fur allgemeine quadratische Matrizen . 65

6.1.1 GMRES Restarted . . . . . . . . . . . . . . . . . . . . 676.1.2 Fehlerabschatzung . . . . . . . . . . . . . . . . . . . . 68

6.2 Prakonditionierung . . . . . . . . . . . . . . . . . . . . . . . . 696.2.1 Wesentliche Prakonditionierer . . . . . . . . . . . . . . 69

6.3 Zusammenfassung . . . . . . . . . . . . . . . . . . . . . . . . . 71

A Literatur 72

3

1 Motivation, Beispiele und Definitionen

Beispiel 1.1

Gegeben sind n + 1 Stutzstellen (xi, yi), 0 ≤ i ≤ n. Gesucht ist ein inter-polierendes Polynom p(x) vom Grad n mit p(xi) = yi, 0 ≤ i ≤ n, alsop(x) = a0 + a1x + . . . + anx

n, das die folgenden Bedingungen erfullt:

a0 + a1x0 + . . . + anxn0 = y0 = p(x0)

......

...a0 + a1xn + . . . + anx

nn = yn = p(xn)

Schreibweise: 1 x0 . . . xn

0

1 x1 . . . xn1

......

1 xn . . . xnn

a0

...an

=

y0...

yn

Das lineare Gleichungssystem hangt von der gewahlten Basis fur Polynomeab, hier: 1, x, x2, . . . xn

Es kann eine beliebige Basis q0(x), . . . , qn(x) fur die Polynome gewahlt wer-den, z.B. 1, 1 + x, 1 + x + x2.

Damit ergibt sich:

p(x) =∑n

i=0 aiqi(x) und p(xk) =∑n

i=0 aiqi(xk), 0 ≤ k ≤ n

In Matrix-Schreibweise: q0(x0) q1(x0) . . . qn(x0)...

...q0(xn) q1(xn) . . . qn(xn)

a0

...an

=

y0...

yn

2

1.1 Lagrange-Polynome

Lagrange-Polynome sind wie folgt definiert:

Lj(x) =(x−x0)...(x−xj−1)(x−xj+1)...(x−xn)

(xj−x0)...(xj−xj−1)(xj−xj+1)...(xj−xn), 0 ≤ j ≤ n, Grad ≤ n

Fur alle Lagrange-Polynome gilt an der Stelle xk:

Lj(xk) =

{0 j 6= k1 j = k

4

Damit ergibt sich fur das Interpolationsproblem: L0(x0) L1(x0) . . . Ln(x0)...

...L0(xn) L1(xn) . . . Ln(xn)

a0

...an

=

y0...

yn

Durch die oben beschriebenen Eigenschaften der Lagrange-Polynome ergibtsich ai = yi als Losung und p(x) =

∑ni=0 yiLi(x) als Polynom.

Im allgemeinen sind derartige Interpolationsprobleme jedoch ungleich kom-plexer. So enthalt ein kompliziertes Problem viele Punkte und fuhrt damitzu einem großen linearen Gleichungssystem. Wird zusatzlich eine hohe Ge-nauigkeit verlangt, werden ebenfalls mehr Stutzstellen benotigt.

Eine n×n-Matrix benotigt O(n2) zur Speicherung und O(n3) fur die Losungdes Gleichungssystems. Ein solches Gleichungssystem ist daher i.a. nicht mitvertretbarem Aufwand losbar.

Auswege:

1. Wahl der Basisfunktion: sin, cos, exp oder geeignet strukturierte Matri-zen, z.B.:

(a) Toeplitz-Matrix t0 t1 tn

. . . . . .

t−1 t1. . . . . .

t−n t−1 t0

Mit Hilfe der schnellen Fourier-Transformation (FFT) lasst sichein solches Gleichungssystem mit Kosten O(n · log n) losen.

(b) Vandermonde-Matrix 1 x0 . . . xn0

......

1 xn . . . xnn

Dieses Gleichungssystem kann mit Kosten O(n) fur die Speiche-rung und O(n2) oder schneller fur die Losung berechnet werden.

5

2. Parallelrechner

3. Reduktion von 2D-Problemen auf 1D-Randprobleme

4. Geeignete Wahl der Basis, so daß die Matrix dunnbesetzt ist. Einesolche Matrix lasst sich mit Kosten O(n) speichern. Der Idealfall isteine Diagonalmatrix, die sich in O(n) losen lasst:

d1 0d2

. . .

0 dn

a1

...an

=

b1...bn

Als Losung ergibt sich trivialerweise ai = bi

di

Die Moglichkeiten 2 - 4 sind sich jedoch in der Praxis schwierig und daherselten.

1.2 Was ist eine dunnbesetzte Matrix?

Der Begriff dunnbesetzte Matrix ist nicht eindeutig definiert. Es gibt diefolgenden teilweise informellen Definitionen:

1. Eine Matrix heißt dunnbesetzt, wenn sie soviele Nullen enthalt, daßman einen wesentlichen Vorteil davon hat, diese Nulleintrage zu beruck-sichtigen

2. Eine Matrix A heißt dunnbesetzt, wenn nnz(A) = O(n) gilt, wobeinnz(A) die Zahl der Eintrage 6= 0 (Number Of Nonzeros) ist

3. Eine Matrix A heißt dunnbesetzt, wenn eine Zahl p � n exitstiert,so daß jede Zeile/Spalte von A hochstens p Elemente 6= 0 (Nonzeros)enthalt

6

Beispiel 1.2

1.

* * . . . ** * 0 . . . 0... 0 *

. . ....

.... . . . . . 0

* 0 . . . 0 *

Nicht dunnbesetzt nach Definition 3, da die Matrix nach dem erstenGauss-Schritt fast vollbesetzt ist.Die Vertauschung der 1. Zeile/Spalte mit der n. Zeile/Spalte ermoglichteine Losung mit Gausselimination in O(n)

2. A =

a1 a2

a3 * * . . . ** * 0 . . . 0 *... 0

. . . . . ....

......

. . . . . . 00 . . . 0 * *

a4 * * . . . * *

= D + a1e

T1 + a2e

Tn + e1a

T3 + ena

T4 , ei = (0 . . . 0, 1, 0 . . . 0)T

= D + U4VT4︸︷︷︸

Rang 4

Sherman-Morrison-Woodbury:

A−1 = D−1 −D−1U4 (I + V T4 D−1U4)

−1︸︷︷︸4 × 4-Matrix

V T4 D−1

3.

* * 0

*. . . . . .. . . . . . *

0 * *

⇒

* * 0

*. . . . . .. . . . . . *

0 * *

* * 0

*. . . . . .. . . . . . *

0 * ** * 0

*. . . . . .. . . . . . *

0 * *

. . .

triagonal blocktriagonal

2

7

1.2.1 Wie kommt man auf ein dunnbesetztes Problem?

Beispiel 1.3

u′′(x) = f(x), u(0) = u(1) = 0

Losung mit Ansatzfunktionen qj(x), 1 ≤ j ≤ n, qj(0) = qj(1) = 0:u′′(x)− f(x) = 0 ⇔ (u′′(x)− f(x))qj(x) = 0, 1 ≤ j ≤ n ⇔∫ 1

0(u′′(x)− f(x))qj(x) dx = 0 ⇔∫ 1

0u′′(x)︸︷︷︸

=Pn

j=1 ajq′′j (x)

qk(x) dx =∫ 1

0f(x)qk(x) dx ⇔

∫ 1

0

∑nj=1 ajq

′′j (x)qk(x) dx =

∫ 1

0f(x)qk(x) dx, 1 ≤ k ≤ n ⇔∑n

j=1 aj[∫ 1

0q′′j (x)qk(x) dx] =

∫ 1

0f(x)qk(x) dx, 1 ≤ k ≤ n

Als Matrix: ( ∫ 1

0q′′1q1dx

∫ 1

0q′′1q2dx . . .

...

) a1...

an

=

b1...bn

Durch eine geeignete Wahl der qj mit kompaktem Trager lasst sich eine dunn-besetzte Matrix erreichen. Werden die qj so gewahlt, daß fur fast alle j, k∫ 1

0q′′j qk = 0 gilt, so ist A dunnbesetzt.

2

1.3 Wiederholung

1.3.1 Vektoren

x ∈ Rn, x =

x1...

xn

xT = (x1, . . . xn)

x + y =

x1 + y1...

xn + yn

8

Skalarprodukt:

(x, y) = x · y = xT y = yT x = (x1, . . . xn)

y1...

yn

=∑n

i=1 xiyi

Vektornormen:‖ x ‖2 =

√(x, x) =

√∑ni=1 x2

i (Euklidsche Norm)‖ x ‖∞ = maxn

i=1|xi|‖ x ‖1 =

∑ni=1 |xi|

Winkel zwischen zwei Vektoren:cos(∠ x, y) = xT y

‖x‖2·‖y‖2≤ 1

x ‖ y: ϕ = 0, cos ϕ = 1 ⇔ xT y = ‖ x ‖2 · ‖ y ‖2

x, y orthgonal: ϕ = π2, cos ϕ = 0 ⇔ xT y = 0

Satz 1.1 (Cauchy-Schwarz)

|xT y| ≤ ‖ x ‖2 · ‖ y ‖2 (ohne Beweis)

1.3.2 Matrizen

A ∈ Rn,m, n×m-Matrix

Produkt AB = C:

1 . . . m

a11 . . . a1m...

...ar1 . . . arm...

...an1 . . . anm

·

1 b11 . . . b1s . . . b1k

......

......

m bm1 . . . bms . . . bmk

=

c11 . . . c1k...

...cn1 . . . cnk

(n×m) · (m× k) = n× k

crs =∑m

j=1 arjbjs

C =

∑m

j=1 a1jbj1 . . .∑m

j=1 a1jbjk

......∑m

j=1 anjbj1 . . .∑m

j=1 anjbjk

9

Produkt zwischen Vektoren:

(a1, . . . am)

b1...

bm

=∑m

j=1 ajbj

x1...

xn

︸︷︷︸

r

(y1, . . . yk)︸︷︷︸s

=

x1y1 . . . x1yk...

...xny1 . . . xnyk

Satz 1.2 (Sherman-Morrison-Woodbury)

(A + uvT︸︷︷︸Rang 1

)−1

= A−1 − 11+vT A−1u

A−1uvT A−1, A ∈ Rn×n

Beweis:A−1A = I(A−1 − 1

1+vT A−1uA−1uvT A−1) · (A + uvT ) =

= I + A−1uvT − 11+vT A−1u

A−1uvT − 11+vT A−1u

A−1u (vT A−1u)︸︷︷︸∈R

vT

= I + A−1uvT − 11+vT A−1u

A−1uvT − vT A−1u1+vT A−1u

A−1uvT

= I + 11+vT A−1u

(1 + vT A−1u− 1− vT A−1u)︸︷︷︸=0

A−1uvT

= I

2

Eine Matrix kann als lineare Abbildung aufgefasst werden:

y = Ax: x ∈ Rn → y = Ax ∈ Rm

Matrixnormen:‖ A ‖2 = supx6=0

‖Ax‖2‖x‖2

‖ A ‖22 = sup

‖Ax‖22‖x‖22

= supx 6=0xT AT Ax

xT x

Frobenius-Norm: ‖ A ‖F =√∑

i,j a2ij

10

Zu einer solchen Abbildung lasst sich ein Bildraum und dessen Dimensiondefinieren:

Bild(A) := {y = Ax|x ∈ Rn} ⊆ Rm

dim(Bild(A)) := rang(A) ≤ n, m

Außerdem gilt in einer n× n-Matrix:det(A) = 0 ⇔ rang(A) < n

Beispiel 1.4

x =

(x1

x2

), A = (1 0)

y = Ax = (1 0)

(x1

x2

)= x1 ∈ R1

In diesem Beispiel wird also der R2 auf den R1 abgebbildet. Außerdem isttrivialerweise rang(A) = 1.

2

Eigenwert:Ausgehend von der Gleichung Ax = λx, heißen die λ, die die Gleichungerfullen Eigenwerte, die dazugehorigen Vektoren x 6= 0 Eigenvektoren. EinEigenvektor ist eine Fixgerade in der durch A beschriebenen Abbildung.

y = Ax = λx ⇔ (A− λI)x = 0 ⇔ det(A− λI) = 0

Aus der Berechnung dieser Gleichung ergeben sich die Eigenwerte als Null-stellen des charakteristischen Polynoms.

Basis:Die Vektoren x1, . . . xn sind eine Basis des Rn, wenn sie folgende Bedingungenerfullen

• lineare Unabhangigkeit:∑n

j=1 αjxj = 0 ⇒ ∀j: αj = 0

• jeder Vektor x ∈ Rn kann durch eine Linearkombination des Basisvek-toren dargestellt werden: x =

∑nj=1 αjxj mit passenden αj

11

Orthonormalbasis:Die Vektoren x1, . . . xn bilden eine Orthonormalbasis des Rn, wenn zusatzlichzu den Bedingungen einer Basis gilt

xTj xk =

{0 j 6= k1 j = k

Orthogonale Matrizen:Seien x1, . . . xn eine Orthonormalbasis und U = (x1| . . . |xn), dann gilt

UT U =

xT1...

xTn

(x1| . . . |xn) =

xT1 x1 . . . xT

1 xn...

...xT

nx1 . . . xTnxn

= I

U heißt orthogonale Matrix mit der Eigenschaft U−1 = UT .

Transponieren von Matrizen:Im Reellen: A = (ajk) → AT = (akj)Im Komplexen: A = (ajk) → AH = conj(AT ), ajk ∈ C, conj(a + ib) = a− ib

Sonderfalle:Hermitsche Matrix: AH = A (komplex)symmetrische Matrix: AT = A (reell)orthogonale Matrix: A−1 = AT (reell)unitare Matrix: A−1 = AH (komplex)

Normalformen von Matrizen:

• Schur-Normalform:A = QR QH︸︷︷︸

=Q−1

,

Q unitar, R =

λ1 * . . . *

0. . . . . .

......

. . . . . . *0 . . . 0 λn

obere Dreiecksmatrix

det(A) = det(QRQH) = det(R) = Πni=1λi

12

Symmetrische Matrix:

A = AT ⇒ RT = R ⇒ R =

λ1 0. . .

0 λn

• Jordan-Normalform:

A = SJS−1, J =

J1 0. . .

0 Jn

, J1 =

λ1 1 0

. . . . . .. . . 1

0 λ1

Die Jordan-Normalform ist jedoch numerisch uninteressant.

1.3.3 Lineare Gleichungssysteme

a11x1+ . . . +a1mxm = b1...

...an1x1+ . . . +anmxm = bn

Ax = b

Fragestellung:

• Liegt b im Bildraum der Abbildung A?

• Erhoht sich der Rang der Matrix A, wenn b dazugenommen wird?

Es gilt: Ax = b losbar ⇔ rang(A) = rang(A|b) ⇒ b ∈ Bild(A)Solche Gleichungssysteme konnen mit der Gauss-Elimination gelost werden.

Dreieckssysteme: a11 . . . a1n

. . ....

0 ann

x1

...xn

=

b1...bn

Dreieckssysteme konnen wie folgt mit Kosten O(n2) gelost werden:annxn = bn ⇒ xn = bn

ann

an−1,n−1xn−1 + an−1,nxn = bn−1

⇒ xn−1 = bn−1−an−1,nxn

an−1,n−1

etc.

13

Algorithmus fur die allgemeine Gauss-Elimination mit Umwandlung in Drei-ecksform und Spaltenpivotsuche:

FOR k = 1: n - 1

alpha = |a(k, k)|; j = k;

FOR s = k + 1: n

IF |a(s, k)| > alpha THEN

alpha = |a(s, k)|; j = s;

ENDIF

END

# Pivotelement a(j, j) und Pivotzeile ist j

FOR i = k: n

alpha = a(k, i); a(k, i) = a(j, i); a(j, i) = alpha;

END

alpha = b(j); b(j) = b(k); b(k) = alpha;

# Eliminationsschritt

FOR s = k + 1: n

I(s, k) = a(s, k) / a(k, k);

b(s) = b(s) - I(s, k)b(k);

FOR i = k + 1: n

a(s, i) = a(s, i) - I(s, k)a(k, i);

END

END

END

Interpolation als LU -Zerlegung ohne Pivotsuche:A = LU , L ist untere, U obere Dreiecksmatrix

b = Ax = (LU)x = L (Ux)︸︷︷︸=y

= Ly

Die Gleichungen Ly = b und Ux = y konnen jeweils mit Kosten O(n2) gelostwerden, die Gauss-Elimination selbst benotigt O(n3).

QR-Zerlegung:A = QR, Q orthogonal, R obere Dreiecksmatrix, Q ist Givens-Rotation

b = Ax = (QR)x = Q(Rx) ⇔ QT b = Rx

14

Da QT b und Rx in O(n2) gelost werden konnen, kann die Gleichung in O(n2)gelost werden.

Least Squares Problem:

minx ‖ Ax− b ‖2 → AT Ax = AT xminx ‖ Ax− b ‖2 = minx ‖ (QR)x− b ‖2

= minx ‖ QT (QRx− b) ‖2

= minx ‖ Rx− QT b︸︷︷︸b′

‖2

= minx ‖ (R′

0)x− ( b′1

b′2) ‖

2

Es ergibt sich: R′x = b′1, b′2 beeinflusst also den Fehler.

Erweiterte Betrachtung der Matrixnormen:

‖ A ‖22 = maxx 6=0

‖Ax‖22‖x‖22

= maxx 6=0xT (AT A)x

xT x=

= λmax(AT A) = σ2

max(A)

‖ A ‖2 =√

λmax(AT A) = σmax(A)

‖ A ‖2F =

∑i,j a2

ij = ‖ UΣV T ‖2

F = ‖ Σ ‖2F =

= ‖

σ1 0. . .

σk

0. . .

0 0

‖2

F , σk: Singularwerte von A

=∑k

j=1 σ2j

‖ A ‖F =√∑k

j=1 σ2j

Singularwertzerlegung:A = UΣV T , U , V T orthogonal, Σ Diagonalmatrix

AT A = (V ΣT UT )(UΣV T )= V ΣT UT U︸︷︷︸

=I

ΣV T

= V ΣT ΣV T

⇒ AAT = U(ΣΣT )UT

15

Fehlerabschatzung:

minx ‖ Ax− b ‖2 = minx ‖ UT (UΣV T x− b) ‖2

= min ‖ Σ (V T x)︸︷︷︸y

−UT b︸︷︷︸b′

‖2

= miny ‖ Σy − b′ ‖2

= miny ‖

σ1y1...

σkyk

0...0

−

b′1...

b′k...

b′n

‖2

Es ergibt sich yj =b′jσj

, 1 ≤ j ≤ k und x = V y.

1.4 Graphen und dunnbesetzte Matrizen

Ein Graph ist ein Tupel G = {E, K} mit Ecken p1, . . . pn ∈ E und Kanten(pi, pj) ∈ K.Im gerichteten Graph gilt: (p1, p2) 6= (p2, p1), im ungerichteten Graphen:(p1, p2) = (p2, p1).

Wege im Graphen:

p1 // p2 oder p11 // p2 ⇔ (p1, p2) ∈ K

p12 // p2 ⇔ (p1, pk) ∈ K und (pk, p2) ∈ K fur passendes k

p1k // p2 oder p1

* // p2 ⇔ es gibt einen Weg von p1 nach p2

Transitiver Abschluss von G:

Beispiel 1.5

Anfang: p5

p1 // p2 //

OO

p3 // p4jj

16

1. Schritt: p5

p1 // 66

>>}}}}}}}}p2 //

OO

((p3 //�� p4^^

2. Schritt: p5

p1 // 66

>>}}}}}}}} ��p2 //

||

OO

((p3 //bb�� p4^^

Abschluß: p5

p1 oo //OO

��

hh

((PPPPPPPPPPPPPPP

>>}}}}}}}}p4OO

��

66

vvnnnnnnnnnnnnnnn

``AAAAAAAA

p2 oo //

11

p3

nn

K4 = {p1, p2, p3, p4} ist ein vollstandiger Teilgraph (Clique) mit vier Kno-ten.

Ein Untergraph ist ein Graph, der aus Teilmengen von E und K besteht.Die Anzahl der von einer Ecke auslaufenden Kanten wird als Außengrad, dieZahl der einlaufenden Kanten als Innengrad der Ecke bezeichnet.

Die Ecke pm heißt von pi aus erreichbar, falls pi* // pm , ein Graph heißt

stark zusammenhangend, wenn ∀i 6= j : pi* // pj .

Bewerteter Graph:

pi(ai)bij // pj(aj) mit Ecken pi, pj und Gewichten ai, aj und bij

17

Quotientengraph:Der Quotientengraph entsteht durch das Bilden von disjunkten Teilmengeny1 . . . yk von p1 . . . pn. Der Quotientengraph enthalt die Kante (yi, yj), wennp ∈ yi, z ∈ yj und (p, z) ∈ K.

p4

��

p1 // p2

>>}}}}}}}}

AAA

AAAA

A

p3

; y1 // y2 mit y1 = {p1, p2}, y2 = {p3, p4}

Adjazenzmatrix:Die Adjazenzmatrix ist eine |E| × |E|-Matrix mit

aij =

{1 falls (pi, pj) ∈ K0 sonst

Inzidenzmatrix:Die Inzidenzmatrix ist eine |K| × |E|-Matrix mit

aij =

1 falls die Kante mi in pj beginnt−1 falls die Kante mi in pj endet0 sonst

Beispiel 1.6

p11// p2

2

__p34 ""

5

__p4

3bb

18

Adjazenz:

p1 p2 p3 p4

p1 0 1 0 0p2 0 0 0 1p3 1 0 0 1p4 0 0 1 0

, Inzidenz:

p1 p2 p3 p4

m1 1 −1 0 0m2 0 1 0 −1m3 0 0 −1 1m4 0 0 1 −1m5 −1 0 1 0

2

Moglichkeiten der Speicherung von Adjazenzmatrizen:

• bei der Speicherung als Bitvektor wird jede Zeile der Matrix als Binardar-stellung einer int-Zahl aufgefasst

• Speicherung von Zeigern auf die 1-Elemente:

p1 : 2p2 : 4p3 : 1 → 4p4 : 3

• Speicherung als Matrix:

p1 2 0p2 4 0p3 1 4p4 3 0

• Speicherung in Cliquen: p6

AAAA

AAAA

p7

}}}}

}}}}

p4 p5

p3

AAAA

AAAA

p1 p2

setzt sich zusammen aus den Cliquen {p1, p2, p3}, {p4, p5, p6, p7}, {p3, p6}und {p2, p7}.

19

Allgemeine Vorgehensweise:Man startet mit einer beliebigen Matrix A und definiert zu A einen Gra-phen G(A). Gemaß der Definition einer Adjazenzmatrix gilt: (pi, pj) ∈ K inG(A) ⇔ aij 6= 0

einfacher:

1, 1 0 3, 5

7, 6√

2 07 0 0

Adjazenzmatrix//

1 0 11 1 01 0 0

Beispiel 1.7

A =

* * 0 0 0* * 0 0 00 0 * * *0 0 * * *0 0 * * *

Adjazenzmatrix//

1 1 0 0 01 1 0 0 00 0 1 1 10 0 1 1 10 0 1 1 1

Daraus ergibt sich folgender Adjazenz-Graph G(A) (ohne Loops):

p1 p2 p3 p4 p5

2

Umbenennung :Durch Umbenennung der Knoten im Adjazenzgraphen entsteht eine Permu-tation der ursprunglichen Matrix

A =

* * 00 * 0* 0 *

G(A) // p1 // p2 p3�� (ohne Loops)

Umbenennung

��

A′ =

* 0 0* * 00 * *

A′oo p2 // p1 p3��

Dies lasst sich aquivalent durch das Matrixprodukt P T AP mit einer passen-den Permutationsmatrix P ausdrucken: 0 1 0

1 0 00 0 1

* * 00 * 0* 0 *

0 1 01 0 00 0 1

=

* 0 0* * 00 * *

20

Reduzible Matrizen:Eine Matrix A heißt reduzibel, wenn es eine Permutation P gibt, so daß gilt

P T AP =

(C 0E D

)Angewandt auf ein lineare Gleichungssystem ergibt sich:

(P T AP )x = b ⇔(

C 0E D

)(x1

x2

)=

(b1

b2

)Fur die Losung ergibt sich:

Cx1 = b1, Ex1 + Dx2 = b2 ⇒ Dx2 = b2 − Ex1

Der Vorteil von reduziblen Matrizen besteht also darin, daß das Gesamt-Gleichungssystem Ax = b auf (in diesem Fall zwei) kleinere lineare Glei-chungssysteme reduziert werden kann.

Allgemein gilt:

A reduzibel ⇒ P T AP =

C1 0. . .

* Ck

, Cj ∈ Rn×n

Die aus der Reduktion entstehenden Matrizen Cj sind dann irreduzibel.

Satz 1.3 Grapheneigenschaften irreduzibler Matrizen

A irreduzibel ⇔ G(A) stark zusammenhangend

Beweis:

”⇒“:

A reduzibel ⇒ P T AP =

(C 0E D

)Als Adjazenzmatrix interpretiert: es gibt keine Kante, die von der erstenGruppe der Indizes in die zweite fuhrt. G(A) zerfallt also in zwei Zusammen-hangskomponenten und ist daher nicht stark zusammenhangend.

21

”⇐“:

Wenn G(A) nicht stark zusammenhangend ist, so gibt es zwei Ecken pi und

pj mit pi* // pj , das ist aber nicht moglich. Werden alle Indizes pr mit

pi* // pr in einer Indexmenge I1 und alle anderen Indizes in der Index-

menge I2 zusammengefasst, ergibt sich nach der symmetrischen PermutationI1 → {1 . . . k}, I2 → {k + 1 . . . n} das folgende Matrixprodukt:

P T AP =

( I2

I1 C 0E D

)A ist also reduzibel.

2

Beispiel 1.8

A =

* 0 0 * 0 0 *0 * 0 0 * 0 00 0 * 0 0 * 0* 0 0 * 0 0 *0 * 0 0 * 0 00 0 * 0 0 * 0* 0 0 * 0 0 *

G(A) : p1 p2 p3 p4 p5 p6 p7 (ohne Loops)

G(A) enthalt folgende Zusammenhangskomponenten:z1 = {p1, p4, p7}, z2 = {p2, p5}, z3 = {p3, p6}

22

Durch Umnummerierung erhalten die Zusammenhangskomponenten aufein-anderfolgende Nummern:

1 → 14 → 27 → 3

z12 → 45 → 5

}z2

6 → 63 → 7

}z3

Durch die Umnummerierung entstehen neue Zusammenhangskomponenten:

z′1 = {p1, p2, p3}, z′2 = {p4, p5}, z′3 = {p6, p7}

Aus diesen neuen Zusammenhangskomponenten ergibt sich folgende Permu-tationsmatrix P :

P =

1 00 0 1 00 0 0 1

0 1 0 0 01 0

1 00 0 1 0 0 0 0

Nach der symmetrischen Permutation P T AP entsteht eine wesentlich uber-sichtlichere und leichter zu losende Form der ursprunglichen Matrix A:

P T AP =

* * ** * ** * *

* ** *

* ** *

Das Gesamtsystem zerfallt also in drei kleinere, irreduzible Teilsysteme, dieden Zusammenhangskomponenten des Graphen entsprechen, nachdem dieseSchritte durchlaufen wurden:

1. Ubergang von A zu G(A)

2. Auflosen von G(A) in Zusammenhangskomponenten

3. Umnummerierung der Zusammenhangskomponenten

4. Symmetrische Permutation

2

23

2 Dunnbesetzte Vektoren und Matrizen

2.1 Dunnbesetzte Vektoren

Dunnbesetzte Vektoren konnen statisch oder dynamisch gespeichert werden.Bei der statischen Speicherung wird der Vektor in einem fest vorgegebenenSpeicherbereich abgelegt, mit dem Nachteil, das dieser Bereich zu klein wer-den kann. Bei der dynamischen Speicherung wird der Speicher dynamischallokiert, dafur ensteht hier Overhead fur Verwaltung und Allokation. Nebender statischen und dynamischen Speicherung gibt es die folgenden Moglich-keiten eine dunnbesetzten Vektor komprimiert oder unkomprimiert zu spei-chern.

Erste Moglichkeit :Der dunnbesetzte Vektor wird vollbesetzt als Array voller Lange n mit expli-ziten Nullen gespeichert. Diese Variante ist relativ einfach in der Handhabungund jeder Eintrag ist direkt ansprechbar, auerdem wird mit O(n) relativ we-nig Speicher verbraucht.

Zweite Moglichkeit :Es werden nur (Wert, Index)-Paare gespeichert. Bei dieser Variante wirdhaufig ein Array V AL fur die Werte und ein Array IND fur die Indizes ver-wendet.

Beispiel 2.1

x ∈ Rn, Dimension nV AL = (1.0, 2.1, 0.3)IND = (1, 2, 5)NZ = 3 (Anzahl der Nonzeros, Lange von V AL und IND)

2

Diese beiden Varianten der Speicherung konnen durch zwei einfache Algo-rithmen ineinander ubergefuhrt werden:

• scatter:Es wird zunachst ein Feld B der Lange n mit Nulleintragen erzeugt.Dann wird IND durchlaufen und die Werte aus V AL an den durchIND angegebenen Positionen in B mit Kosten O(NZ) eingefugt.

24

• gather:Es wird der volle Vektor durchlaufen und die Nonzero-Eintrage werdenin V AL und IND gespeichert.

2.1.1 Das SDOT-Problem

Das SDOT-Problem ist die Berechnung des inneren Prodults zweier Vektorenx und y, also xT y = Σn

i=1xiyi. In Abhangigkeit von der Form, in der x und yvorliegen, gibt es die folgenden Algorithmen.

Erster Algorithmus :x ist in V AL, IND und NZ gepackt, y ist voller Vektor der Lange n.

prod = 0;

FOR k = 1 : NZ

prod += VAL(k) * y(IND(k));

END

Dieser Algorithmus hat Kosten von O(NZ).

Zweiter Algorithmus :x und y liegen gepackt vor: x = (V ALx, INDx, NZx), y = (V ALy, INDy,NZy).

• Erste Moglichkeit: scattern von y, dann Anwenden des ersten Algorith-mus (SDOT)

• Zweite Moglichkeit unter der Annahme, daß INDx und INDy aufstei-gend sortiert sind:

Kx = Ky = 1;

DO UNTIL Kx > NZx OR Ky > NZy

INDx(Kx) > INDy(Ky): Ky++;

< Kx++;

= prod += VALx(Kx) * VALy(Ky);

Kx++; Ky++;

END

Kosten O(NZx + NZy).

25

Dritter Algorithmus :x und y sind gepackt und unsortiert.

prod = 0;

FOR k = 1 : NZx

FOR j = 1 : NZy

IF INDx(k) == INDy(j) THEN

prod += VALx(k) * VALy(j);

ENDIF

END

END

Kosten: O(NZx ·NZy).

2.1.2 Das SAXPY-Problem

Das SAXPY-Problem ist die Berechnung der Summe x = x+ ay zweier Vek-toren x und y. Wie beim SDOT-Problem gibt es auch hier in Abhangigkeitvom Speicherformat verschiedene Algorithmen.

Erster Algorithmus :x und y liegen gepackt vor.

# entpacke y in vollen Vektor w

w = scatter(y);

# durchlaufe x und pruefe w(INDx(i)) fuer jeden Eintrag VALx(i)

FOR i = 1: NZx

IF w(INDx(i)) != 0 THEN

VALx(i) += a * w(INDx(i));

w(INDx(i)) = 0;

ENDIF

END

26

# durchlaufe y und pruefe w(INDy(j)) fuer jeden Eintrag VALy(j)

FOR j = 1: NZy

IF w(INDy(j)) != 0 THEN

NZx++;

VALx(NZx) = a * w(INDy(j));

INDx(NZx) = INDy(j);

w(INDy(j)) = 0;

ENDIF

END

Kosten: O(NZx + NZy).

Zweiter Algorithmus :x und y sind sortiert und gepackt.Der Algorithmus verlauft analog zum ersten Fall.Kosten: O(NZx + NZy).

Dritter Algorithmus :x und y sind gepackt und unsortiert.

durchlaufe INDx:durchlaufe INDy:

gleicher Index: xi += ayi

durchlaufe INDy:durchlaufe INDx:

y-Index nicht vorhanden in x: xi = ayi

Kosten: O(NZx ·NZy).

2.2 Speicherung von dunnbesetzten Matrizen

Der folgende Abschnitt zeigt einige wichtige Verfahren, eine n×n-Matrix derForm

A =

1 0 0 2 03 4 0 5 06 0 7 8 90 0 10 11 00 0 0 0 12

27

zu speichern. Abhangig von der Gestalt der Matrix sind die einzelnen Vefah-ren dabei unterschiedlich effizient.

2.2.1 Koordinatenform

Bei diesem Verfahren werden alle Nonzeros der Matrix unkomprimiert ge-speichert. Dazu werden drei Vektoren verwendet: der Vektor AA enthalt dieMatrixelemente, der Vektor JR die Zeilen- und der Vektor JC die Spaltenin-dizes:

Values: AA 12 9 7 5 1 2 11 3 6 4 8 10Row-Indices: JR 5 3 3 2 1 1 4 2 3 2 3 4Column-Indices: JC 5 5 3 4 1 4 4 1 1 2 4 3

Die Lange der Vektoren betragt nnz(A), nnz(A) ist die Anzahl der Non-zeros in der Matrix.Die Speicherung in diesem Format hat Kosten O(3 ·nnz(A)), es entsteht alsoRedundanz.

Algorithmus fur das Matrix/Vektor-Produkt c = Ab:

c = 0;

FOR j = 1: NNZ(A)

c(JR(j)) += AA(j) * b(JC(j));

END

Dieser Algorithmus hat Kosten O(NNZ(A)). Ein Nachteil dieses Verfahrensist die hohere Laufzeit bedingt durch die indirekte Adressierung in c und b.

2.2.2 Compressed Sparse Row Format (CSR)

Dieses Verfahren speichert Nonzeros der Matrix zeilenweise in den VektorenAA fur die Werte und JA fur den Spaltenindex:

Values AA 2 1 5 3 4 6 7 8 9 10 11 12Column-Indices JA 4 1 4 1 2 1 3 4 5 3 4 5

Ein Querstrich bezeichnet jeweils den Beginn einer Zeile.

28

Zur Adressierung der Zeilen wird ein zusatzlicher Vektor IA verwendet, indem an der i. Position ein Zeiger auf den Beginn der i. Zeile gespeichert ist:

IA 1 3 6 10 12 13Fur die Lange von IA ergibt sich |IA| = n + 1, dieses Verfahren hat alsoKosten O(2 · nnz(A) + n + 1).

Algorithmus fur c = Ab:

c = 0;

FOR i = 1: n # durchlaufe IA

FOR j = IA(i): IA(i + 1) - 1

c(i) += AA(j) * b(JA(j)); # indirekte Adressierung

END

END

Entsprechend arbeitet das Compressed Sparse Column Format (CSC). DasSpeicherschema ist hier ahnlich, der Algorithmus arbeitet allerdings anders.

2.2.3 CSR mit extrahierter Diagonale

Bei diesem Verfahren kommt der Hauptdiagonalen eine Sonderrolle zu, dadie Diagonalelemente der zu speichernden Matrix i.d.R. Nonzeros sind. Eswird folgendes Speicherschema verwendet:

AA Hauptdiagonale Nicht-Diagonale in CSRJA Zeiger auf den Beginn der i. Zeile in AA Column-Indices

Mit den Werten der Beispielmatrix:AA 1 4 7 11 12 * 2 3 5 6 8 9 10JA 7 8 10 13 14 14 4 2 4 1 4 5 3

Der in CSR fur die Spaltenindizes der Hauptdiagonalen vorgesehene Spei-cherplatz wird hier verwendet, um die Zeiger auf den jeweiligen Zeilenanfangzu speichern.Kosten: O(2(nnz(A) + 1)).

29


FOR i = 1: n

c(i) = AA(i) * b(i);

FOR j = JA(i): JA(i + 1) - 1

c(i) += AA(j) * b(JA(j));

END

END

2.2.4 Diagonalweise Speicherung

Dieses Verfahren ist v.a. fur Bandmatrizen anwendbar. Dabei werden dieeinzelnen Diagonalen mit Ganzzahlen nummeriert, so daß die Hauptdiagonalevon A die Nummer 0 erhalt:

A =

−1 0 1 2 . . .

1 0 2 0 03 4 0 5 00 6 7 0 80 0 9 0 00 0 0 11 12

Anhand dieser Nummerierung wird jetzt jede Diagonale von A Spalte einerneuen Matrix COEF:

COEF =

* 1 23 4 56 7 89 0 *11 12 *

, IOFF = (-1, 0, 2)

Im Vektor IOFF werden die zu den Spalten von COEF gehorenden Nummernder Diagonalen gespeichert.

Fur Bandmatrizen entsteht bei diesem Verfahren nahezu kein Overhead anuberflussigem Speicherplatz.Kosten: O(n · nd), nd ist die Zahl der Diagonalen.

Allgemeiner fur nicht-Bandmatrizen:Fur nicht Band-Matrizen lassen sich

”kunstliche“ Diagonalen erzeugen, in-

dem zeilenweise jeweils das am weitesten links stehende Element, das nochnicht gespeichert wurde an einen Vektor angefugt wird. Diese Vektoren sinddann die Spalten der neuen Matrix COEF:

30

A =

1 0 2 0 03 4 0 5 00 6 7 0 80 0 9 10 00 0 0 11 12

, COEF =

1 2 03 4 56 7 89 10 011 12 0

Da es sich um keine echten Diagonalen handelt, mussen zusatzlich die Spal-tenindizes der Eintrage von COEF im Vektor JCOEF gespeichert werden:

JCOEF =

1 3 ?1 2 42 3 53 4 ?4 5 ?

Kosten: O(2n · nl), nl := Lange der langsten Zeile von A. Dieses Verfahrenist jedoch nur effizient fur nl � n.


c = 0;

FOR i = 1: n

FOR j = 1: nl

c(i) += COEF(i, j) * b(JCOEF(i, j));

END

END

2.2.5 Jagged Diagonalformat

Bei diesem Verfahren werden die Zeilen der Matrix zunachst so umgeordnet,daß die Zeilen mit den meisten Nonzeros oben sind, die Zeilen werden alsonach unten monoton kurzer.

A =

1 0 2 0 03 4 0 5 00 6 7 0 80 0 9 10 00 0 0 11 12

Nach der Umordnung der Zeilen werden analog zur bandweisen Speiche-rung fur Nicht-Band-Matrizen

”kunstliche“ Diagonalen erzeugt, die monoton

kurzer werden. Im Beispiel ergeben sich die folgenden Jagged-Diagonalen:

31

A =

3

====

=== 4

====

=== 0 5

BBBB

BBBB

B 0

0 6

��

��

7 0 8

1

NNNNNNNNNNNNNN 0 2

@@@@

@@@@ 0 0

0 0 9

@@@@

@@@@ 10

BBBB

BBBB

0

0 0 0 11 12

Diese Jagged-Diagonalen werden in den Vektoren AA fur die Werte und JAfur die Spaltenindizes gespeichert:

Values AA 3 6 1 9 11 4 7 2 10 12 5 8Column-Indices JA 1 2 1 3 4 2 3 3 4 5 4 5

Zusatzlich wird ein Vektor JDIAG verwendet, der an der i. Stelle einen Zei-ger auf den Beginn der i. Jagged-Diagonale enthalt:

JDIAG: 1 6 11 13Lange von JDIAG: |JDIAG| = Anzahl der Jagged-Diagonalen︸︷︷︸

NDIAG

+ 1.

Kosten: O(2nnz(A) + |JDIAG|).


c = 0;

FOR i = 1: NDIAG

FOR j = JDIAG(i): JDIAG(i + 1) - 1

c(i) += AA(j) * b(JA(j));

END

END

2.2.6 Zusammenfassung

Allgemein lasst sich also festestellen, daß es Speicherverfahren gibt die zei-lenweise, spaltenweise oder auch diagonalweise arbeiten.

32

Bemerkung:

• Zu jeder n× n-Matrix existiert ein Graph G(A)

• Zu jedem linearen Gleichungssystem existiert ein bewerteter Graph:

Ax = b: p1, b1

a13

!!

a12 99 p2 p3 p4, b4

a42

``

Dieser Graph enthalt Kanten, die mit dem Wert des Matrixelementsbewertet sind und Ecken, die mit dem Element der rechten Seite ge-wichtet sind.

2.3 Elementare Matrixoperationen

2.3.1 Das Matrix/Vektor-Produkt

Das Matrix/Vektorprodukt kann zuruckgefuhrt werden auf ein inneres Pro-dukt zwischen Vektoren (SDOT) oder auf die Addition von Vektoren (SAX-PY).

Erste Moglichkeit :

y = Ax =

a1...

an

︸︷︷︸

A zeilenweise

·x =

a1x...

anx

: n-mal SDOT

Ist x vollbesetzt, enstehen n-mal Kosten O(p) fur das Produkt eines dunn-besetzten mit einem vollbesetzten Vektor, also insgesamt O(np).

Zweite Moglichkeit :

y = Ax = (α1 . . . αn)︸︷︷︸A spaltenweise

·

x1...

xn

= x1α1 + . . . + xnαn

Fur n Additionen eines vollbesetzten mit einem dunnbesetzten Vektor enste-hen jeweils Kosten O(p), somit ergeben sich Gesamtkosten von O(np).

33

In seltenen Fallen kann es vorkommen, daß x selbst dunnbesetzt ist. Auchin diesem Fall gibt es eine zeilenweise und eine spaltenweise Variante derProblemlosung.

Spaltenweise Variante:

y =∑

j∈J xjαj, xj 6= 0 fur j ∈ J , |J | ≤ p.

Hier mussen also wenige SAXPYs zwischen dunnbesetzten Vektoren berech-net werden, fur jedes dieser SAXPYs ergeben sich Kosten O(p), also Gesamt-kosten O(p2).

Zeilenweise Variante:

y =

a1x...

anx

, eigentlich n-mal SDOT

Durch unnotiges Berechnen von Nullen (da die meisten aix = 0) entste-hen Kosten O(n). Dieser unnotige Aufwand lasst sich durch Betrachten desAdjazenzgraphen vermeiden, indem anhand des Patterns eine

”Vorhersage“

getroffen wird, welche aix 6= 0 sind. Gegeben sei das folgende Produkt:

y = Ax =

j1 j2

* *...

...* *

·

0

j1 *0

j2 *0

Die Indexmengen Ij1 und Ij2 bezeichnen die Nichtnulleintrage der Spalten j1

und j2 von A. Ist der Index j /∈ Ij1 ∪ Ij2 , so ist yj = 0, es mussen also nurdie yj = ajx fur die Indizes j mit j ∈ Ij1 ∪ Ij2 berechnet werden. Fur dasPattern und den zugehorigen Graphen ergibt sich:

0 * 0 *0 0 * 0* 0 0 *0 0 * 0

*00*

=

*0*0

G(A) : p1

**// p2 // p3hh p4

Fur die Berechnung sind die Ecken pi relevant, an denen xi 6= 0 ist, hier alsop1 und p4. Im resultierenden Vektor y ist yj 6= 0, wenn eine Kante von pj

nach pi fuhrt.

34

2.3.2 Das Matrix/Matrix-Produkt

Standardverfahren:Das zu berechnende Matrix/Matrix-Produkt liegt in der normalen innerenProduktform vor:

C = AB =

ai * . . . *

·

bj

*...*

In dieser Form lasst sich das Problem mit Hilfe des Standardverfahrens losen.Dabei ist die Matrixkomponente cij = aibj, also ein SDOT zwischen dendunnbesetzten Vektoren ai und bj.Es ergibt sich jedoch das Problem, daß die meisten der cij = 0 sind, zusatzlichentstehen Kosten O(n2). Hier lasst sich analog zum Matrix/Vektor-Produktdas Pattern der resultierenden Matrix C aus den Graphen der beiden Fak-toren A und B ablesen. Es gilt: cij = Σn

k=1aikbkj 6= 0, wenn

• fur ein k aik 6= 0 und bik 6= 0

• ein k existiert, so daß piG(A) // pk

G(B) // pj

Außere Produktform von AB:Als bessere Losung bietet sich die Berechnung des Produkts in der folgendenDarstellung an:

C = AB =(

α1 . . . αn

) b1

...bn

=(α1

10...0

T

+ . . . + αn

0...01

T)(

10...0

b1 + . . . +

0...01

bn

)= α1b1 + . . . + αnbn

Die aus den Produkten αibi entstehenden Matrizen haben jeweils Rang 1.

Aus der außeren Multiplikation dunnbesetzter Vektoren entsteht also einedunnbesetzte Matrix mit ≤ p Elementen. Dazu sind n Additionen einer Ma-trix mit einer dunnbesetzten Matrix notig, in Abhangigkeit vom Speicherfor-mat konnen Kosten O(np2) erreicht werden.

35

Beispiel 2.2 (0 ** *

)(* *0 *

)=

(0 ** *

)Als Kombination der beiden Adjazenzgraphen G(A) und G(B) ergibt sich(ohne Loops):

p1""r _ Lp2//oo

G(A)___ G(B)

Wie oben beschrieben, lasst sich die Besetztheit von AB aus den Adjazenz-graphen von A und B ableiten. Am Beispiel c11:

c11 6= 0 ⇔ ∃k: p1G(A) // pk

G(B) //___ p1

Da aber ein kein Weg der Lange 2 existiert, dessen erster Teil in G(A) unddessen zweiter Teil in G(B) liegt, es also kein passendes k gibt, ist c11 = 0

2

Im Spezialfall B = A erhalt man das Pattern von C = AB = A2 analogdurch das Hinzunehmen aller Wege der Lange 2, wobei in diesem Fall alleTeile des Weges in G(A) liegen, da G(A) = G(B).Dieses Verfahren lasst sich fur beliebige Potenzen Ai bis zur Stagnation beiAk forsetzen, G(Ak) ist dann die transitive Hulle von G(A).

2.4 Algorithmen

2.4.1 Matrix/Vektor-Produkt

FOR i = 1: nFOR j = 1: n

ci += aijbj

END

SDOT

END

Der Algorithmus arbeitet also nach dem folgenden Schema:

c = Ab =

a1

...an

b1

...bn

=

a1b...

anb

36

2.4.2 Matrix/Vektor-Produkt mit vertauschtem i und j

FOR j = 1: nFOR i = 1: n

ci+ = aijbj

ENDEND

GAXPY (Folge von SAXPYs mit gleichem c)

Der Algorithmus arbeitet nach folgendem Schema:

c = Ab =(

a1 . . . an

) b1...bn

= b1a1 + . . . + bnan

2.4.3 Matrix/Matrix-Produkt

c∗j =∑n

k=1a∗kbkj

FOR i = 1: nFOR j = 1: n

FOR k = 1: ncij += aik︸︷︷︸

zeilenweise

· bkj︸︷︷︸spaltenweise

END

SDOT

ENDEND

2.4.4 Matrix/Matrix-Produkt mit vertauschtem k und j

ci∗ =∑n

k=1aikbk∗

FOR i = 1: nFOR k = 1: n

FOR j = 1: ncij += aikbkj

END

SAXPY

END

GAXPY

END

Es werden also n-mal dunnbesetzte Vektoren mitKosten O(p2) addiert, waszu Gesamtkosten O(np2) fuhrt.

37

3 Gauss-Elimination

3.1 kij-Standardform

Beim diesem Standardverfahren der Gauss-Elimination wird eine n × n-Matrix A mit Rang(A) = n nach dem folgenden Verfahren eliminiert, bissie obere Dreiecksgestalt erreicht hat:

A =

a11

��

��

. . . a1n

.... . .

...

an1 . . . ann

;

a11 . . . a1n

0 a22

��

��

. . . a2n

......

. . ....

0 an2 . . . ann

; . . . ;

a11 * . . . *

0 a′22 . . . *

.... . .

...

0 . . . 0 a′nn

38

Nach der Termination des folgenden Algorithmus, der das oben beschriebeneVerfahren umsetzt, liegt die Matrix A in eine untere Dreiecksmatrix L undeine obere Dreiecksmatrix U zerlegt vor, es gilt also A = LU .

FOR k = 1: n - 1FOR i = k + 1: n

lik = aik

akk

END

l∗k = 1akk

a∗k

FOR i = k + 1: nFOR j = k + 1: n

aij = aij − likakj

END

SAXPY

END

kein GAXPY

END

Dieses Verfahren wird auch als”Right-Looking-Gausselimination“ bezeich-

net. L wird dabei spaltenweise und U zeilenweise berechnet, wobei 1 ≤ k < i,j ≤ n gilt. Schematisch lasst sich das Verfahren wie folgt veranschaulichen:

. . . U. . .

* * ** 2 2

L * 2 2

* aktuelle Berechnung2 upgedatet

L, U bereits berechnet

3.2 ikj-Form

FOR i = 2: nFOR k = 1: i - 1

lik = aik

akk

FOR j = k + 1: naij− = likakj

END

SAXPY

END

GAXPY zeilenweise

END

39

Bei diesem Verfahren werden sowohl L als auch U zeilenweise von oben nachunten berechnet. Schematisch dargestellt arbeitet der Algorithmus wie folgt:

. . . U. . .

L. . .

* * * * *

A

* gerade berechnet

L, U bereits berechnetA original

3.3 jki-Form

FOR j = 2: nFOR m = j: n

lm,j−1 =am,j−1

aj−1,j−1

ENDFOR k = 1: j - 1

FOR i = k + 1: naij− = likakj

ENDEND

GAXPY spaltenweise

END

Dieses Verfahren wird auch als”Left Looking Gauss-Elimination“ bezeich-

net. L und U werden dabei spaltenweise von links nach rechts berechnet.Schematisch dargestellt arbeitet das Verfahren wie folgt:

. . . U *. . . *

* A*

L *

* gerade berechnet

L, U bereits berechnetA original

3.4 Crout-Form

Ausgehend von dem Produkt LU = A werden L und U ausgehend von demfolgenden Ansatz berechnet, zur Vereinfachung seien L und U dabei jeweils3× 3-Blockmatrizen:

40

bekannt l11 0 0wird berechnet l21 l22 0ignorieren l31 l32 l33

·

bekannt wird berechnet ignorieren

u11 u12 u13

0 u22 u23

0 0 u33

=

=

l11u11 l11u12 l11u13

l21u11 l21u12 + l22u22 ** * *

!=

a11 a12 *a21 a22 ** * *

Dieses Matrixprodukt fuhrt auf die folgenden Gleichungen:

• l11u11!= a11: muß durch vorhergehende Rechnung erfullt sein

• l11u12 = a12

l21u11 = a21

l21u12 + l22u22 = a22

Die unbekannten Blockmatrizen werden dann durch die folgenden Schrittebestimmt:

1. l11u12 = a12: Losen eines Dreiecksgleichungssystems in l11

2. l21u11 = a21: Losen eines Dreiecksgleichungssystems in u11

3. l21u12 + l22u22 = a22

l22u22 = a22 − l21u12 = a′22 (l21 und u12 sind jetzt bekannt)

Das Problem ist also auf eine kleinere LU -Zerlegung von a′22 und damit einekleinere Gauss-Elimination reduziert. Schematisch lasst sich die Crout-Formder Gauss-Zerlegung wie folgt darstellen:

2 2 2

2 2 2

2 2 2

2 2 2 * * *2 2 2 *2 2 2 *

* wird berechnet2 wird verwendetRest wird nicht verwendet

3.5 Verfahren fur symmetrisch positiv definite Matri-zen

Die folgenden Verfahren sind nur fur symmetrisch positiv definite Matrizen,d.h. Matrizen der Form

A = AT = LU = LLT

anwendbar. Die Zerlegung A = LLT wird als Cholesky-Zerlegung bezeichnet.

41

3.5.1 Borderingform

Ausgehend von der Cholesky-Zerlegung A = LLT soll die i. Spalte von Aberechnet werden. li1 bis li,i−1 seien bereits bekannt, zu berechnen ist also lii: l11 0 0

.... . . 0

li−1,1 . . . li−1,i−1

li1

...li,i−1

=

ai1...

ai,i−1

Dieser Ansatz fuhrt auf die folgende Gleichung:aii = lil

Ti (Produkt der i. Zeile von L mit der i. Spalte von LT )

= Σij=1l

2ij

= Σi−1j=1l

2ij + l2i,i

⇒ lii =√

aii − Σi−1j=1l

2ij︸︷︷︸

>0

Schematisch lasst sich die Bordering-Form der Gauss-Elimination fur L wiefolgt darstellen (die Darstellung fur LT ist wegen A = LLT analog):

. . .

2. . .

2 2. . .

* * *. . .

. . .. . .

* wird berechnet2 verandertRest unverandert

3.5.2 Innere Produktform

ajj = Σjk=1l

2jk ⇒ ljj =

√ajj − Σj−1

k=1l2jk, 1 ≤ j ≤ n

aij = Σjk=1ljklik = Σj−1

k=1likljk︸︷︷︸bekannt

+lijljj

lij = 1ljj

(aij − Σj−1k=1likljk), j + 1 ≤ i ≤ n

42

Schematische arbeitsweise des Algorithmus fur L (LT analog):

. . .. . .

L. . .

2 2 2. . .

2 2 2 *. . .

2 2 2 *. . .

* wird berechnetL bereits berechnet2 verandertRest unverandert

43

4 Band- und Envelope-Methoden fur sym-

metrisch positiv definite Matrizen

4.1 Band

Definition 4.1 Kleinster Nonzero-Index, Diagonalenabstand, Bandmetho-de, Band(A)

• fi(A) := min{j|aij 6= 0}, 1 ≤ i ≤ n, ist der kleinste Spaltenindex inZeile i mit Nonzero-Eintrag

• βi(A) := i − fi(A) ist der Abstand von der Diagonalen zum kleinstenNonzero-Eintrag

• Band(A) := {(i, j)|0 < i− j ≤ β(A)}

• β(A) := maxni=1βi(A) heißt die Bandmethode

Beispiel 4.1

* * 0 * 0* * 0 0 00 0 * 0 * ** 0 0 * * 0 0

0 * * * 0 ** 0 0 * *

0 0 * * *

Zeile fi(A) βi(A) = i− fi(A)1 1 02 1 13 3 04 1 35 3 26 3 37 5 2

Der umrandete Bereich bezeichnet Band(A). Alle in Band(A) enthaltenenNullen werden explizit gespeichert, alle Nullen, die außerhalb des Bandes lie-gen jedoch nicht.

A kann also dargestellt werden durch eine Bandmatrix mit der Bandbreiteβ. Da A symmetrisch positiv definit ist, gilt fi ≤ i, da die Diagonalelemente6= 0 sind.

2

Satz 4.1 Bandstruktur bei Gauss-Elimination ohne Pivotsuche (ohne Be-weis)

Eine Gauss-Elimination ohne Pivotsuche liefert nach Terminierung A = LLT

mit Band(A) = Band(L+LT ), die Bandstruktur der ursprunglichen Matrixbleibt also erhalten.

2

44

4.2 Envelope

Definition 4.2 Envelope

Env(A) := {(i, j)|0 < i− j ≤ βi(A)} = {(i, j)|fi(A) ≤ j < i}, jede Zeile hatdabei die lokale Bandbreite βi.

Beispiel 4.2

** *

** 0 0 *

* * ** 0 0 *

0 * * *

Der umrandete Bereich bezeichnet die Envelope von A. Alle in der Enve-lope enthaltenen Nullen sind explizit gespeichert, alle, die außerhalb liegendagegen nicht.

2

Die Envelope hat die folgenden Eigenschaften:

• Env(A) ⊆ Band(A)

• |Env(A)| = Σni=1βi(A)

Satz 4.2 Envelope bei der Gauss-Elimination ohne Pivotsuche (ohne Be-weis)

Eine Gauss-Elimination ohne Pivot-Suche liefert nach Terminierung A =LLT mit Env(A) = Env(L + LT ), die Envelope der ursprunglichen Matrixbleibt also erhalten.

2

45

Definition 4.3 Frontbreite, maximale Frontbreite

• ωi(A) := |{k|k > i und akl 6= 0 fur ein l ≤ i}| heißt Frontbreite derMatrix A und ist die Zahl der im i. Schritt der Gauss-Eliminationbetroffenen Zeilen

• ω(A) := maxni=1ωi(A) heißt maximale Frontbreite von A

Beispiel 4.3

* * 0 * 0* * 0 0 00 0 * 0 * ** 0 0 * * 0 0

0 * * * 0 ** 0 0 * *

0 0 * * *

Zeile ωi(A) βi(A)1 2 02 1 13 3 04 2 35 2 26 1 37 0 2

2

Den Zusammenhang zwischen den Kosten der Gauss-Elimination und derFrontbreite zeigt die folgendene Matrix beim i. Schritt der Gauss-Elimination:

* . . .

... *

��

��

0 * 0 *

0 . . .

*

0

*

Im i. Schritt entstehen also Kosten: ω2

i (A) fur die Gauss-Elimination undsomit Gesamtkosten von

∑ni=1 ω2

i (A).

46

Beispiel 4.4

A =

* . . . . . . *... * 0 . . . 0

0. . .

......

... 0* 0 . . . 0 *

Eine Standard-Gauss-Elimination ist hier ungunstig, da A nach direkterGauss-Elimination vollbesetzt ist.

Durch Permutation soll eine moglichst kleine Bandbreite erreicht werden,dazu werden die erste Zeile und die erste Spalte in die Mitte gebracht:

A′ =

* . . . *...

. . ....

. . .

* . . . * . . . *. . .

.... . .

...* . . . *

, β = n2

Durch eine weitere Permutation soll eine moglichst kleine Enevelope erreichtwerden, die erste Zeile und die erste Spalte werden dazu nach hinten bzw.unten gebracht:

A′′ =

* *. . .

...* . . . *

Env(A′′) enthalt jetzt n − 1 Elemente, in dieser Form sind Speicher undRechenzeit optimal.

2

Ist j naher an der Diagonalen als fi(A), so existiert ein k mit pi → pk inG(A), 1 ≤ k ≤ j, genau dann wenn in G(Env(A)) ein Weg von pi nach pj

existiert.A soll also so permutiert werden, daß Band(A) oder Env(A) verkleinert wer-den. Dieses Ziel lasst sich mit dem Cuthill-McKee-Algorithmus erreichen.

47

4.3 Der Cuthill-McKee-Algorithmus

Durch die symmetrische Permutation PAP T bzw. eine Umnummerierung imGraphen soll die Band- und Envelope-Struktur von A verbessert werden,die Eintrage sollen also moglichst nahe an der Diagonalen liegen, so daßpi → pj ∈ K nur fur moglichst kleines |j − i|.

Idee:

1. Starte mit einer Ecke, der die Nummer 1 zugeordnet wird

2. Untersuche alle Nachbarn von Nummer 1 Diese Nachbarn erhalten dienachsten Nummern 2, 3 . . . mit steigendem Grad, Die Nummer 2 erhaltdabei die Ecke mit dem kleinsten Grad

3. Wiederhole mit Nummer 2

Algorithmus:

1. Bestimme Startecke r und setze j1 = r und k = 1, k ist dabei die Zahlder mit neuen Nummern versehenen Indizes

2. Fur i = 1 . . . n: bestimme alle bis dahin umnummerierten Nachbarnvon ji und gib ihnen die nachsten freien Nummern jk+1, jk+2 . . . nachwachsendem Grad (also zuerst die mit wenigen Nachbarn)

3. Update k

4. Gehe zum nachsten i

Beispiel 4.5

c(3) → 6

MMMMMMMMMMb(3) → 4 j(1) → 7

e(2) → 3

MMMMMMMMMM

i(1) → 10 d(4) → 9

pppppppppppg(4) → 1

hhhhhhhhhhhhhhhhhhhhh

a(2) → 8 f(4) → 5 h(2) → 2

Der Grad der Ecke ist in Klammern angegeben.

48

Anwendung des Cuthill-McKee-Algorithmus:

1. Start mit g: j1 = g, j2 = h, j3 = e, j4 = b, j5 = f

2. Weiter mit j2 = h: alle Nachbarn schon nummeriert

3. j3 = e: j6 = c

4. j4 = b: j7 = j

5. j5 = f : j8 = a, j9 = d

6. j6 = c: nichts zu tun, keine weiteren Nachbarn

7. j7 = j: nichts zu tun

8. j8 = a: nichts zu tun

9. j9 = d: j10 = i

10. Fertig

Nach der Umnummerierung ergibt sich dieses Pattern:

g h e b f c j a d ig *h * *e * *b * * *f * * *c 0 * * *j 0 * *a 0 * *d * * * *i 0 0 0 * *

Der umrandete Bereich bezeichnet die Envelope.

Die Matrix hat Bandbreite β = 4, die Envelope enthalt 24 Elemente.

2

Es lasst sich beobachten, daß die Umkehrung des Cuthill-McKee-Verfahrensmeist bessere Resultate mit kleinerer Enevelope liefert. Diese Erkenntnisfuhrt auf das Reverse Cuthill-McKee-Verfahren.

49

4.4 Der Reverse Cuthill-McKee-Algorithmus

Beim Reverse Cuthill-McKee-Algorithmus wird die Nummerierung j1 . . . jn

des ursprunglichen Algorithmus umgekehrt, so daß i1 = jn, . . . in = j1. DiesesVorgehen beruht auf der Erkenntnis, daß die Endecke des Cuthill-McKee-Algorithmus auch eine gunstige Wahl fur die Startecke sein kann. Nach derAnwendung des Reverse Cuthill-McKee-Verfahrens ergibt sich mit dem obi-gen Beispiel folgendes Pattern:

i d a j c f b e h gi *d * *a * *j *c * *f * * *b * * *e * * *h * *g * * * * *

Der umrandete Bereich bezeichnet die Envelope.

Die Matrix hat die Bandbreite β = 4, die Envelope enthalt 22 Elemente. Un-ter der Voraussetzung, daß A irreduzibel bzw. G(A) zusammenhangend ist,hat dieser Algorithmus Kosten O(n ·maximaler Eckengrad). Nicht eindeutigbeantwortet bleibt jedoch die Frage nach der Wahl des Startknotens.

Definition 4.4 Abstand, Exzentrizitat, Durchmesser, periphere Ecken

• d(pi, pj) heißt Abstand zwischen den Knoten pi und pj, also die Langedes kurzesten Weges von pi nach pj

• l(pi) := maxj {d(pi, pj)} heißt Exzentrizitat von pi

• δ(G) := maxi l(pi) = maxi,j{d(pi, pj)} heißt Durchmesser von G

• pi heißt peripher, wenn l(pi) = δ(pi)

Anhand dieser Definition bietet sich also die Wahl einer peripheren Ecke alsStartecke an.

50

Die Reihenfolge im Cuthill-McKee-Algorithmus durch kann durch Level-Sets,also Mengen von Ecken mit gleichem Abstand zur Startecke, beeinflusst wer-den:

2

��

��

5 8

BBBB

BBBB

1

====

=== 3 6 9 11

||||

||||

4 7 10

In diesem Beispiel ist d(p1, p11) = 4 und somit maximal, wahle daher 1 alsStartecke. Damit ergeben sich fur diesen Graphen die folgenden Levelsets:

S1 = {1}, S2 = {2, 3, 4}, S3 = {5, 6, 7}, S4 = {8, 9, 10}, S5 = {11}

51

5 Gauss-Elimination und der Eliminations-

graph

Der folgende Abschnitt beschrankt sich auf die Behandlung symmetrisch po-sitiv definiter Matrizen bzw. ungerichteter Graphen. Sei S(A) die Menge derNonzero-Eintrage von A:

S(A) := {(i, j)|aij 6= 0}

Liegt A in der Cholesky-Zerlegung A = LLT vor, so ergibt sich:

S(L) = S(L + LT ) = {(i, j)|lij 6= 0 oder lji 6= 0}

Außerdem wird folgendes angenommen:

• S(A) ⊆ S(L)

• fill(A) = S(L) \ S(A)

Durch die Anwendung der Right-Looking Gauss-Elimination wird der be-arbeitete Bereich der Matrix A in jedem Schritt um eine Zeile und Spalteverkurzt, aus der n × n-Matrix A0 = A wird also nach dem ersten Elimina-tionsschritt die zu bearbeitende (n− 1)× (n− 1)-Matrix A1 usw.Beim Eliminieren der ersten Spalte entstehen in den betroffenen Zeilen neueEintrage an den Stellen, an denen die erste Zeile Nonzero-Eintrage hat. DerGraph von A1 ensteht nach dem folgenden Algorithmus:

1. Entferne Ecke p1

2. Fuge neue Kanten hinzu, so daß die ursprunglich mit p1 verbundenenEcken eine Clique bilden

Es ergibt sich eine Folge von Eliminationsgraphen:

• G(A0) = G(A) mit den Ecken p1 . . . pn und der Kantenmenge K0

• G(A1) mit den Ecken p2 . . . pn und der Kantenmenge K1

• . . .

• G(An−1) mit der Ecke pn und der Kantenmenge Kn−1

• G(An) = G(L + LT ) mit der Kantenmenge KL ∪ {K0, . . . Kn−1}

52

Satz 5.1

pi → pj ∈ KL zu S(L) ⇔pi → pj ∈ KA oder ∃k < min{i, j}: pi → pk → pj ∈ KL

Beweisidee: im k. Schritt der Gauss-Elimination wird aij 6= 0 erzeugt.

2

Definition 5.1 Reachable Set

Sei S die Eckenmenge S = {pi1 . . . pik}, dann ist der Reachable Set reach(pj, S)wie folgt definiert:

reach(pj, S) := {pk /∈ S | pk ist von pj aus uber S erreichbar}

Beispiel 5.1

a S4 S3 y S1

b c S2

d

S = {S1, S2, S3, S4}, reach(y, S) = {b, a, c}.

2

Satz 5.2 Kantenmenge von G(L + LT )

Wird eine Gauss-Elimination ohne Pivotsuche in der Reihenfolge 1, 2 . . .ndurchgefuhrt, so gilt fur KL:

KL = {pi → pj | pj ∈ reach(pi, {p1 . . . pi−1})}

53

Beweisidee:nach i − 1 Eliminationsschritten sind alle Nachbarn von p1 . . . pi−1 in einerClique, beim nachsten Eliminationsschritt mit pj entsteht der neue Eintragaij 6= 0 und damit die Kante pi → pj.

2

Beispiel 5.2

1 2 3

4

6 5

;

2

====

=== 3

4

6 5

;

4 3

��

��

6 5

;

4

��

��

====

===

6 5

; 6 5 ; 6

Fur KL ergibt sich:

1 2�

��

��

3

��

��

�

��

�t

lf

4

��

��

==

==

6 5

bestehende Kante

___ eingefugte Kante

p2 → p6: p6 ∈ reach(p2, {p1})p4 → p6: p6 ∈ reach(p4, {p1, p2, p3})

2

Definition 5.2 Quotientengraph

Sei G = (E, K) und S die Menge der bereits im Verlauf der Gauss-Eliminationeliminierten Ecken. Dann Lasst sich S in Zusammenhangskomponenten c1 . . . ck

aufteilen, so daß gilt S = c1 ∪ . . . ∪ ck.

54

Sei weiterhin C(S) = {c1, . . . ck} die Menge der Zusammenhangskomponentenvon S und C′(S) =

⋃ki=1ci∪{{y} | y ∈ E \S}, dann ist der Quotientengraph

G/C′(S) der Graph, der durch die Verschmelzung der Zusammenhangskom-ponenten c1 . . . ck zu Einzelknoten entsteht.

2

Beispiel 5.3

Sei S = {1, 2, 3, 4, 5} und der folgende Graph gegeben:

7 4 2 5 1

6 8 3

9

Seien c1 = {1, 2, 4, 5} und c2 = {3} die Zusammenhangskomponenten, soergibt sich fur G/C′(S):

7 c1

��

��

????

????

6 8 c2

9

Dieses Verfahren fuhrt schließlich auf KL.

2

Aus dem Graphen entsteht also der Quotientengraph und aus diesem wie-derum nach folgendem Algorithmus der Eliminationsgraph:

1. Entferne die Superknoten in G/C′(Si)

2. Verbinde alle Nachbarn dieser Superknoten vollstandig

55

Nach der Elimination von a11 . . . a55 ergibt sich mit obigem Beispiel fur A5

die folgende Struktur:

7

====

===

8

��

��

6 9

KL = reachG(A)(y, {p1, . . . pi}) = reachG/C(S1)(y, C(Si)):

7 2 4 5 1

6 8 3

9

;

7 5

====

===

6 8

9

5.1 Quotient-Tree-Methoden

Idee:

• Eine blockweise Betrachtung der Matrix fuhrt auf den Quotientengraph

• Welche Graphen fuhren offensichtlich zu keinem Fill-In durch die Gauss-Elimination?

* * * ** * * *

* * 0 0 0* * 0 0 0 *

0 0 * * * 00 0 * * 0

* 0 0 * * 0* * 0 0 0 *

6 5 7 1 2 3

8 4

Fuhrt auf den Quotientengraph:

(5, 6, 7) (1, 2) (3, 4)

wwww

wwww

w

(8)

Hier bietet sich eine blockweise Gauss-Elimination an, die den Vorteil hat,daß der zugehorige Graph einfach ist. Der entscheidende Nachteil ist jedoch,daß in einigen Blocken Fill-In entsteht,

”ungeschickt“ ware z.B. der Start mit

der Elimination von (1, 2):

56

(5, 6, 7)

GGGG

GGGG

G(3, 4)

zzzz

zzzz

(8)

Besser ist der Start mit (5, 6, 7), hier entsteht kein Fill-In.

Allgemein stellt sich die Frage, welche Graphen garantiert zu keinem Fill-Infuhren. Auf eine Antwort fuhrt der folgende Graph G(A):

8

��

��

====

===

7

��

��

====

=== 6

4 5 3 2

1

G(A) ist ein Baum, der so nummeriert ist, daß immer Kind < Eltern gilt undin dem alle Ecken keine Nachbarn außer den Eltern haben. Die Wurzel ist dieletzte von der Gauss-Elimination betroffene Ecke, daher entsteht kein Fill-In.Wenn es also gelingt, eine Quotientengraphen zu finden, der ein Baum ist, soentsteht Fill-In nur innerhalb der bereits besetzten Blocke, es entstehen alsokeine neuen Nonzero-Blocke.

Beispiel 5.4

* * * 2 * 0 0 0 0 0* * * 2 2 0 0 0 0 0* * * * 0 0 2 2 2

2 2 * * * 0 0 2 * 2

* 2 * * 0 0 * 2 2

0 0 0 0 0 * * * 2 2

0 0 0 0 0 * * 2 2 *0 0 2 2 * * 2 * 2 2

0 0 2 * 2 2 2 2 * 2

0 0 2 2 2 2 * 2 2 *

In den 2-Feldern entsteht Fill-In

57

Graph: Quotientengraph:

1 2

��

��

3 4

��

��

5 6

��

��

7

8 9 10

;

Gausschritt I (1, 2)

II (3, 4, 5)

SSSSSSSSSSSSSSIII (6, 7)

ooooooooooo

IV (8, 9, 10)

2

Es stellt sich also die Frage, wie man von Anfang an durch Nummerierungeinen Quotient-Tree erreichen kann. Dies ist mit Hilfe von Dissection moglich.

5.2 Dissection

Dissection erlaubt die Nummerierung des folgenden Graphen so, daß einQuotient-Tree entsteht:

��

��

1 2 19

�� 7 8 22

�� 13 14

3 4 20

�� 9 10 23

�� 15 16

5 6 21

�� 11 12 24

�� 17 18

Der Graph wird dazu gruppenweise nummeriert, so daß der jeweils am wei-testen links liegende Knoten, der nicht auf einer Trennlinie die niedrigsteNummer erhalt bis alle Knoten, die nicht auf einer Trennlinie liegen numme-riert sind. Anschlieend werden die Knoten auf den Trennlinien ebenfalls vonlinks nach rechts nummeriert. Es entsteht der folgende Graph:

(1 . . . 6)

MMMMMMMMMM(7 . . . 12) (13 . . . 18)

ppppppppppp

(19 . . . 24)

58

Daraus ergibt sich die folgende Matrix mit optimaler Envelope, in der keinFill-In in nicht-besetzten Blocken entsteht:

6× 6 0 0 2

0 6× 6 0 2

0 0 6× 6 2

2 2 2 6× 6

Durch Nested Dissection, also eine rekursive Wiederholung der Dissection-Methode, lassen sich großere Graphen rekursiv bearbeiten:

3

1

• •5�

�� • •

8�� • • • • •

4_____ _____• •6�

�� •

7______ •

9�� • • • • •

2• •12�

�� •

10

•15�

�� • 17• • • •

11_____ _____• •13�

�� •

14______ •

16�� • • • • •

• • • • • • • • •

Daraus entsteht nach obigem Schema die folgende rekursiv verschachtelteMatrix:

Dissection 17 0 2

0 Dissection 2 2

2 2 2︸︷︷︸Dissection 1

59

Der Quotientengraph hat die folgende Struktur (unvollstandige Darstellung):

1

~~~~

~~~~

BBBB

BBBB

2

��

��

@@@@

@@@@ 17

3

pppppppppppppp

OOOOOOOOOOOOOOO 10

4

��

��

====

=== 7

~~~~

~~~~

BBBB

BBBB

B

5 6 8 9

Der Quotientengraph ist also ein Baum mit zusatzlichen Großelternverbin-dungen, Urgroßelternverbindungen, usw., daher entsteht kein Fill-In.Bei der Gauss-Elimination wird so nummeriert, daß zunachst die großen unddann die kleinen Indizes bearbeitet werden, die Gauss-Elimination beginntalso immer in den Blattern.

5.3 Lokale Pivotstrategien

In der Praxis empfiehlt sich die Anwendung von Pivotstrategien, da das Star-ten der Gauss-Elimination mit einem Index von minimalem Grad meist keinoptimales Ergebnis liefert.

Beispiel 5.5

2

��

��

====

=== 7

��

��

====

===

1

====

=== 4

��

��

5 6

====

=== 8

��

��

3 9

60

Als Matrix:

* * * ** * * ** * * ** * * * * 2

* * *2 * * * * *

* * * ** * * ** * * *

Die 2-Felder bezeichnen den Fill-In, der beim Start mit 5 entsteht.

Ohne Permutation wurde hier kein Fill-In entstehen, das Minimum-Degree-Verfahren fuhrt jedoch zum Start mit 5 als Pivotelement.

2

5.3.1 Hauptkriterien fur lokale Pivot-Strategien

Ein wichtiges Kriterium ist das Markovitz-Kriterium fur allgemeine Matrizen.Wie oben beschrieben, kann der nach dem i. Schritt der Gauss-Eliminationzu bearbeitende Bereich durch eine (n − i) × (n − i)-Matrix beschriebenwerden. Es ensteht also eine Folge von Matrizen A0 = A, A1, . . . An−1. Essoll nun fur jeden Teilschritt eine Permutation PK gesucht werden, so daßP T

KAKPK moglichst wenig Fill-In erzeugt und das Problem numerisch stabilbleibt. Das Gesamtproblem ist damit auf das Teilproblem AK reduziert, PK

muß bestimmt werden.

Definition 5.3 Nonzeros pro Zeile und Spalte

• r(K)i := Anzahl der Nonzeros in Zeile i

• c(K)j := Anzahl der Nonzeros in Spalte j

Im Zuge der Gauss-Elimination mussen c(K)j − 1 Eintrage eliminiert werden,

dazu werden r(K)i − 1 Eintrage zu diesen Zeilen dazuaddiert. Das fuhrt auf

eine vollbesetzte Untermatrix, die minimiert werden muß: (cKj −1)·(r(K)

i −1).

Nach dem Markovitz-Kriterium:a

(K)ij wird als Pivotelement gewahlt, dabei gilt (c

(K)j − 1) · (r(K)

i )!= min.

61

5.3.2 Die algebraische Struktur der Inversen

Gegeben seien eine Matrix A und der zugehorige Adjazenzgraph G(A). Seiaußerdem b ein dunnbesetzter Vektor und J das Pattern S(b).

Definition 5.4 Transitive Hulle, Column Intersection Graph, Closure

• G∗(A) := {(i, j) | pi*→ pj in G(A)}

G∗(A) heißt die transitive Hulle von G(A)

• G+(A) := {(i, j) | i 6= j und pi*→ pj in G(A) uber Ecken < min{i, j}}

• G∩(A) := {(i, j) | i 6= j und ∃k: Aki 6= 0, Akj 6= 0}G∩(A) heißt Column Intersection Graph

• closure(b) :=⋂{I | J ⊆ I und I abgeschlossen}

I heißt abgeschlossen, wenn in G(A) aus I Kanten herausfuhren, aberkeine hineinfuhren

Ausgehend von dieser Definition gelten die folgenden Eigenschaften:

• |AT∗i| · |A∗j| 6= 0 ⇔ (i, j) ∈ G∩(A)

• G∩(A) = G(|AT | · |A|) ≈ G(AT A)

Beispiel 5.6

1 //

��

2 //

��

5

||

4

]]

��

7oo

3!!6aa

J = {5}, closure(5) = {5, 2, 7, 1, 4, 7}.

2

62

Satz 5.3

Wenn Aii 6= 0 fur alle i, dann gilt G(|L|+ |U |) ⊆ G+(A) (ohne Beweis).

2

Satz 5.4

Bei der Gauss-Elimination mit Pivotsuche gilt G(|L| + |U |) ⊆ G+∩ (A) ≈

G+(AT A) (ohne Beweis).

2

Satz 5.5

S(A−1b) ⊆ closure(b). Die Closure entsteht also durch das Hinzunehmen al-ler Indizes zu Wegen, die auf S(b) fuhren.

Beweis:

b =

*0*0*

⇒ closure(b) =

*0*0*0*

closure(b) ist weniger dunn besetzt

Das Gleichungssystem wird nun so permutiert, daß zunachst S(b), dannclosure(b) und dann der Rest kommt:

Ax = b:

(B DC E

)(yz

)=

(d0

)d sind die Indizes zu closure(b), B ∈ Rk×k

C = 0, denn closure(b) = {1 . . . k} und es fuhren keine Kanten von außennach closure(b). Jeder Eintrag cij 6= 0 wurde bedeuten, daß eine Kante voneinem Index /∈ closure(b) zu einem Index ∈ closure(b) existiert, was nichtsein kann. (

B D0 E

)(yz

)=

(d0

)

63

Ez = 0 ⇒ z = 0

⇒ x =

(y0

)y: closure(b)

⇒ S(x) ⊆ closure(b)

2

Satz 5.6

G(A−1) ⊆ G∗(A)

Beweis:Wenn die Kante i → j in G(A−1) enthalten ist, dann hat A−1 in der j. Spalteder i. Zeile einen Nonzero-Eintrag.⇒ (A−1 ej︸︷︷︸

b

)i 6= 0 ⇒ i ∈ S(A−1ej), ej ist Einheitsvektor

⇒ i ∈ closure(ej) nach dem vorhergehenden Satz⇒ es fuhrt ein Weg von i nach j⇒ i → j ∈ G∗(A)

2

Satz 5.7 Vereinfachter Satz

Sei A eine n× n-Matrix, dann existiert das charakteristische Polynom∑ni=0 αiA

i = 0

Mit A−1 multipliziert ergibt sich:αnA

n−1 + αn−1An−2 + . . . + α1I + α0A

−1 = 0 ⇔α0A

−1 = −αnAn−1 − . . .− α2A− α1I ⇔

A−1 = − 1α0

(αnAn−1 + . . . + α1I)

⇒ S(A−1) ⊆⋃n−1

j=0 S(|Aj|)

0. Naherung an S(A−1): Diagonale1. Naherung: S(A)2. Naherung: S(A), S(A2)

2

64

6 Iterative Verfahren zur Losung linearer Glei-

chungssysteme

Der Vorteil dieser Verfahren liegt darin, daß sie nur (Matrix ·Vektor)-Komplexitathaben, die betrachtete Matrix bleibt dabei außerdem konstant. Es stellt sichjedoch die Frage nach der Konvergenz.

6.1 Das GMRES-Verfahren fur allgemeine quadrati-sche Matrizen

Beim GMRES-Verfahren wird ein kleiner Unterraum Um betrachtet und indiesem eine moglichst gute Naherungslosung fur Ax = b gesucht.

Ansatz fur x:

x = Umy =∑m

i=1 yiui, Um ist die Matrix der Basisvektoren

Ziel des Verfahrens:

minx∈Um‖ Ax− b ‖2 = miny‖ A(Umy)− b ‖2

Zur Losung wird hier die Normalgleichung (AUm)T (AUm)y = (AUm)T b her-angezogen. Diese fuhrt auf y und damit auf x. Insgesamt soll also u1 auf x1

fuhren, x1 auf u2, u2 auf x2, usw.Es stellt sich die Frage, welcher Unterraum Um gewahlt werden soll. Es bietetsich eine wachsende Folge von Unterraumen bezogen auf A und b an:

Um = span{b, Ab, A2b, . . . Am−1b} = Km(A, b)

Hierbei entsteht jedoch das Problem, daß die Akb stark linear abhangig wer-den. b, Ab, A2b . . . sind daher eine schlechte Basis.

Erster Schritt: Finden einer Orthonormalbasis fur Um:

u1 := 1‖b‖2

b

u′2 := Au1 − (uT1 Au1) · u1

u2 := 1‖u′2‖2

u′2

u′j := Auj−1 −∑j−1

k=1(uTk Auj−1) · uk

uj := 1‖u′j‖2

u′j

65

Die so gefundene Basis Um erfullt die Gleichung

uT1 u′2 = uT

1 Au1 − (uT1 Au1)(u

T1 u1) = 0

daher ist Um = (u1 . . . um) , span{b, Ab, . . . Am−1b} eine Orthonormalbasis.

Fur das Produkt von A und dem j − 1. Basisvektor ergibt sich:

Auj−1 = u′j︸︷︷︸uj ·‖ u′j ‖2︸︷︷︸

hj,j−1

+∑j−1

k=1 (uTk Auj−1)︸︷︷︸hk,j−1

·uk, j = 1, 2, . . .

=∑k

j=1 hk,j−1 · uk

Fur das Produkt AUm ergibt sich:

A(u1 . . . um) = (Au1 . . . Aum)= (u1 . . . um+1) ·Hm+1,m

= (u1 . . . um+1) ·

h11 . . . . . . h1m

h21...

0...

.... . . hm,m−1 hmm

0 . . . 0 hm+1,m

︸︷︷︸

Hessenbergmatrix

Es ergibt sich also die Gleichung

AUm = Um+1 ·Hm+1,m

Dieses sog. Arnoldi-Verfahren fuhrt zur Orthonormierung von {b, Ab, . . .}:

minx∈Um‖ Ax− b ‖2 = miny‖ A (Umy)︸︷︷︸=x

−b ‖2

= miny‖ (AUm)y − b ‖2

= miny‖ (Um+1 ·Hm+1,m)y − b ‖2

= miny‖ (Um+1 ·Hm+1,m)y − u1‖ b ‖2 ‖2

= miny‖ Um+1(Hm+1,m · y − e1‖ b ‖2) ‖2

= miny‖ Hm+1,m · y − e1‖ b ‖2 ‖2

Zu Losen ist also ein Least-Squares-Problem in einer Hessenbergmatrix, zurLosung wird die QR-Zerlegung Hm+1,m = QR verwendet. Nach m Multipli-kationen mit Givensmatrizen G12, G23, . . .Gm−1,m, Gm,m+1 ergibt sich:

66

(Gm,m+1 · . . . ·G23 ·G12)︸︷︷︸Q

·Hm+1,m =

* *

. . .

*0 . . . 0

=

(R0

)

Hm+1,m hat also die QR-Zerlegung

Hm+1,m = QT

(R0

)Ausgehend vom oben formulierten Ziel des Verfahrens ergeben sich die fol-genden Gleichungen:

minx∈Um ‖ Ax− b‖2 = miny ‖ Hm+1,m · y − e1‖ b ‖2‖2

= miny ‖ QT

(R0

)y − e1‖ b ‖2‖2

= miny ‖(

Ry0

)− (Qe1)‖ b ‖2‖2

= miny ‖ Ry−‖b‖2(Qe1)′

−‖b‖2(Qe1)m‖2

Die Losung der Gleichung Ry = ‖ b ‖2(Qe1)′ fuhrt also auf y und damit auf

x = Umy.

GMRES ist eine geschickte Implementierung der folgenden Schritte:

1. Berechnen von Hm+1,m und der Orthonormalbasis u1 . . . um

2. QR-Zerlegung von Hm+1,m in Q und R

3. Auflosen des Dreieckssystems Ry, das eine Naherungslosung fur Um

ergibt

6.1.1 GMRES Restarted

GMRES Restarted bzw. GMRES(m) ist eine wiederholte Anwendung desGMRES-Verfahrens, dabei wird als Startwert das Ergebnis xm der letztenNaherung verwendet. Dieses Verfahren ist nur bis m = 20 effizient, fur große-re m wird es zu teuer.

67

6.1.2 Fehlerabschatzung

Zur Abschatzung des beim GMRES-Verfahren entstehenden Fehlers wirduber einem Krylov-Raum minimiert. Daraus ergeben sich die folgenden Glei-chungen:

‖ rm‖2 = ‖ Axm − b ‖2, xm ist Naherung nach m Schritten GMRES= minx∈Um ‖ Ax− b‖2

= minαj‖ A(

∑m−1j=0 αj AjB︸︷︷︸

Basis von Um

)− b‖2

= minαj‖ (A ·

∑m−1j=0 αjA

j) · b− b‖2

= minpm−1 ‖ Apm−1(A)b− b‖2

= minqm(0)=1 ‖ qm(A)b‖2, b =∑m

j=1 βjuj

Eigenwertzerlegung A = UDUH :Die Matrix U = (u1 . . . um) besteht aus den Eigenvektoren ui von A. DieDiagonalmatrix

D =

d1

. . .

dm

enthalt auf der Diagonalen die Eigenwerte di von A. Es gilt daher

Aui = diui, b =∑m

j=1 βjuj

Ausgehend von der obigen Minimierungsgleichung ergibt sich:

‖ rm‖2 = minqm(0)=1 ‖ qm(A)b‖2

= minqm(0)=1 ‖ qm(A) ·∑m

j=1 βjuj‖2

= minqm(0)=1 ‖∑m

j=1 βjqm(A)uj‖2

= minqm(0)=1 ‖∑m

j=1 βjqm(dj)uj‖2

≤ Ax− b≤ minqm(0)=1 maxm

j=1 ‖ qm(dj)‖2 · ‖ b ‖2︸︷︷︸‖r0‖2

Beispiel 6.1

A habe nur die zwei Eigenwerte λ1 und λ2. Durch die Wahl von

q2(x) = (x−λ1)(x−λ2)(−λ1)(−λ2)

68

ergibt sich fur die geforderten Bedingungen:

q2(0) = −λ1

−λ1· −λ2

−λ2= 1, q2(λ1) = q2(λ2) = 0

Dieses Polynom ist also fur m = 2 zulassig, GMRES konvergiert nach zweiSchritten. Falls die Eigenwerte von A nahe beieinander und getrennt von 0(qm(0) = 1) sind, konvergiert GMRES schnell.

Ungunstig ist dagegen folgender Fall:0 11 000 0 1 0

Hier benotigt GMRES m Schritte, um zu konvergieren.

2

6.2 Prakonditionierung

Eine iterative Losung ist nur dann sinnvoll, wenn das angewandte Verfahren

”schnell“ konvergiert. Die Konvergenz ist dabei von der Lage der Eigenwerte

abhangig. Fur regulares M gilt die Gleichung

Ax = b ⇔ MAx = MbX

durch die die Lage der Eigenwerte verandert wird, MA hat also andere Ei-genwerte als A.Das Ziel ist also, M und Mb schnell zu berechnen, außerdem soll MA gunstigeEigenwerte haben. Die Verbesserung der Kondition (also der Lage der Eigen-werte) fuhrt somit auf die Notwendigkeit einer Prakonditionierung.

6.2.1 Wesentliche Prakonditionierer

1. Jacobi-PrakonditionierungAusgehend von der Diagonalmatrix D = diag(A) wird hier M = D−1

verwendet

Ax = b ⇔ D−1Ax = D−1b ⇔(D− 1

2 AD− 12 )D

12 x = D− 1

2 b ⇔(AD−1) (Dx)︸︷︷︸

y

= b

69

2. Gauss-Seidel-PrakonditionierungHier wird nach der abgewandelten LU -Zerlegung A = D + L + U alsFaktor M = (D + L)−1 verwendet

(D + L)−1 · Ax = (D + L)−1 · b

3. Unvollstandige Berechnung der LU-Gauss-EliminationAm einfachsten ist die reduzierte Berechnung der LU -Zerlegung, z.B.auf das Pattern von A. Diese Methode wird auch als Incomplete LU-Decomposition (ILU) bezeichnet.

Ausgehend von dem Ansatz

MAx = Mb

ware M = A−1 ideal. Leichter zu erreichen ist jedoch M ≈ A−1, gesucht istalso min ‖ AM − I‖F . Die Frobenius-Norm bietet sich hier an, da sie leichtberechen- und differenzierbar ist.Betrachtet werden also Kandidaten zur Minimierung mit vorgegebenem Pat-tern z.B. G(A) = G(M), gesucht wird minM,G(M)=G(A) ‖ AM − I‖F . Nachder Anwendung vom Jacobi, Gauss-Seidel bzw. ILU gilt M ≈ A, es muß alsonoch invertiert werden.

‖(

a11 a12

a21 a22

)‖2

F = (a211 + a2

12) + (a221 + a2

22)

= ‖(

a11

a21

)‖2

2+ ‖(

a12

a22

)‖2

2

Fur die Minimierung ergibt sich:

‖ AM − I‖2F = ‖ (AM1, . . . AMn)− (e1, . . . en)‖2

F

= ‖ (AM1 − e1, AM2 − e2, . . . AMn − en)‖2F

= ‖ AM1 − e1‖22 + . . . + ‖ AMn − en‖2

2

=∑n

i=1 ‖ AMi − ei‖22

Fur M ergibt sich also die folgende spaltenweise Darstellung:

M = (M1, . . . Mn)

M wird somit bestimmt durch n unabhangige Optimierungsprobleme derForm min ‖ AMi − ei‖2. Mi hat dabei die Besetztheitsstruktur der i. Spaltevon A.

70

Beispiel 6.2

‖

2 −1−1 2 −1

−1 2. . .

−1−1 2

·

000

mi−1

mi

mi+1

0...0

−

0001000

‖2 =

= ‖

−12 −1−1 2 −1

−1 2−1

·

mi−1

mi

mi+1

−

010

‖2

Die Losung des Teilproblems ‖ A′m′− b′‖2 liefert mi−1, mi und mi+1. Darauskann dann Mi und schließlich M berechnet werden.

2

6.3 Zusammenfassung

Fur allgemeines A mit A = AT > 0 existieren verschiedene iterative Verfah-ren, das wichtigste ist das oben beschriebene GMRES bzw. GMRES Restar-ted. Desweiteren gibt es die Verfahren BCG, BiCGSTAB und QMR. DieseVerfahren sind zwar nicht optimal bzgl. einer Norm, dafur in jedem Schrittbilliger.

71

A Literatur

• George Liu: Computer Solution Of Large Sparse Positive Definite Sy-stems

• Duff, Erisman, Reid: Direct Methods For Sparse Matrices

• Yousef Saad: Iterative Methods For Sparse Linear Systems

• Barrett, Berry, Chan, Demmel et al.: Templates For The Solution OfLinear Systems

• Hackbusch: Iterative Losung großer linearer Gleichungssysteme

72

Date post:	12-May-2019
Category:	Documents
Upload:	dangcong
View:	222 times
Download:	0 times

Dunnbesetzte Matrizen: Algorithmen und¨ Datenstrukturen · 1 Motivation, Beispiele und...

Documents