Frobenius Matrix

\[\begin{split}L_k := \begin{pmatrix} 1 & \\ & \ddots & & & 0\\ & & 1 \\ & & & 1 \\ & & & -\ell_{k+1,k} & 1\\ & 0 & & \vdots & & \ddots\\ & & & -\ell_{n,k} & 0 & & 1 \end{pmatrix}\end{split}\]

Betrachte die Aufgabe aus dem Notebook Beispiel LU-Zerlegung.

[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]])
A
[2]:
array([[ 2, -1, -3,  3],
       [ 4,  0, -3,  1],
       [ 6,  1, -1,  6],
       [-2, -5,  4,  1]])

1. Spalte

Aus der Gauss-Elimintation folgt für \(L_1\)

[3]:
L1 = np.eye(N)
L1[1:,0] = -np.array([2,3,-1])
L1
[3]:
array([[ 1.,  0.,  0.,  0.],
       [-2.,  1.,  0.,  0.],
       [-3.,  0.,  1.,  0.],
       [ 1.,  0.,  0.,  1.]])
[4]:
L1@A
[4]:
array([[ 2., -1., -3.,  3.],
       [ 0.,  2.,  3., -5.],
       [ 0.,  4.,  8., -3.],
       [ 0., -6.,  1.,  4.]])

Die Inverse der Frobenius Matrix ist gegeben durch

\[\begin{split}L_k^{-1} := \begin{pmatrix} 1 & \\ & \ddots & & & 0\\ & & 1 \\ & & & 1 \\ & & & \ell_{k+1,k} & 1\\ & 0 & & \vdots & & \ddots\\ & & & \ell_{n,k} & 0 & & 1 \end{pmatrix}\end{split}\]

in dem Fall somit:

[5]:
invL1 = np.eye(N)
invL1[1:,0] = np.array([2,3,-1])
invL1
[5]:
array([[ 1.,  0.,  0.,  0.],
       [ 2.,  1.,  0.,  0.],
       [ 3.,  0.,  1.,  0.],
       [-1.,  0.,  0.,  1.]])

Kontrolle:

[6]:
L1@invL1
[6]:
array([[1., 0., 0., 0.],
       [0., 1., 0., 0.],
       [0., 0., 1., 0.],
       [0., 0., 0., 1.]])

2. Spalte

[7]:
L2 = np.eye(N)
L2[2:,1] = -np.array([2,-3])
L2
[7]:
array([[ 1.,  0.,  0.,  0.],
       [ 0.,  1.,  0.,  0.],
       [ 0., -2.,  1.,  0.],
       [ 0.,  3.,  0.,  1.]])

Für die Matrix \(A\) folgt:

[8]:
L2@L1@A
[8]:
array([[  2.,  -1.,  -3.,   3.],
       [  0.,   2.,   3.,  -5.],
       [  0.,   0.,   2.,   7.],
       [  0.,   0.,  10., -11.]])

Für die Inverse folgt wiederum:

[9]:
invL2 = np.eye(N)
invL2[2:,1] = np.array([2,-3])
invL2
[9]:
array([[ 1.,  0.,  0.,  0.],
       [ 0.,  1.,  0.,  0.],
       [ 0.,  2.,  1.,  0.],
       [ 0., -3.,  0.,  1.]])

Kontrolle

[10]:
L2@invL2
[10]:
array([[1., 0., 0., 0.],
       [0., 1., 0., 0.],
       [0., 0., 1., 0.],
       [0., 0., 0., 1.]])

3. Spalte

[11]:
L3 = np.eye(N)
L3[3:,2] = -np.array([5])
L3
[11]:
array([[ 1.,  0.,  0.,  0.],
       [ 0.,  1.,  0.,  0.],
       [ 0.,  0.,  1.,  0.],
       [ 0.,  0., -5.,  1.]])

Für die Matrix \(A\) folgt:

[12]:
L3@L2@L1@A
[12]:
array([[  2.,  -1.,  -3.,   3.],
       [  0.,   2.,   3.,  -5.],
       [  0.,   0.,   2.,   7.],
       [  0.,   0.,   0., -46.]])

Für die Inverse folgt wiederum:

[13]:
invL3 = np.eye(N)
invL3[3:,2] = np.array([5])
invL3
[13]:
array([[1., 0., 0., 0.],
       [0., 1., 0., 0.],
       [0., 0., 1., 0.],
       [0., 0., 5., 1.]])

Kontrolle

[14]:
L3@invL3
[14]:
array([[1., 0., 0., 0.],
       [0., 1., 0., 0.],
       [0., 0., 1., 0.],
       [0., 0., 0., 1.]])

Resultat

\[L_3\cdot L_2\cdot L_1 \cdot A = R\]

Gesucht ist

\[(L_3\cdot L_2\cdot L_1)^{-1}\]

so, dass die Matrix Faktorisierung

\[A = L\cdot R\]

folgt. Das Berechnen der Inversen von \(L_j\) ist einfach (vgl. oben). Für die Inverse gilt

\[(L_3\cdot L_2\cdot L_1)^{-1} = L_1^{-1}\cdot L_2^{-1}\cdot L_3^{-1}\]
[15]:
invL1@invL2@invL3
[15]:
array([[ 1.,  0.,  0.,  0.],
       [ 2.,  1.,  0.,  0.],
       [ 3.,  2.,  1.,  0.],
       [-1., -3.,  5.,  1.]])

Kontrolle:

[16]:
L3@L2@L1@invL1@invL2@invL3
[16]:
array([[1., 0., 0., 0.],
       [0., 1., 0., 0.],
       [0., 0., 1., 0.],
       [0., 0., 0., 1.]])