LU decomposition

lu.py

__author__ = 'Giorgi Tkeshelashvili'
from pprint import pprint

def mult_matrix(A, B):
    """ Multiply square matrices A and B of the same dimension (nxn)"""
    # Converts N into a list of tuples of columns                                                                                                                                                                                                      
    col_tuples = zip(*B)
    # Multiply A and B
    return [[sum(ea*eb for ea,eb in zip(a,b)) for b in col_tuples] for a in A]

def create_pivot_matrix(M):
    """Creates the permutation matrix for M."""
    n = len(M)
    # Create identity matrix
    ID = [[float(i == j) for i in range(n)] for j in range(n)]
    # Place largest element of each column of matrix M on the diagonal
    for j in range(n):
        row = max(range(j, n), key=lambda i: abs(M[i][j]))
        if j != row:
            #swap the rows
            ID[j], ID[row] = ID[row], ID[j]
    return ID

def lu(A):
    """Decomposes any square (nxn) matrix A into PA=LU and returns P,L and U"""
    n = len(A)
    # initialize L and U to zero matrices
    L = [[0.0] * n for i in range(n)]
    U = [[0.0] * n for i in range(n)]
    # Create permutation matrix of A
    P = create_pivot_matrix(A)
    # Multiply matrix A and its permutation matrix P
    PA = mult_matrix(P, A)
    # Decompose into LU 
    for j in range(n):
        # Diagonal of L is set to 1
        L[j][j] = 1.0
        # Create U matrix
        for i in range(j+1):
            s1 = sum(U[k][j] * L[i][k] for k in range(i))
            U[i][j] = PA[i][j] - s1
        # Check the diagonal of U and see if U is nonsingular
        for i in range(j+1):
            if(U[j][j]==0):
                print "Matrix is Singular"
                # Exit because the matrix is singular
                return
        # Create L matrix
        for i in range(j, n):
            s2 = sum(U[k][j] * L[i][k] for k in range(j))
            L[i][j] = (PA[i][j] - s2) / U[j][j]
    # First Matrix returned is L, second is U and the third is P
    return (L, U, P)

#################################################
# example 1: Matrix a:
#    [4 3
#     6 3]
#
a = [[4,3], [6,3]]
for part in lu(a):
    pprint(part, width=19)
    print
print
# Output for example 1:

# L matrix:
#
#                [[1.0,               0.0],
#                [0.6666666666666666, 1.0]]
#
# U matrix:
#               [[6.0, 3.0],
#                [0.0, 1.0]]
# 
# P matrix:
#               [[0.0, 1.0],
#                [1.0, 0.0]]
#
#
#
#################################################
# example 2: Matrix b: 
#
#       [1,  2,  3
#        2, -4,  6
#        3, -9, -3]
#
b = [[1,2,3], [2,-4,6],[3,-9,-3]]
for part in lu(b):
    pprint(part)
    print
# Output for example 2:

# L matrix:   
#           [[1.0,                0.0,    0.0],
#           [0.3333333333333333, 1.0,    0.0],
#           [0.6666666666666666, 0.4,    1.0]]
#                                            
#
# U matrix:
#           [[3.0,  -9.0,   -3.0],
#            [0.0,   5.0,    4.0],
#            [0.0,   0.0,    6.4]]
# 
# P matrix:
#       [[0.0, 1.0],
#       [1.0, 0.0]]
#
#
#
#################################################
# example 3: Matrix c: Singular Matrix (det(c)=0)
#       2 6
#       1 3 
c = [[2,6],[1,3]]
lu(c)
# Output for example 3:
#    "Matrix is Singular"