Cholesky Algorithm

choleskyAlgorithm.py

# ------------------------------------------------------------------------
# The following Python code is implemented by Professor Terje Haukaas at
# the University of British Columbia in Vancouver, Canada. It is made
# freely available online at terje.civil.ubc.ca together with notes,
# examples, and additional Python code. Please be cautious when using
# this code; it may contain bugs and comes without any form of warranty.
# Please see the Programming note for how to get started, and notice
# that you must copy certain functions into the file terjesfunctions.py
#
# The Cholesky algorithm implemented below is named after Andre-Louis
# Cholesky (1857-1918) and takes a symmetric and positive definite matrix
# A as input. The objective is to decompose A into a lower and upper
# triangular matrix where one equals the transpose of the other:
#
#                     A = L * L^T
#
# The Cholesky algorithm is used in solving linear systems of equations
# and also in reliability analysis when random variables are transformed
# into standard normal uncorrelated variables. See the Mathematics note
# and the document on MVFOSM, FORM, and SORM for further details.
#
# ------------------------------------------------------------------------

import numpy as np

def choleskyDecomposition(A):

    # Get the size of the input matrix
    n = A.shape[0]

    # Create output matrices of the same size
    L = np.zeros((n, n))
    inverseL = np.zeros((n, n))

    # Loop over rows and columns
    for i in range(n):

        for j in range(n):

            if (i < j):

                L.itemset((i, j), 0.0)

            elif (i == j):

                sum = 0.0

                for k in range(i):

                    sum += L[i, k] * L[i, k]

                if (A[i, j] - sum > 0.0):

                    L.itemset((i, j), np.sqrt(A[i, j] - sum))

                else:

                    print("The Cholesky decomposition failed at Operation 1")
                    return [L, inverseL]

            elif (i > j):

                sum = 0.0
                for k in range(j):

                    sum += L[i, k] * L[j, k]

                if (L[j, j] != 0.0):

                    L.itemset((i, j), ((A[i, j] - sum) / L[j, j]))


                else:

                    print("The Cholesky decomposition failed at Operation 2")
                    return [L, inverseL]

    # Compute the inverse Cholesky matrix
    for i in range(n):

        for j in range(n-1, -1, -1):

            if (i < j):

                inverseL.itemset((i, j), 0.0)

            elif (i == j):

                if (L[i, j] != 0.0):

                    inverseL.itemset((i, j), 1.0 / L[i, j])

                else:

                    print("The Cholesky inversion failed at Operation 3")
                    return [L, inverseL]

            elif (i > j):

                sum = 0.0
                for k in range(j, i, 1):

                    sum += L[i, k] * inverseL[k, j]

                if (L[i, i] != 0.0):

                    inverseL.itemset((i, j), (-sum/L[i, i]))

                else:

                    print("The Cholesky inversion failed at Operation 4")
                    return [L, inverseL]

    return [L, inverseL]

# Test of the Cholesky algorithm
A = np.array([(2.0, 1.0),
              (1.0, 3.0)])
[L, inverseL] = choleskyDecomposition(A)
print("The Cholesky decomposition is:", L)
print("  ")
print("The inverse is:", inverseL)