Gauss Elimination mit Pivotisierung¶
Lernziel:
Algorithmus der Gauss Elimination schrittweise zu verstehen
Kennenlernen der Idee wie das Kopieren der Zeilen vermieden werden kann
Wir betrachten wieder dasselbe lineare Gleichungssystem:
[1]:
import numpy as np
[2]:
n = 4
A = np.array([[2, -1, -3, 3],
[4, 0, -3, 1],
[6, 1, -1, 6],
[-2, -5, 4, 1]],dtype=float)
R = np.array(A) # R-Matrix
L = np.zeros((n,n)) # L-Matrix
A
[2]:
array([[ 2., -1., -3., 3.],
[ 4., 0., -3., 1.],
[ 6., 1., -1., 6.],
[-2., -5., 4., 1.]])
In der Gauss-Elimination soll das Tauschen der Matrixeinträge vermieden werden. Dazu benutzt man die ind-Liste, in der die Reihenfolge der Zeilen gespeichert wird. Initital ist diese gegeben durch:
[3]:
ind = np.arange(n)
ind
[3]:
array([0, 1, 2, 3])
LR - Faktorisierung¶
Erste Spalte¶
Pivotelement ist die grösste Zahl einer Spalte. Wir suchen daher die Zeile mit dem grössten Betrag, wobei als Reihenfolge der Zeilen der Matrix der Vektor ind benutzt wird.
Der Index j geht über die Spalten der Matrix, daher
[4]:
j = 0
Gesucht ist daher die Zeile mit dem grössten absoluten Wert der ersten Spalte:
[5]:
indPivot = j
for i in range(j,n):
if np.abs(R[ind[i],j]) >= np.abs(R[ind[indPivot],j]):
indPivot = i
indPivot
[5]:
2
das geht auch kurz und viel schneller mit:
[6]:
j+np.argmax(abs(R[ind[j:],j]))
[6]:
2
Die Reihenfolge wird nun aktualisiert:
[7]:
i = ind[indPivot] # Wert sichern
ind[indPivot] = ind[j] # neue Reihenfolge
ind[j] = i
ind
[7]:
array([2, 1, 0, 3])
Über die restlichen Zeilen (bezogen auf die Einträge im Vektor ind) wird die Gauss Elimination (Algorithmus 2.2) angewandt:
[8]:
for i in range(j+1,n):
lij = R[ind[i],j]/R[ind[j],j]
L[ind[i],j] = lij # speichern der unteren Dreiecksmatrix
R[ind[i],j:] -= lij*R[ind[j],j:] # bearbeitet die Zeile ind(i)
R
[8]:
array([[ 0. , -1.33333333, -2.66666667, 1. ],
[ 0. , -0.66666667, -2.33333333, -3. ],
[ 6. , 1. , -1. , 6. ],
[ 0. , -4.66666667, 3.66666667, 3. ]])
In der neuen Ordnung:
[9]:
R[ind,:]
[9]:
array([[ 6. , 1. , -1. , 6. ],
[ 0. , -0.66666667, -2.33333333, -3. ],
[ 0. , -1.33333333, -2.66666667, 1. ],
[ 0. , -4.66666667, 3.66666667, 3. ]])
Zweite Spalte¶
Procedere für die zweite Spalte
[10]:
j = 1
Pivotelement bestimmen:
[11]:
ind
[11]:
array([2, 1, 0, 3])
[12]:
indPivot = j+np.argmax(abs(R[ind[j:],j]))
indPivot
[12]:
3
Die Reihenfolge wird nun aktualisiert (kompakt geschrieben):
[13]:
ind[indPivot], ind[j] = ind[j], ind[indPivot] # neue Reihenfolge
ind
[13]:
array([2, 3, 0, 1])
Gauss Elimination (Algorithmus 2.2):
[14]:
for i in range(j+1,n):
lij = R[ind[i],j]/R[ind[j],j]
L[ind[i],j] = lij # speichern der unteren Dreiecksmatrix
R[ind[i],j:] -= lij*R[ind[j],j:] # bearbeitet die Zeile ind(i)
R
[14]:
array([[ 0. , 0. , -3.71428571, 0.14285714],
[ 0. , 0. , -2.85714286, -3.42857143],
[ 6. , 1. , -1. , 6. ],
[ 0. , -4.66666667, 3.66666667, 3. ]])
In der neuen Ordnung:
[15]:
R[ind,:]
[15]:
array([[ 6. , 1. , -1. , 6. ],
[ 0. , -4.66666667, 3.66666667, 3. ],
[ 0. , 0. , -3.71428571, 0.14285714],
[ 0. , 0. , -2.85714286, -3.42857143]])
Dritte Spalte¶
Procedere für die dritte Spalte
[16]:
j = 2
[17]:
indPivot = j+np.argmax(abs(R[ind[j:],j]))
indPivot
[17]:
2
Die Reihenfolge wird nun aktualisiert:
[18]:
ind
[18]:
array([2, 3, 0, 1])
[19]:
ind[j]
[19]:
0
[20]:
ind[indPivot], ind[j] = ind[j], ind[indPivot] # neue Reihenfolge
ind
[20]:
array([2, 3, 0, 1])
Gauss Elimination (Algorithmus 2.2):
[21]:
for i in range(j+1,n):
lij = R[ind[i],j]/R[ind[j],j]
L[ind[i],j] = lij # speichern der unteren Dreiecksmatrix
R[ind[i],j:] -= lij*R[ind[j],j:] # bearbeitet die Zeile ind(i)
R
[21]:
array([[ 0. , 0. , -3.71428571, 0.14285714],
[ 0. , 0. , 0. , -3.53846154],
[ 6. , 1. , -1. , 6. ],
[ 0. , -4.66666667, 3.66666667, 3. ]])
In der neuen Ordnung:
[22]:
R[ind,:]
[22]:
array([[ 6. , 1. , -1. , 6. ],
[ 0. , -4.66666667, 3.66666667, 3. ],
[ 0. , 0. , -3.71428571, 0.14285714],
[ 0. , 0. , 0. , -3.53846154]])
Damit ist die Faktorisierung abgeschlossen. Die aus der Pivotisierung erhaltene Reihenfolge der Zeilen ist gegeben durch
[23]:
ind
[23]:
array([2, 3, 0, 1])
Damit erhalten wir eine obere Dreiecksmatrix:
[24]:
R[ind,:]
[24]:
array([[ 6. , 1. , -1. , 6. ],
[ 0. , -4.66666667, 3.66666667, 3. ],
[ 0. , 0. , -3.71428571, 0.14285714],
[ 0. , 0. , 0. , -3.53846154]])
und für die untere Dreiecksmatrix folgt nach Addition der Einheitsmatrix:
[25]:
L[ind,:]+= np.eye(n)
Kontrolle¶
[26]:
(L[ind]@(R[ind,:]))-A[ind]
[26]:
array([[ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00],
[ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00],
[ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00],
[ 0.00000000e+00, -4.22942105e-17, 0.00000000e+00,
0.00000000e+00]])
Damit das Handling praktisch bleibt, speichern wir \(L,R\) als untere und obere Dreiecksmatrix:
[27]:
L = L[ind,:]
R = R[ind,:]
Die Permutationsmatrix kann einfach erstellt werden:
[28]:
P = np.eye(n)[ind]
P
[28]:
array([[0., 0., 1., 0.],
[0., 0., 0., 1.],
[1., 0., 0., 0.],
[0., 1., 0., 0.]])
Damit ist das Resultat gegeben durch:
[29]:
P, L, R
[29]:
(array([[0., 0., 1., 0.],
[0., 0., 0., 1.],
[1., 0., 0., 0.],
[0., 1., 0., 0.]]),
array([[ 1. , 0. , 0. , 0. ],
[-0.33333333, 1. , 0. , 0. ],
[ 0.33333333, 0.28571429, 1. , 0. ],
[ 0.66666667, 0.14285714, 0.76923077, 1. ]]),
array([[ 6. , 1. , -1. , 6. ],
[ 0. , -4.66666667, 3.66666667, 3. ],
[ 0. , 0. , -3.71428571, 0.14285714],
[ 0. , 0. , 0. , -3.53846154]]))
Lösen der linearen Gleichung¶
Das Lösen der linearen Gleichung
erfolgt nun mit Hilfe der unteren und oberen Dreickmatrix durch Vorwärts- bzw. Rückwärtseinsetzen (Algorithmus 3.1).
Es gilt
somit folgt
Mit \(z = R\cdot x\) erhalten wir für (*) die beiden linearen Gleichungen
Rechte Seite für das Problem:
[30]:
b = np.array([1, -8, -16, -12])
b
[30]:
array([ 1, -8, -16, -12])
Vorwärtseinsetzen¶
[31]:
z = np.zeros(n)
z[0] = b[ind[0]]
for i in range(1,n):
z[i] = b[ind[i]] - L[i,:i]@z[:i]
z
[31]:
array([-16. , -17.33333333, 11.28571429, -3.53846154])
Vorsicht: Das Vorgehen setzt eine normierte Matrix \(L\) voraus!
Rückwärtseinsetzen¶
[32]:
x = np.zeros(n)
x[n-1] = z[n-1]/R[n-1,n-1]
for i in range(2,-1,-1):
x[i] = (z[i]-R[i,i+1:]@x[i+1:])/R[i,i]
x
[32]:
array([-4.5, 2. , -3. , 1. ])
Bemerkung: Die Reihenfolge im Vektor z entspricht nicht der Reihenfolge der Zeilen! Da der Vektor nur ein Zwischenresultat enthält, hat dies keinen Einfluss.
Residuum: \(\|A\cdot x -b\|\)
[33]:
np.linalg.norm(A@x-b)
[33]:
1.831026719408895e-15
Python lu Methode in Scipy¶
In Python steht uns die LR-Zerlegung in scipy der Erweiterung von numpy zur Verfügung:
[34]:
from scipy.linalg import lu
[35]:
help(lu)
Help on function lu in module scipy.linalg.decomp_lu:
lu(a, permute_l=False, overwrite_a=False, check_finite=True)
Compute pivoted LU decomposition of a matrix.
The decomposition is::
A = P L U
where P is a permutation matrix, L lower triangular with unit
diagonal elements, and U upper triangular.
Parameters
----------
a : (M, N) array_like
Array to decompose
permute_l : bool, optional
Perform the multiplication P*L (Default: do not permute)
overwrite_a : bool, optional
Whether to overwrite data in a (may improve performance)
check_finite : bool, optional
Whether to check that the input matrix contains only finite numbers.
Disabling may give a performance gain, but may result in problems
(crashes, non-termination) if the inputs do contain infinities or NaNs.
Returns
-------
**(If permute_l == False)**
p : (M, M) ndarray
Permutation matrix
l : (M, K) ndarray
Lower triangular or trapezoidal matrix with unit diagonal.
K = min(M, N)
u : (K, N) ndarray
Upper triangular or trapezoidal matrix
**(If permute_l == True)**
pl : (M, K) ndarray
Permuted L matrix.
K = min(M, N)
u : (K, N) ndarray
Upper triangular or trapezoidal matrix
Notes
-----
This is a LU factorization routine written for SciPy.
Examples
--------
>>> from scipy.linalg import lu
>>> A = np.array([[2, 5, 8, 7], [5, 2, 2, 8], [7, 5, 6, 6], [5, 4, 4, 8]])
>>> p, l, u = lu(A)
>>> np.allclose(A - p @ l @ u, np.zeros((4, 4)))
True
Wir erhalten hier exakt obiges Resultat mit der Permutationsmatrix:
[36]:
Ps, Ls, Rs = lu(A)
Ps, Ls, Rs
[36]:
(array([[0., 0., 1., 0.],
[0., 0., 0., 1.],
[1., 0., 0., 0.],
[0., 1., 0., 0.]]),
array([[ 1. , 0. , 0. , 0. ],
[-0.33333333, 1. , 0. , 0. ],
[ 0.33333333, 0.28571429, 1. , 0. ],
[ 0.66666667, 0.14285714, 0.76923077, 1. ]]),
array([[ 6. , 1. , -1. , 6. ],
[ 0. , -4.66666667, 3.66666667, 3. ],
[ 0. , 0. , -3.71428571, 0.14285714],
[ 0. , 0. , 0. , -3.53846154]]))
[37]:
ind
[37]:
array([2, 3, 0, 1])
[ ]: