Adds matrix inverse operation

Author: DrPaulSharp

sasdata/quantities/quantity.py

@@ -52,6 +52,18 @@ def trace(a: Union["Quantity[ArrayLike]", ArrayLike], offset: int = 0, axis1: in
         return np.trace(a, offset, axis1, axis2)
 
 
+def matinv(a: Union["Quantity[ArrayLike]", ArrayLike]):
+    """Find the inverse of a matrix."""
+    if isinstance(a, Quantity):
+        return DerivedQuantity(
+            value=np.linalg.inv(a.value),
+            units=a.units,
+            history=QuantityHistory.apply_operation(MatInv, a.history),
+        )
+    else:
+        return np.linalg.inv(a)
+
+
 def determinant(a: Union["Quantity[ArrayLike]", ArrayLike]):
     """Find the determinant of an array or an array based quantity."""
     if isinstance(a, Quantity):
@@ -1180,6 +1192,32 @@ def _clean_ab(self, a, b):
         return MatMul(a, b)
 
 
+class MatInv(UnaryOperation):
+    """Matrix inversion, using numpy"""
+
+    serialisation_name = "matinv"
+
+    def evaluate(self, variables: dict[int, T]) -> T:
+        return np.linalg.inv(self.a.evaluate(variables))
+
+    def _derivative(self, hash_value: int) -> Operation:
+        return Neg(Matmul(MatMul(MatInv(self.a), self.a._derivative(hash_value)), MatInv(self.a)))
+
+    def _clean(self):
+        clean_a = self.a._clean()
+
+        if isinstance(clean_a, MatInv):
+            # Removes double inversions
+            return clean_a.a
+
+        elif isinstance(clean_a, Determinant):
+            # Inverse of determinant is determinant of inverse
+            return Determinant(MatInv(clean_a.a))
+
+        else:
+            return MatInv(clean_a)
+
+
 class TensorDot(Operation):
     serialisation_name = "tensor_product"
 
@@ -1236,6 +1274,7 @@ def _deserialise(parameters: dict) -> "Operation":
     Determinant,
     Dot,
     MatMul,
+    MatInv,
     TensorDot,
 ]
 

test/quantities/utest_math_operations.py

@@ -4,7 +4,7 @@
 import pytest
 
 from sasdata.quantities import units
-from sasdata.quantities.quantity import NamedQuantity, tensordot, trace, transpose
+from sasdata.quantities.quantity import NamedQuantity, matinv, tensordot, trace, transpose
 
 order_list = [[0, 1, 2, 3], [0, 2, 1], [1, 0], [0, 1], [2, 0, 1], [3, 1, 2, 0]]
 
@@ -57,6 +57,21 @@ def test_trace_axes(matrix, axis1, axis2, expected_trace):
     ).all()
 
 
+@pytest.mark.parametrize(
+    "matrix, expected_inverse",
+    [
+        (np.array([[1]]), np.array([[1]])),
+        (np.array([[-2.0, 1.0], [1.5, -0.5]]), np.array([[1, 2], [3, 4]])),
+    ],
+)
+def test_inverse(matrix, expected_inverse):
+    """Check that the matinv operation correctly inverse for raw data and quantities."""
+    print(matinv(matrix))
+    print(expected_inverse)
+    assert (matinv(matrix) == expected_inverse).all()
+    assert (matinv(NamedQuantity("testmat", matrix, units=units.none)).value == expected_inverse).all()
+
+
 rng_seed = 1979
 tensor_product_with_identity_sizes = (4, 6, 5)
 

test/quantities/utest_operations.py

@@ -18,6 +18,7 @@
     Inv,
     Ln,
     Log,
+    MatInv,
     MatMul,
     Mul,
     MultiplicativeIdentity,
@@ -48,7 +49,7 @@
 )
 
 
-@pytest.fixture(params=[Determinant, Inv, Exp, Ln, Neg, Sin, ArcSin, Cos, ArcCos, Tan, ArcTan, Transpose])
+@pytest.fixture(params=[Determinant, Inv, Exp, Ln, MatInv, Neg, Sin, ArcSin, Cos, ArcCos, Tan, ArcTan, Transpose])
 def unary_operation(request):
     return request.param(x)
 
@@ -232,6 +233,7 @@ def test_derivative(op, result):
     [
         (Neg(Neg(x))),
         (Inv(Inv(x))),
+        (MatInv(MatInv(x))),
     ],
 )
 def test_clean_double_applications(op):