Bemerkungen zur Cholesky Zerlegung¶
[1]:
import numpy as np
from scipy.linalg import lu
from numpy.linalg import norm
Implementierung in NumPy und SciPy¶
[2]:
from numpy.linalg import cholesky
Konstruktion einer symmetrischen positiv definiten Matrix \(B = A^T\cdot A\) mit Zufallszahlen für die Matrix \(A\):
[3]:
A = np.random.rand(5,5)
B=A.T@A
Berechnung der Matrix-Faktorisierung
[4]:
L = cholesky(B)
[5]:
L
[5]:
array([[0.62323157, 0. , 0. , 0. , 0. ],
[1.00244892, 0.88547428, 0. , 0. , 0. ],
[0.65200281, 0.83659992, 1.04851317, 0. , 0. ],
[1.25927485, 0.81392201, 0.61344245, 0.7100628 , 0. ],
[1.04752713, 0.55628067, 0.3017882 , 0.17667294, 0.05119053]])
In scipy ist ebenfalls eine Cholesky Zerlegung vorhanden, bei welcher die Funktionalität analog zu Matlab gegeben ist.
[6]:
from scipy.linalg import cholesky as choleskyScipy
[7]:
choleskyScipy(B,lower=True)
[7]:
array([[0.62323157, 0. , 0. , 0. , 0. ],
[1.00244892, 0.88547428, 0. , 0. , 0. ],
[0.65200281, 0.83659992, 1.04851317, 0. , 0. ],
[1.25927485, 0.81392201, 0.61344245, 0.7100628 , 0. ],
[1.04752713, 0.55628067, 0.3017882 , 0.17667294, 0.05119053]])
Darstellung der Cholesky Zerlegung¶
Die Darstellungen sind äquivalent.
Man kann jeweilen die andere leicht berechnen.
[8]:
invsqrtD = np.diag(1/np.diag(L))
Normierte L Matrix:
[9]:
L1 = L@invsqrtD
[10]:
L1
[10]:
array([[1. , 0. , 0. , 0. , 0. ],
[1.60846941, 1. , 0. , 0. , 0. ],
[1.0461646 , 0.94480432, 1. , 0. , 0. ],
[2.0205569 , 0.91919329, 0.58505937, 1. , 0. ],
[1.68079921, 0.62822905, 0.2878249 , 0.24881311, 1. ]])
Diagonalmatrix muss Diagonale im quadrat sein! Folgt direkt aus \(L\cdot L^T\)
[11]:
D = np.diag(np.diag(L)**2)
D
[11]:
array([[0.38841759, 0. , 0. , 0. , 0. ],
[0. , 0.7840647 , 0. , 0. , 0. ],
[0. , 0. , 1.09937986, 0. , 0. ],
[0. , 0. , 0. , 0.50418917, 0. ],
[0. , 0. , 0. , 0. , 0.00262047]])
[12]:
L1@D@L1.T-B
[12]:
array([[-5.55111512e-17, 0.00000000e+00, 5.55111512e-17,
2.22044605e-16, 1.11022302e-16],
[ 0.00000000e+00, -2.22044605e-16, 0.00000000e+00,
2.22044605e-16, 0.00000000e+00],
[ 5.55111512e-17, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 2.22044605e-16],
[ 2.22044605e-16, 4.44089210e-16, 0.00000000e+00,
8.88178420e-16, 8.88178420e-16],
[ 1.11022302e-16, 2.22044605e-16, 2.22044605e-16,
8.88178420e-16, 4.44089210e-16]])
[13]:
norm(L1@D@L1.T-B,np.inf)
[13]:
np.float64(2.4424906541753444e-15)
Rechenzeit¶
[14]:
A = np.random.rand(5000,5000)
B = A.T@A
Rechenzeit für die LR-Zerlegung:
[15]:
%timeit lu(B)
937 ms ± 11.2 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
Rechenzeit für die Cholesky-Zerlegung:
NumPy Variante:
[16]:
%timeit cholesky(B)
545 ms ± 5.71 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
SciPy Variante:
[17]:
%timeit choleskyScipy(B)
455 ms ± 1.4 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
Lösen einer linearen Gleichung¶
[18]:
from scipy.linalg import solve_triangular
Als Beispiel betrachten wir das lineare Gleichungssystem
[19]:
n = 500
A = np.random.rand(n,n)
B = A.T@A
b = np.random.rand(n)
Darstellung \(B = L\cdot L^T\)
Vorwärtseinsetzen:
Rückwärtseinsetzen:
[20]:
L = cholesky(B)
[21]:
z = solve_triangular(L,b,lower=True) # Achtung: Vorwärtseinsetzen MIT Division durch Diagonale (ist hier der Fall)
x = solve_triangular(L.T,z,lower=False)
Kontrolle:
[22]:
norm(B@x-b,np.inf)
[22]:
np.float64(1.2863987652877995e-10)
Darstellung \(B = L\cdot \underbrace{D\cdot L^T}_{=R} = L\cdot R\)
Vorwärtseinsetzen:
Rückwärtseinsetzen:
[23]:
D = np.diag(np.diag(L)**2)
L1 = L@np.diag(1/np.diag(L))
[24]:
z = solve_triangular(L1,b,lower=True) # Achtung: Vorwärtseinsetzen MIT Division durch Diagonale (ist hier der Fall)
x = solve_triangular(D@L1.T,z,lower=False)
Kontrolle:
[25]:
norm(B@x-b,np.inf)
[25]:
np.float64(9.953604607204625e-11)