2. Direkte Verfahren zur L¨ osung
linearer Gleichungssysteme
Einleitung (1)
Eine zentrale Rolle bei numerischen Berechnungen spielen lineare Glei- chungssysteme
• Es sind die am h¨aufigsten auftretenden numerischen Probleme
• Anwendungsgebiete sind z.B.
* fast alle naturwissenschaftlich-technischen Problemstellungen vom Wetterbericht bis zur W¨armeentwicklung auf einer Koch- platte oder der Planung der Leiterbahnen auf Mikrochips,
* Bildverarbeitung oder z.B Beleuchtungsprobleme in der Com- putergrafik,
* wirtschaftlichen Fragestellungen wie Versicherungskosten oder B¨orsenkursvorhersage.
Einleitung (2)
L¨osungsverfahren
• Die direkten Verfahren liefern eine mit Rundungsfehlern behaftete L¨osung nach endlich vielen Schritten.
• Die iterativen Verfahren beginnen mit einer Anfangsn¨aherung und produzieren eine verbesserte N¨aherungsl¨osung nach endlich vielen Schritten.
• Falls m¨oglich wird das Problem mit einem direkten Verfahren be- rechnet und anschließend werden die Rundungsfehler mit einem iterativen Verfahren verringert.
Einleitung (3)
Problemstellung: Berechne den Vektor x = (x1, x2, . . . xn) aus
a1,1x1 + a1,2x2 + · · ·a1,nxn = b1 a2,1x1 + a2,2x2 + · · ·a2,nxn = b2
·
·
an,1x1 + an,2x2 + · · ·an,nxn = bn
oder in Matrix-Schreibweise
Ax = b
Bemerkung: Vektoren werden hier ohne Vektorpfeil geschrieben
Einleitung (4)
• Es existieren entweder keine, eine oder unendlich viele L¨osungen.
• Das Gleichungssystem hat eine L¨osung, wenn die Inverse Matrix A−1 zu A existiert.
• Die L¨osung kann mit Hilfe von A−1: A−1Ax = A−1b = x.
• Dieses ist der Fall, wenn die Matrix nicht singul¨ar ist, d.h. die Determinante ungleich Null ist (Leibniz-Formel)
detA = X
π∈Sn
(signπ)a1,π(1)a2,π(2). . . an,π(n).
* π(1), . . . π(n) bedeutet eine Permutation der Zahlen 1 bis n.
* signπ ist der Vorzeichen der Permutation
* Pπ∈Sn bedeutet die Summe ¨uber alle Permutationen
• Es wird im folgenden vorausgesetzt, dass die Matrizen nicht sin-
Gauß-Verfahren (1)
Das Gauß-Eliminationsverfahren ist bekannt aus der Mathema- tikvorlesung
Dieses direkte L¨osungsverfahren bringt das Gleichungssystem in Drei- ecksform und berechnet den L¨osungsvektor.
Schulbeispiel
E1 : x1 + x2 + 3x4 = 4 E2 : − x2 − x3 − 5x4 = −7 E3 : − 4x2 − x3 − 7x4 = −15 E4 : 3x2 + 3x3 + 2x4 = 8
Gauß-Verfahren (2)
Matrix A:
A =
1 1 0 3
0 −1 −1 −5 0 −4 −1 −7
0 3 3 2
,
Determinante:
detA = a1,1 ∗ a2,2 ∗ a3,3 ∗ a4,4 − a1,1 ∗ a2,2 ∗ a3,4 ∗ a4,3
− a1,1 ∗ a2,3 ∗ a3,2 ∗ a4,4 + a1,1 ∗ a2,3 ∗ a3,4 ∗ a4,2 + a1,1 ∗ a2,4 ∗ a3,2 ∗ a4,3 − a1,1 ∗ a2,4 ∗ a3,3 ∗ a4,2
= 2 − 21 − 8 + 21 + 60 − 15 = 39
Gauß-Verfahren (3)
Erlaubte Transformationen zur L¨osung des Gleichungssystems
• Multiplizieren einer Zeile (Gleichung) mit einer Zahl verschieden von Null
• Addieren eines Vielfachen einer Zeile zu einer anderen Zeile
• Vertauschen von Zeilen (Gleichungen) bzw. Spalten (Unbekann- ten, entspricht Umnummerierung)
Mit Hilfe der Transformationen reduziere Gleichungssystem auf ein Dreieckssystem.
Gauß-Verfahren (4)
• Durchlaufe alle Zeilen.
• Bei jeder Zeile j unterhalb der aktuellen Zeile i ersetzt die Ele- menten durch
Zeile j → Zeile j - Zeile i ·(aj,i/ai,i)
• Ersetze bj ebenfalls durch bj → bj − bi · aj,i/ai,i
• Dieser Schritt ¨andert den L¨osungsvektor nicht!
• Es f¨uhrt dazu, dass alle Elemente der i-ten Spalte unterhalb der i-ten Zeile zu Null werden.
aj,i = aj,i − ai,i · (aj,i/ai,i) = 0 f¨ur j > i
Gauß-Verfahren (5)
Reihe i
Spalte i Reihe j
0 a
j,i
Gauß-Verfahren (6)
Zweite Spalte des Schulbeispiels:
1 1 0 3
0 −1 −1 −5 0 −4 −1 −7
0 3 3 2
·
x1 x2 x3 x4
=
4
−7
−15 8
·
a3,i = a3,i − a3,2/a2,2 · a2,i b3 = b3 − a3,2/a2,2 · b2 a4,i = a4,i − a4,2/a2,2 · a2,i b4 = b4 − a4,2/a2,2 · b2
Gauß-Verfahren (7)
Allgemein:
1. Die Matrix A ¨andert sich mit jedem Eliminationsschritt 2. Starte mit der Ursprungsmatrix A(1) = A
3. F¨uhre die nachfolgenden Schritte f¨ur alle Zeilen k = 1, . . . , n−1 zur Berechnung einer ver¨anderten Matrix A(k+1) aus A(k) durch, die in der Spalte unter dem Diagonalelement ak,k nur Nullen stehen hat
4. Berechnet die Zahlen li,k aus dem Quotienten der Elemente a(k)i,k und a(k)k,k nach
li,k = a(i,kk) a(k)k,k
, i = k + 1, . . . , n
Gauß-Verfahren (8)
5. Bestimme die ver¨anderte Matrizen A(k+1) = a(i,jk+1), die sich aus dem k-ten Schritt des Gauß-Verfahrens ergibt durch
a(k+1)i,j =
a(k)i,j − li,ka(k)k,j , i = k + 1, . . . , n; j = k, . . . , n a(i,jk) sonst
6. Ver¨andere die Vektor b(k)i zu b(k+1)i , so dass der Ergebnisvektor unver¨andert bleibt, gem¨aß
b(ik+1) =
b(k)i − li,kb(k)k , i = k + 1, . . . , n b(k)i sonst
Gauß-Verfahren (9)
C-Code
// for each row "k" except the last for (k=0; k<n-1; k++) {
// step through subsequent rows "i"
for (i=k+1; i<n; i++)
l[i][k] = a[i][k]/a[k][k];
// step through subsequent rows "i"
for (i=k+1; i<n; i++) {
// step through subsequents cols "j" of "i"
for (j=k; j<n; j++)
// modified n-k-1 elements of row "i"
a[i][j] = a[i][j] - l[i][k]*a[k][j];
// modify right side
b[i] = b[i] - l[i][k]*b[k];}}
F¨ur j = k gilt
a[i][k] = a[i][k] - (a[i][k]/a[k][k]) * a[k][k] = 0;
d.h. die Spalte a[i][k] mit i>k wird auf 0 gesetzt.
Gauß-Verfahren (10)
Schulbeispiel wird zu
E1 : x1 + x2 + 3x4 = 4
E2 : − x2 − x3 − 5x4 = −7
E3 : 3x3 + 13x4 = 13
E4 : − 13x4 = −13
oder allgemein
a1,1 a1,2 · · · a1,n a2,2 · · · a2,n
. . . ...
an,n
·
x1 x2 ...
xn
=
b1 b2 ...
bn
·
Gauß-Verfahren (11)
Ein Dreieckssystem ist leicht zu l¨osen. Aus E4: x4 = 1
x4 einsetzten in E3:
3x3 + 13 = 13 → x3 = 0 x3, x4 einsetzten in E2:
−x2 − 5 = −7 → x2 = 2 x2, x3, x4 einsetzten in E1:
x1 + 2 + 3 = 4 → x1 = −1
Gauß-Verfahren (12)
Allgemein: Ausgehend von einer Dreiecksmatrix A gilt xn = 1
an,n · bn xn−1 = 1
an−1,n−1 · bn−1 − an−1,nxn, · · · bzw.
xk = 1 ak,k ·
bk −
n X j=k+1
ak,jxj
f¨ur k = n, . . . , 1.
Wichtig: Die Diagonalelemente akk, die Pivotelemente m¨ussen 6= 0 sein bzw.
det(A) = a11 · a22 · . . . · ann 6= 0.
Gauß-Verfahren (13)
Ordnung des Gauß-Algorithmus
• Eliminationsschritt: Es m¨ussen f¨ur k = n,· · · ,1 neue Werte f¨ur k2 Matrixeintr¨age berechnet werden, insgesamt
n X k=1
k2 = n(n + 1)(2n + 1)/6 ≈ n3/3
• L¨osung des Dreiecks-Gleichungssystem: Es m¨ussen zur Berech- nung der k-ten Unbekannten k Terme zusammengefasst werden, insgesamt
n X k=1
k = n(n + 1)/2 ≈ n2/2
• Ordnung insgesamt O(n3).
Pivotsuche (1)
Problem:
1. Die “Pivotelemente” a(k,kk) k¨onnen in jedem Schritt gleich Null sein.
2. Ist das Pivotelemente a(k)k,k viel kleiner als a(k)i,k , wird li,k = a
(k) i,k
a(k)k,k
sehr groß und Rundungsfehler zerst¨oren die L¨osung Beispiel:
• Berechne mit 4-stelliger Gleitpunktarithmetik
−10−5 1
·
x1
=
1
·
Pivotsuche (2)
• L¨osung mit Gauß-Algorithmus:
−10−5 1 1
2 1 0
→
−10−5 1 1
0 200000 200000
mit der L¨osung x = (0,1)
• Berechne das gleiche System mit vorheriger Vertauschung der Zei- len 1 und 2:
2 1 0
−10−5 1 1
→
2 1 0 0 1 1
mit der L¨osung x = (−0.5,1)
• Exakte L¨osung: x = (−0.4999975. . . ,0.999995 . . .)
Pivotsuche (3)
Ausweg aus beiden Problemen:
1. Spaltenpivotsuche:
Bestimme das betragsm¨aßig gr¨oßte Element a(k)r,k , k ≤ r ≤ n und vertausche Zeile r mit Zeile k und b(rk) mit b(kk), falls r 6= k.
2. Zeilenpivotsuche:
Bestimme das betragsm¨aßig gr¨oßte Element a(k,rk), k ≤ r ≤ n und vertausche Spalte r mit Spalte k, falls r 6= k. Dies entspricht einer Umnummerierung des L¨osungsvektors.
3. Totalpivotsuche:
Das ist die Kombination aus Zeilen- und Spaltenpivotsuche.
Pivotsuche (4)
Algorithmus f¨ur Spaltenpivotsuche:
for k = 1 to n − 1 do
bestimme s mit |as,k| = max{|ai,k|, i = k, . . . , n}
vertausche Zeilen s und k von A.
· · ·
Praxis: Die Zeilen werden nicht vertauscht, sondern ¨uber einen In- dexvektor angesprochen oder die Zeiger auf die Zeilen “umgeh¨angt”, um Speicherzugriffe und eventuell Kommunikation zu sparen.
Pivotsuche (5)
Algorithmus mit Indexvektor:
Indexvektor ind(i) = i, i = 1, . . . , n for k = 1 to n − 1 do
bestimme s mit |aind(s),k| = max{|aind(i),k|, i = k, . . . , n}
vertausche ind(s) mit ind(k)
· · ·
for i = k + 1 to n do for j = k to n do
aind(i),j = aind(i),j − lind(i),kaind(k),j
· · ·
In C einfach a[k] ↔ a[s]
Es gibt viele L¨osungen im Netz, z.B. unter
Kondition einer Matrix (1)
Der Einfachheit halber sei A exakt gegeben Ungenaue Eingabe in B: ∆b
Ungenaue Ausgabe in x: ∆x Zu berechnen: ||∆x||/||x||
Aus Ax = b folgt A(x + ∆x) = b + ∆b oder A∆x = ∆b und aus
||∆x|| = ||A−1∆b|| ≤ ||A−1|| · ||∆b||
||b|| = ||Ax|| ≤ ||A||||x||
folgt
||∆x||
||x|| ≤ ||A−1|| · ||A|| · ||∆b||
||b||
κ(A) = ||A−1|| · ||A|| muss klein sein. Was ist unter ||A|| zu verstehen?
Kondition einer Matrix (2)
Vektornorm: Unterschiedliche Definitionen von Normen sind m¨oglich.
• 1-Norm:
||x||1 =
n X i=1
|xi|
• 2-Norm oder euklidische Norm:
||x||2 =
v u u t
n X i=1
x2i
• ∞-Norm:
||x||∞ = max
i=1,...n |xi|
Analog f¨ur Matrizen, ¨uber Spaltensumme, Spektralradius bzw. Zeilen- summe. Diese Gr¨oßen sollten vor der L¨osung des Systems berechnet
Kondition einer Matrix (3)
Analog ohne Beweise:
• 1-Norm, Spaltennorm:
||A||1 = max
j=1,...n n X i=1
|ai,j|
• 2-Norm, Spektralnorm:
||A||2 = maxx6=0||Ax||
||x|| =
q
ρ(ATA)
• ∞-Norm, Zeilensummennorm:
||A||∞ = max
i=1,...n n X j=1
||ai,j||
Gauß-Elimination mit Pivotsuche gilt als numerisch stabil!
LU (LR)-Zerlegung (1)
Die Matrix A l¨asst sich als Produkt zweier Dreiecksmatrizen schreiben:
A = LU
Das Verfahren wird als LU-Zerlegung bezeichnet und entspricht der Gauß-Elimination. Es gilt
L =
1 0 0 · · · 0 l2,1 1 0 · · · 0 l3,1 l3,2 1 · · · 0 ... ... ... ... ...
ln,1 ln,2 ln,3 · · · 1
U =
a(1,1n) a(1,2n) a(1,3n) · · · a(1,nn) 0 a(n)2,2 a(n)2,3 · · · a(n)2,n 0 0 a(n)2,3 · · · a(n)2,n ... ... . . . ... ...
0 0 0 · · · a(n)n,n
LU-Zerlegung (2)
Ohne Beweis: Die Matrixeintr¨age li,j von L sind die Gewichte der Gauß-Elimination und Einsen auf der Diagonalen.
Unser Schulbeispiel:
1 0 0 0 0 1 0 0 0 4 1 0 0 −3 0 1
·
1 1 0 3
0 −1 −1 −5
0 0 3 13
0 0 0 −13
=
1 1 0 3
0 −1 −1 −5 0 −4 −1 −7
0 3 3 2
.
LU-Zerlegung (3)
Vorteil:
Ist A einmal zerlegt, l¨asst sich das System Ax = b f¨ur alle b schnell l¨osen, ohne das noch einmal eine Gauß-Elimination durchgef¨uhrt wer- den muss.
Ax = LUx = b
L¨ose zuerst Ly = b und anschließend Ux = y. Beides sind Probleme mit Dreiecksmatrizen und damit von der Ordnung
O(n2).
HPL
Der HPL-Benchmark der TOP500 Liste f¨ur die weltweit leis- tungsf¨ahigsten Rechner ist L¨osung eines linearen Gleichungs- systems ¨uber eine LU-Zerlegung mit Spaltenpivotsuche.
1. Das Programm ist mit MPI (siehe Masterveranstaltung Parallel Computing) geschrieben.
2. Es fasst zur Beschleunigung, d.h. zur Optimierung des Speicher- zugriffs mehrere Spalten in Bl¨ocken zusammen.
3. Es verwendet BLAS Bibliotheksfunktionen (Basic Linear Algebra Subprograms).
Cholesky-Zerlegung
F¨ur sogenannte symmetrische und positiv definite Matrizen, also Ma- trizen, f¨ur die gilt
AT = A und xTAx > 0
gibt es eine rechts-obere Dreiecksmatrix R mit rii < 0 und A = RTR
Der zugeh¨orige Algorithmus ist ca. doppelt so schnell wie die Gauß- Elimination.
Solche Matrizen kommen in der Praxis sehr h¨aufig vor, z.B. bei im Bereich der Computergraphik.
QR-Zerlegung (1)
Bei der QR-Zerlegung A = QR handelt es sich um eine Zerlegung in 2 Matrizen, von denen eine eine sogenannte orthogonale Matrix ist, also Q−1 = QT und die andere eine rechts-obere Dreiecksmatrix R Vorteil:
• Die Gauß-Elimination bzw. LR-Zerlegung kann die Kondition der Matrix stark ¨andern, so dass es trotz Pivot-Suche zu hohen Run- dungsfehlern kommen kann.
• Die QR-Zerlegung kann auch auf nicht-quadratische Matrizen an- gewandt werden, so wie sie bei Ausgleichsrechnungen vorkommen.
• Es lassen sich ¨uber die QR-Zerlegung die Eigenwerte der Matrix berechnen.
QR-Zerlegung (1)
Nachteil:
• Das Verfahren ist aufw¨andiger als die LR-Zerlegung
Ist das System einmal zerlegt
Ax = QRx = b,
dann l¨ose zuerst Qy = b ¨uber y = Q−1b = QTb und anschließend Rx = y. Beides sind Probleme von der Ordnung
O(n2).
Beispiel einer großen Matrix: Das Radiosity-Verfahren (1)
• Globales Beleuchtungsmodell
• Berechnung der diffus abgestrahlten Energie zwischen Objekten einer Szene
• Energieerhaltung: Summen der abgestrahlten und aufgenommenen Energien sind gleich
• Lichtquellen sind auch Objekte der Szene
Das entstehende Gleichungssystem wird aus Zeitgr¨unden meist iterativ und nicht mit dem Gauß-Verfahren gel¨ost.
siehe http://de.wikipedia.org/wiki/Radiosity_(Computergrafik)
Beispiel: Radiosity-Verfahren (2)
• N Objekte in einer Szene Gr¨oße Bedeutung
Ai Oberfl¨ache von Objekt i
Bi Von i abgestrahlte Leistung/Oberfl¨ache Fi,j Formfaktor: Anteil, der von Bi auf Aj trifft Ei von i erzeugte Strahlungsleistung/Fl¨ache Ri Reflektionskoeffizient von i
• Gesamte Strahlungsleistung von Objekt j, die das Objekt i trifft ist damit BjAjFj,i
• Gesamte abgestrahlte Energie pro Oberfl¨ache ist die Summe aus der erzeugten und der reflektierten Energie.
Beispiel: Radiosity-Verfahren (3)
• Energiebilanz:
Bi = Ei +
n X j=1,j6=i
Ri(BjAjFj,i)/Ai
• Der Lichtweg ist umkehrbar → Fj,iAj = Fi,jAi
• Es findet keine Absorption statt → Fi,i = 0
• Energiebilanz:
Bi = Ei +
n X j=1
Ri(BjAjFj,i)/Ai
• Gleichungssystem:
Ei =
n X j=1
(δi,j − RiFi,j)Bj oder Matrix: Ai,j = δi,j − RiFi,j
Erg¨ anzung: Eigenwerte (1)
Ein Gleichungssystem hat eine eindeutige L¨osung, wenn die Matrix A nicht singul¨ar ist, d.h. A−1 existiert. Dies ist der Fall, wenn
• die Determinante ungleich Null oder
• alle Eigenwerte der Matrix ungleich Null sind.
Eine n×n Matrix A f¨uhrt bei Multiplikation mit einem Vektor x diesen in einen Vektor y ¨uber. Sind x und y parallel ist x ein Eigenvektor von A zum Eigenwert λ
Ax = y = λx
Eigenwerte (2)
Die Eigenwerte k¨onnen bestimmt werden ¨uber das charakteristische Polynom, da (A − λI)x = 0
P(λ) = det(A − λI) = det
a1,1 − λ a1,2 · · · a1,n a2,1 a2,2 − λ · · · a2,n
... . . . ...
an,1 an,n − λ
= 0,
wobei I die Einheitsmatrix ist.
• x ist nur bis auf eine Konstante definiert
• Die Eigenwerte k¨onnen komplex sein
• Die Matrix (A − λI) ist singul¨ar
Eigenwerte (3)
Eigenwerte und Eigenvektoren haben vielf¨altige Anwendungen in z.B.
• Physik (Schwingungen, Drehungen, Quantenmechanik)
• Maschinenbau (Festigkeitslehre und Knicklasten),
• Biologie und Wirtschaftswissenschaften (Entwicklung eines biolo- gischen bzw. wirtschaftlichen Systems ¨uber Wahrscheinlichkeits- matrizen bzw. Markov-Ketten),
• Bildbearbeitung, z.B. Objektausrichtung
• PageRank einer Homepage als Eigenvektor der Google-Matrix
• Regelungstechnik und vielem mehr,
Eigenwerte (4)
Drei einfache L¨osungsans¨atze:
• Bei kleinen Systemen kann die Eigenwertgleichung direkt berech- net werden.
• Große Systeme sind meist symmetrisch, d.h. A = AT. Dann exis- tieren verschiedene schnelle Verfahren zur Berechnung der Eigen- werte ¨uber die QR-Zerlegung.
• Wird nur eine wesentliche Eigenschaft ben¨otigt, z.B. der gr¨oßte Eigenwert, kann dieser leicht durch die Potenzmethode bestimmt werden.
Eigenwerte (5)
Potenzmethode:
• Voraussetzung (normalerweise g¨ultig): Alle Eigenvektoren sind li- near unabh¨angig.
• Dann l¨asst sich ein beliebiger Vektor x schreiben als x =
n X j=1
ajv(j)
• Annahme: Die Eigenvektoren v(j) geh¨oren zu den der Gr¨oße nach sortierten Eigenwerten |λ1| > |λ2|. . . ≥ |λj| ≥ . . . ≥ |λn| ≥ 0 und der gr¨oßte Eigenwert kommt nur einmal vor.
• Dann gilt
Ax =
n
X ajAv(j) =
n
X ajλjv(j)
Eigenwerte (6)
Dieser Ausdruck wird nun immer wieder mit A multipliziert:
A2x =
n X j=1
ajA2v(j) =
n X j=1
ajλ2jv(j)
Akx =
n X j=1
ajAkv(j) =
n X j=1
ajλkjv(j)
= λk1
n X j=1
ajλkj
λk1v(j)
k→∞lim Akx = lim
k→∞λk1a1v(1)
Daraus kann der gr¨oßte Eigenwert mit dem zugeh¨origen Eigenvektor gewonnen werden.
