cholesky

import numpy as np
A = np.array([[0.1665, 0.6561, 1.1457, 1.6353],
             [0.6561, 12.9283, 18.7171, 24.5059],
             [1.1457, 18.7171, 63.6745, 79.5721],
             [1.6353, 24.5059, 79.5721, 177.7939]])

def L(A):
    L = np.zeros(A.shape)
    n = A.shape[0]

    for i in range(n):
        for j in range(i + 1):
            tmp_sum = sum(L[i][k] * L[j][k] for k in range(j))

            if (i == j):
                    L[i][j] = (A[i][i] - tmp_sum) ** 0.5
            else:
                    L[i][j] = (1.0 / L[j][j] * (A[i][j] - tmp_sum))
    return L

def inverse_mat(A):
    l = L(A)
    P = np.zeros(A.shape)
    n = A.shape[0]
    for i in range(n):
        P[i][i] = 1 / l[i, i]
        for j in range(i+1, n):
            tmp = 0
            for k in range(j):
                tmp += l[j, k] * P[k][i]
            P[i][j] = - tmp / l[j, j]
    return P.T @ P
ans = inverse_mat(A)
print(ans)