Implementing the choice of the matrix inversion algorithm for the user (#51)

* Added a block inversion algorithm for numerical case

* Matrix inversion algorithm is now configurable

* Preparing for v1.0.6 beta release

Author: jodemaey

documentation/source/conf.py

@@ -19,7 +19,7 @@
 project = 'LayerCake'
 copyright = '2025-2026, Jonathan Demaeyer and Oisín Hamilton'
 author = 'Jonathan Demaeyer and Oisín Hamilton'
-release = 'v1.0.5-beta'
+release = 'v1.0.6-beta'
 version = release
 
 # -- General configuration ---------------------------------------------------

layercake/__init__.py

@@ -13,4 +13,4 @@
            'OperatorTerm', 'ProductOfTerms', 'AdditionOfTerms', 'LinearTerm', 'ConstantTerm', 'Equation', 'Laplacian', 'D',
            'Layer', 'Cake']
 
-__version__ = '1.0.5b0'
+__version__ = '1.0.6b0'

layercake/bakery/cake.py

@@ -55,6 +55,7 @@ def __init__(self):
         self.layers = list()
         self._lhs_inversion_in_layer = True
         self._simplify_after_LHS_inversion = True
+        self._matrix_inv = np.linalg.inv
 
     def add_layer(self, layer):
         """Add a layer object to the cake.
@@ -65,6 +66,7 @@ def add_layer(self, layer):
             Layer object to add to the cake.
         """
         layer._cake_order = len(self.layers)
+        layer._matrix_inv = self._matrix_inv
         self.layers.append(layer)
         layer._cake = self
         for equation in layer.equations:
@@ -262,7 +264,7 @@ def tensor(self):
             if not self._lhs_inversion_in_layer:
                 try:
                     lhs_mat_inverted = np.zeros((self.ndim + 1, self.ndim + 1))
-                    lhs_mat_inverted[1:, 1:] = np.linalg.inv(lhs_mat.todense()[1:, 1:])
+                    lhs_mat_inverted[1:, 1:] = self._matrix_inv(lhs_mat.todense()[1:, 1:])
                     tensor = sp.COO(np.tensordot(lhs_mat_inverted, tensor.to_coo(), 1))
                 except LinAlgError:
                     raise LinAlgError(f'The left-hand side of the cake is not invertible with the provided basis.')

layercake/bakery/layers.py

@@ -63,6 +63,7 @@ def __init__(self, name=''):
         self._lhs_inverted = False
         self._lhs_mat = None
         self._simplify_after_LHS_inversion = True
+        self._matrix_inv = np.linalg.inv
 
     @property
     def _cake_first_index(self):
