cholesky decomposition floating point error

Question

Implementing the Choleky decomposition in Python without the use of numpy.linalg.* (yes this is an assignment) I encountered a Problem regarding the use of floating point numbers. My Algorithm works fine with regualar integers:

    for i in range(n):
        for k in range(i + 1):
            sum = 0
            for j in range(k):
                sum += L[i, j] * L[k, j]
            if i == k:
                if M[i, i] - sum < 0:
                    raise ValueError("Matrix is non-positive definite")
                L[i, i] = np.sqrt(M[i, i] - sum)
            else:
                if np.isclose(L[k, k] * (M[i, k] - sum), 0):
                    raise ValueError("Matrix is non-positive definite")
                L[i, k] = (L[k, k] * (M[i, k] - sum)) ** (-1)

I tested the matrix for symmetry beforehand; n is the dimension and L becomes the lower triangular Cholesky Factor.

Using random floating point nxn matrices multiplied with their transpose (in order to get a positive definite matrix) both ValueErrors are raised, i.e w/out raising the ValueErrors the L output matrix is partly filled with NaN and inf values. How can I work with floating point numbers in python?

Edit: Minimal Reproducible Example:

M = np.random.randn(2, 2)
M = np.dot(M, M.transpose())
# produces for example [[0.68283219, -0.33497034], [-0.33497034, 0.72113541]]
run_cholesky(M)

Answer 1

Saving M[i, k] in a variable and then subtracting instead of summing up is fixing the problem:

for i in range(n):
    for k in range(i + 1):
        val = M[i, k]
        for j in range(k):
            val -= L[i, j] * L[k, j]
        if i == k:
            if val < 0:
                raise ValueError("Matrix is non-positive definite")
            L[i, k] = np.sqrt(val)
        else:
            if np.isclose(L[k, k], 0):
                raise ValueError("Matrix is non-positive definite")
            L[i, k] = val / L[k, k]