def cholesky(A): L = [[0.0] * len(A) for _ in range(len(A))] for i, (Ai, Li) in enumerate(zip(A, L)): for j, Lj in enumerate(L[:i+1]): s = sum(Li[k] * Lj[k] for k in range(j)) Li[j] = sqrt(Ai[i] - s) if (i == j) else \ (1.0 / Lj[j] * (Ai[j] - s)) return L