@@ -298,7 +299,7 @@ def compute_tensor(self, numerical=True, compute_inner_products=False, compute_i
                                 self._lhs_mat[lhs_order:lhs_order + ndim, ofield_order:ofield_order + ondim] + lhs_term.inner_products.todense()
                 else:
                     try:
-                        lhs_mat_inverted[lhs_order:lhs_order + ndim, lhs_order:lhs_order + ndim] = np.linalg.inv(eq.lhs_inner_products_addition.todense())
+                        lhs_mat_inverted[lhs_order:lhs_order + ndim, lhs_order:lhs_order + ndim] = self._matrix_inv(eq.lhs_inner_products_addition.todense())
                         self._lhs_inverted = True
                     except LinAlgError:
                         raise LinAlgError(f'The left-hand side of the equation {eq} is not invertible with the provided basis.')
@@ -345,7 +346,7 @@ def compute_tensor(self, numerical=True, compute_inner_products=False, compute_i
                 lhs_order += ndim
             if self._lhs_inversion and not self._lhs_inverted:
                 try:
-                    lhs_mat_inverted[1:, 1:] = np.linalg.inv(self._lhs_mat.todense()[1:, 1:])
+                    lhs_mat_inverted[1:, 1:] = self._matrix_inv(self._lhs_mat.todense()[1:, 1:])
                     self._lhs_inverted = True
                 except LinAlgError:
                     raise LinAlgError(f'The left-hand side of the layer {self} is not invertible with the provided basis.')

layercake/utils/matrix.py

@@ -9,67 +9,105 @@
 """
 import warnings
 
+import numpy as np
+from numpy.linalg import LinAlgError
+
 from sympy import MutableSparseMatrix
 from sympy.matrices.exceptions import NonInvertibleMatrixError
 
 def block_matrix_inverse(P, blocks_extent, simplify=True):
-    """Function to invert a symbolic matrix :math:`P` divided by blocks.
+    """Function to invert a symbolic or numeric matrix :math:`P` divided by blocks.
 
     Parameters
     ----------
-    P: ~sympy.matrices.immutable.ImmutableSparseMatrix or ~sympy.matrices.mutable.MutableSparseMatrix
+    P: ~sympy.matrices.immutable.ImmutableSparseMatrix or ~sympy.matrices.mutable.MutableSparseMatrix or ~numpy.ndarray
         The block matrix to invert.
     blocks_extent: list(tuple)
         The extent of each block, as a list of 2-tuple.
     simplify: bool, optional
         Try to simplify the inverse.
+        Only applies to symbolic matrices.
         Default to `True`.
 
     Warnings
     --------
-    If the `simplify` argument is set to `False`, this function will not check if the
+    In the symbolic case, if the `simplify` argument is set to `False`, this function will not check if the
     matrix :math:`P` is invertible !
 
     """
-    if simplify:
-        if P.det().simplify() == 0:
-            raise NonInvertibleMatrixError('block_matrix_inverse: The provided matrix is not invertible.')
-    else:
-        warnings.warn(f'Inverting a symbolic matrix without checking that it is invertible. '
-                      'Be cautious about the result.')
-    be = blocks_extent.copy()
-    PP = P
-    B_list = list()
-    C_list = list()
-    Am1_list = list()
-    while len(be) >= 2:
-        rest_block_extent = (be[1][0], be[-1][1])
-        Ainv, B, C, D, PPP = _block_matrix_inverse_2x2(PP, (be[0], rest_block_extent))
-
-        B_list.append(B)
-        C_list.append(C)
-        Am1_list.append(Ainv)
-
-        PP = PPP
-        ret = be[0][1]
-        be = [(p[0]-ret, p[1]-ret) for p in be[1:]]
-
-    PP = PP.adjugate() / PP.det()
-    for Am1, B, C in zip(Am1_list[::-1], B_list[::-1], C_list[::-1]):
-        Ashape = Am1.shape[0]
-        Pshape = PP.shape[0]
-        mat_shape = Ashape + Pshape
-        Pp1 = MutableSparseMatrix(mat_shape, mat_shape, {})
-        Pp1[slice(0, Ashape), slice(0, Ashape)] = Am1 + Am1 @ B @ PP @ C @ Am1
-        Pp1[slice(0, Ashape), slice(Ashape, mat_shape)] = - Am1 @ B @ PP
-        Pp1[slice(Ashape, mat_shape), slice(0, Ashape)] = - PP @ C @ Am1
-        Pp1[slice(Ashape, mat_shape), slice(Ashape, mat_shape)] = PP
-        PP = Pp1
-
-    if simplify:
-        return PP.simplify()
-    else:
+    if isinstance(P, np.ndarray):
+        if np.linalg.det(P) == 0:
+            raise LinAlgError('block_matrix_inverse: The provided matrix is not invertible.')
+        be = blocks_extent.copy()
+        PP = P
+        B_list = list()
+        C_list = list()
+        Am1_list = list()
+        while len(be) >= 2:
+            rest_block_extent = (be[1][0], be[-1][1])
+            Ainv, B, C, D, PPP = _block_matrix_inverse_2x2_num(PP, (be[0], rest_block_extent))
+
+            B_list.append(B)
+            C_list.append(C)
+            Am1_list.append(Ainv)
+
+            PP = PPP
+            ret = be[0][1]
+            be = [(p[0] - ret, p[1] - ret) for p in be[1:]]
+
+        PP = np.linalg.inv(PP)
+        for Am1, B, C in zip(Am1_list[::-1], B_list[::-1], C_list[::-1]):
+            Ashape = Am1.shape[0]
+            Pshape = PP.shape[0]
+            mat_shape = Ashape + Pshape
+            Pp1 = np.zeros((mat_shape, mat_shape))
+            Pp1[slice(0, Ashape), slice(0, Ashape)] = Am1 + Am1 @ B @ PP @ C @ Am1
+            Pp1[slice(0, Ashape), slice(Ashape, mat_shape)] = - Am1 @ B @ PP
+            Pp1[slice(Ashape, mat_shape), slice(0, Ashape)] = - PP @ C @ Am1
+            Pp1[slice(Ashape, mat_shape), slice(Ashape, mat_shape)] = PP
+            PP = Pp1
+
         return PP
+    else:
+        if simplify:
+            if P.det().simplify() == 0:
+                raise NonInvertibleMatrixError('block_matrix_inverse: The provided matrix is not invertible.')
+        else:
+            warnings.warn(f'Inverting a symbolic matrix without checking that it is invertible. '
+                          'Be cautious about the result.')
+        be = blocks_extent.copy()
+        PP = P
+        B_list = list()
+        C_list = list()
+        Am1_list = list()
+        while len(be) >= 2:
+            rest_block_extent = (be[1][0], be[-1][1])
+            Ainv, B, C, D, PPP = _block_matrix_inverse_2x2(PP, (be[0], rest_block_extent))
+
+            B_list.append(B)
+            C_list.append(C)
+            Am1_list.append(Ainv)
+
+            PP = PPP
+            ret = be[0][1]
+            be = [(p[0]-ret, p[1]-ret) for p in be[1:]]
+
+        PP = PP.adjugate() / PP.det()
+        for Am1, B, C in zip(Am1_list[::-1], B_list[::-1], C_list[::-1]):
+            Ashape = Am1.shape[0]
+            Pshape = PP.shape[0]
+            mat_shape = Ashape + Pshape
+            Pp1 = MutableSparseMatrix(mat_shape, mat_shape, {})
+            Pp1[slice(0, Ashape), slice(0, Ashape)] = Am1 + Am1 @ B @ PP @ C @ Am1
+            Pp1[slice(0, Ashape), slice(Ashape, mat_shape)] = - Am1 @ B @ PP
+            Pp1[slice(Ashape, mat_shape), slice(0, Ashape)] = - PP @ C @ Am1
+            Pp1[slice(Ashape, mat_shape), slice(Ashape, mat_shape)] = PP
+            PP = Pp1
+
+        if simplify:
+            return PP.simplify()
+        else:
+            return PP
 
 def _block_matrix_inverse_2x2(P, blocks_extent):
     be = blocks_extent
@@ -81,6 +119,16 @@ def _block_matrix_inverse_2x2(P, blocks_extent):
 
     return Ainv, B, C, D, D - C @ Ainv @ B
 
+def _block_matrix_inverse_2x2_num(P, blocks_extent):
+    be = blocks_extent
+    A = P[slice(*be[0]), slice(*be[0])]
+    B = P[slice(*be[0]), slice(*be[1])]
+    C = P[slice(*be[1]), slice(*be[0])]
+    D = P[slice(*be[1]), slice(*be[1])]
+    Ainv = np.linalg.inv(A)
+
+    return Ainv, B, C, D, D - C @ Ainv @ B
+
 
 if __name__ == '__main__':
     A = MutableSparseMatrix(3, 3, {})
@@ -102,3 +150,22 @@ def _block_matrix_inverse_2x2(P, blocks_extent):
     bl = [(0, 3), (3, 6), (6, 9)]
     Ginv = block_matrix_inverse(G, bl)
 
+    A = np.zeros((3, 3))
+    A[0,0] = 1
+    A[1,1] = 3
+    A[2,2] = 2
+    C = A /2
+    D = 3 * A
+    G =  np.zeros((9, 9))
+    G[:3, :3] = A
+    G[:3, 3:6] = -C
+    G[:3, 6:] = D
+    G[3:6, :3] = C
+    G[3:6, 3:6] = C
+    G[3:6, 6:] = - D
+    G[6:, :3] = - C
+    G[6:, 3:6] = A
+    G[6:, 6:] = D
+    bl = [(0, 3), (3, 6), (6, 9)]
+    Ginv_num = block_matrix_inverse(G, bl)
+

pyproject.toml

@@ -7,7 +7,7 @@ requires = ["setuptools", "wheel"]
 [project]
 requires-python = ">=3.10"
 name = "layercake_model"
-version = "1.0.5-beta"
+version = "1.0.6-beta"
 description = "A framework to design systems of partial differential equations (PDEs), and convert them to ordinary differential equations (ODEs) via Galerkin-type expansions. "
 readme = "README.md"
 authors = [