CODEWITHSUNDEEP
pythonmath
user54677
user54677

Reputation: 11

Cholesky decomposition in python, recursive

I am trying to make a recursive program to compute the Cholesky factorization, but the output is not in lower triangular form. How can I change this to compute it correctly?

def cholesky(A):

    n = np.shape(A)[0]
    A[0,0] = math.sqrt(abs(A[0,0]))

    if n == 1:
        return A
    else:
        A[1:,0] /= A[0,0]
        A[1:,1:] -= np.dot(A[1:,0], (A[1:,0]).T)
        cholesky(A[1:,1:])

Upvotes: 1

Views: 2012

Answers (1)

colcarroll
colcarroll

Reputation: 3682

I'm sure there's a prettier way (in particular, numpy.linalg.cholesky), but this looks like it does ok on reasonable test data. I follow the notation from the wikipedia article, and use their example as the test data.

import numpy as np

def cholesky_reduce(A):
    pivot = A[0, 0]
    b = np.mat(A[1:, 0])
    B = A[1:, 1:]
    return B - (b.T * b) / pivot


def L(A):
    n = A.shape[0]
    if n == 1:
        return np.sqrt(A)
    b = np.mat(A[1:, 0])
    pivot = np.sqrt(A[0, 0])
    return np.bmat([
        [np.mat(pivot), np.zeros((1, n - 1))],
        [b.T / pivot, L(cholesky_reduce(A))]
    ])


def __main():
    A = np.array([[4, 12, -16], [12, 37, -43], [-16, -43, 98]])
    print(L(A))


if __name__ == '__main__':
    __main()

Upvotes: 1

Related Questions