Added faster method for inverting symmetric matrices

Fixes #55.
The implementation was taken from here:
https://stackoverflow.com/a/58719188
with some minor improvements.

fitk/derivatives.py

@@ -19,7 +19,13 @@
 
 # first party imports
 from fitk.tensors import FisherMatrix
-from fitk.utilities import P, ValidationError, find_diff_weights, is_iterable
+from fitk.utilities import (
+    P,
+    ValidationError,
+    fast_positive_definite_inverse,
+    find_diff_weights,
+    is_iterable,
+)
 
 
 def _zero_out(array, threshold: float):
@@ -636,7 +642,7 @@ def fisher_matrix(
                 *zip(names, fiducials), **kwargs_covariance
             )
 
-        inverse_covariance_matrix = np.linalg.inv(covariance_matrix)
+        inverse_covariance_matrix = fast_positive_definite_inverse(covariance_matrix)
 
         covariance_shape = np.shape(inverse_covariance_matrix)
 

fitk/utilities.py

@@ -12,11 +12,13 @@
 import math
 from collections.abc import Collection, Sequence
 from dataclasses import dataclass
+from functools import lru_cache
 from math import factorial
 from typing import Optional, Union
 
 # third party imports
 import numpy as np
+from scipy.linalg import lapack
 
 
 @dataclass
@@ -396,3 +398,33 @@ def math_mode(
         raise TypeError(err) from err
 
     return [f"${_}$" for _ in arg]
+
+
+@lru_cache(maxsize=None)
+def _get_inds(size: int):
+    return np.tri(size, k=-1, dtype=bool)
+
+
+def upper_triangular_to_symmetric(upper_triangular):
+    r"""
+    Convert an upper triangular matrix to a full, symmetric matrix
+    """
+    size = upper_triangular.shape[0]
+    inds = _get_inds(size)
+    upper_triangular[inds] = upper_triangular.T[inds]
+
+
+def fast_positive_definite_inverse(matrix):
+    r"""
+    Compute the fast inverse of a positive-definite symmetric NxN matrix
+    """
+    cholesky, info = lapack.dpotrf(matrix)
+    if info != 0:
+        raise ValueError(f"dpotrf failed on input {matrix}")
+
+    inv, info = lapack.dpotri(cholesky)
+    if info != 0:
+        raise ValueError(f"dpotri failed on input {cholesky}")
+
+    upper_triangular_to_symmetric(inv)
+    return inv