Added faster method for inverting symmetric matrices

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):
@@ -570,6 +576,9 @@ def fisher_matrix(
         ValueError
             if the argument `external_covariance` is not a square matrix
 
+        LinAlgError
+            if the covariance matrix is singular
+
         Notes
         -----
         The `order` parameter is ignored if passed to `D`.
@@ -636,7 +645,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/tensors.py

@@ -29,6 +29,7 @@
     MismatchingValuesError,
     P,
     ParameterNotFoundError,
+    fast_positive_definite_inverse,
     get_index_of_other_array,
     is_iterable,
     is_positive_semidefinite,
@@ -1413,6 +1414,11 @@ def inverse(self):
         array_like : float
             the covariance matrix as a numpy array
 
+        Raises
+        ------
+        LinAlgError
+            if the Fisher matrix is singular
+
         Examples
         --------
         >>> fm = FisherMatrix(np.diag([1, 2, 5]))
@@ -1421,7 +1427,7 @@ def inverse(self):
                [0. , 0.5, 0. ],
                [0. , 0. , 0.2]])
         """
-        return np.linalg.inv(self.values)
+        return fast_positive_definite_inverse(self.values)
 
     def determinant(self):
         """

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.lapack import dpotrf, dpotri  # pylint: disable=no-name-in-module
 
 
 @dataclass
@@ -396,3 +398,29 @@ 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 = dpotrf(matrix)
+    if info != 0:
+        raise np.linalg.LinAlgError(f"dpotrf failed on input {matrix}")
+
+    inv, info = dpotri(cholesky)
+    if info != 0:
+        raise np.linalg.LinAlgError(f"dpotri failed on input {cholesky}")
+
+    upper_triangular_to_symmetric(inv)
+    return inv