Instantly share code, notes, and snippets.

Embed
What would you like to do?
Implementation of the algorithm of Banachiewicz for Cholesky-Decomposition
#!/usr/bin/python
# -*- coding: utf-8 -*-
def choleskyBanachiewicz(A):
"""
@param A: a nested list which represents a symmetric,
positiv definit n x n matrix.
@return: False if it didn't work, otherwise the matrix G
"""
n = len(A)
# Initialisation of G with A
G = []
for i in xrange(n):
line = []
for j in xrange(n):
line.append(float(A[i][j]))
G.append(line)
# Computation of G
for j in xrange(n):
# doesn't need the diagonals
tmp = (A[j][j] - sum([G[j][k]**2 for k in xrange(0, j)]))
if tmp < 0:
return False
G[j][j] = tmp**0.5
for i in xrange(n):
G[i][j] = 1/G[j][j] * (A[i][j] - sum([G[i][k]*G[j][k] for k in xrange(0, j)]))
return G
if __name__ == "__main__":
print choleskyBanachiewicz([[1,2],[2,1]]) # Nicht positiv definit
print choleskyBanachiewicz([[5,2],[1,1]]) # not Hermitian
print choleskyBanachiewicz([[5,2],[2,1]]) # should be [[sqrt(5), 2/sqrt(5)],[0, 1/sqrt(5)]]
print choleskyBanachiewicz([[1,0,0],[0,1,0],[0,0,1]])
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment