Merge pull request #338 from hugohadfield/hitzer-inverse

Add the closed form inverse for total dimension <= 5

clifford/_layout.py

@@ -458,6 +458,18 @@ def dual_func(Xval):
                 return gmt_func(Xval, Iinv)
             return dual_func
 
+    @_cached_property
+    def _grade_invol(self):
+        """
+        Generates the grade involution function
+        """
+        signs = np.power(-1, self._basis_blade_order.grades)
+        @_numba_utils.njit
+        def grade_inv_func(mv):
+            newValue = signs * mv.value
+            return self.MultiVector(newValue)
+        return grade_inv_func
+
     @_cached_property
     def vee_func(self):
         """
@@ -568,6 +580,50 @@ def comp_func(Xval):
             return Yval
         return comp_func
 
+    @_cached_property
+    def _hitzer_inverse(self):
+        """
+        Performs the inversion operation as described in the paper :cite:`Hitzer_Sangwine_2017`
+        """
+        tot = np.sum(self.sig != 0)
+        @_numba_utils.njit
+        def hitzer_inverse(operand):
+            if tot == 0:
+                numerator = operand.layout.scalar
+            elif tot == 1:
+                # Equation 4.3
+                mv_invol = operand.gradeInvol()
+                numerator = mv_invol
+            elif tot == 2:
+                # Equation 5.5
+                mv_conj = operand.conjugate()
+                numerator = mv_conj
+            elif tot == 3:
+                # Equation 6.5  without the rearrangement from 6.4
+                mv_conj = operand.conjugate()
+                mv_mul_mv_conj = operand * mv_conj
+                numerator = (mv_conj * ~mv_mul_mv_conj)
+            elif tot == 4:
+                # Equation 7.7
+                mv_conj = operand.conjugate()
+                mv_mul_mv_conj = operand * mv_conj
+                numerator = mv_conj * (mv_mul_mv_conj - 2 * mv_mul_mv_conj(3, 4))
+            elif tot == 5:
+                # Equation 8.22 without the rearrangement from 8.21
+                mv_conj = operand.conjugate()
+                mv_mul_mv_conj = operand * mv_conj
+                combo_op = mv_conj * ~mv_mul_mv_conj
+                mv_combo_op = operand * combo_op
+                numerator = combo_op * (mv_combo_op - 2 * mv_combo_op(1, 4))
+            else:
+                raise NotImplementedError(
+                    'Closed form inverses for algebras with more than 5 dimensions are not implemented')
+            denominator = (operand * numerator).value[0]
+            if denominator == 0:
+                raise ValueError('Multivector has no inverse')
+            return numerator / denominator
+        return hitzer_inverse
+
     @_cached_property
     def gmt_func(self):
         return get_mult_function(self.gmt, self._basis_blade_order.grades)

clifford/_multivector.py

@@ -716,6 +716,9 @@ def normal(self) -> 'MultiVector':
 
         return self / abs(self)
 
+    def hitzer_inverse(self):
+        return self.layout._hitzer_inverse(self)
+
     def leftLaInv(self) -> 'MultiVector':
         """Return left-inverse using a computational linear algebra method
         proposed by Christian Perwass.
@@ -729,14 +732,17 @@ def _pick_inv(self, fallback):
         ----------
         fallback : bool, optional
             If `None`, perform no checks on whether normal inv is appropriate.
-            If `True`, fallback to a linalg approach if necessary.
+            If `True`, fallback to a Hitzer and Sangwine's method :cite:`Hitzer_Sangwine_2017` if possible and a linalg approach if not.
             If `False`, raise an error if normal inv is not appropriate.
         """
         Madjoint = ~self
         MadjointM = (Madjoint * self)
         if fallback is not None and not MadjointM.isScalar():
             if fallback:
-                return self.leftLaInv()
+                try:
+                    return self.hitzer_inverse()
+                except NotImplementedError:
+                    return self.leftLaInv()
             else:
                 raise ValueError("no inverse exists for this multivector")
 
@@ -801,12 +807,7 @@ def gradeInvol(self) -> 'MultiVector':
             M^* = \sum_{i=0}^{\text{dims}}
                   {(-1)^i \left<M\right>_i}
         """
-
-        signs = np.power(-1, self.layout._basis_blade_order.grades)
-
-        newValue = signs * self.value
-
-        return self._newMV(newValue)
+        return self.layout._grade_invol(self)
 
     @property
     def even(self) -> 'MultiVector':

clifford/numba/_multivector.py

@@ -337,3 +337,16 @@ def MultiVector___abs__(self):
 @numba.extending.overload_method(MultiVectorType, 'normal')
 def MultiVector_normal(self):
     return MultiVector.normal
+
+
+@numba.extending.overload_method(MultiVectorType, 'gradeInvol')
+def MultiVector_gradeInvol(self):
+    g_func = self.layout_type.obj._grade_invol
+    def impl(self):
+        return g_func(self)
+    return impl
+
+
+@numba.extending.overload_method(MultiVectorType, 'conjugate')
+def MultiVector_conjugate(self):
+    return MultiVector.conjugate

clifford/test/test_clifford.py

@@ -60,7 +60,7 @@ def test_inverse(self, algebra):
             1 / a
         for i in range(10):
             a = randomMV(layout, grades=[0, 1])
-            denominator = float(a(1)**2-a(0)**2)
+            denominator = (a(1)**2)[()]-(a[()]**2)
             if abs(denominator) > 1.e-5:
                 a_inv = (-a(0)/denominator) + ((1./denominator) * a(1))
                 assert abs((a * a_inv)-1.) < 1.e-11

clifford/test/test_multivector_inverse.py

@@ -0,0 +1,21 @@
+import numpy as np
+import pytest
+
+import clifford as cf
+
+
+class TestClosedForm:
+
+    @pytest.mark.parametrize('p, q', [
+        (p, total_dims - p)
+        for total_dims in [1, 2, 3, 4, 5]
+        for p in range(total_dims + 1)
+    ])
+    def test_hitzer_inverse(self, p, q):
+        Ntests = 100
+        layout, blades = cf.Cl(p, q)
+        for i in range(Ntests):
+            mv = layout.randomMV()
+            mv_inv = mv.hitzer_inverse()
+            np.testing.assert_almost_equal((mv * mv_inv).value,
+                                           layout.scalar.value)

docs/refs.bib

@@ -62,3 +62,16 @@ @article{direct-linear-interpolation
     journal = {Advances in Applied Clifford Algebras},
     doi = {10.1007/s00006-019-1003-y}
 }
+
+@article{Hitzer_Sangwine_2017,
+    title = {Multivector and multivector matrix inverses in real Clifford algebras},
+    volume = {311},
+    ISSN = {0096-3003},
+    DOI = {10.1016/j.amc.2017.05.027},
+    abstractNote = {We show how to compute the inverse of multivectors in finite dimensional real Clifford algebras Cl(p, q). For algebras over vector spaces of fewer than six dimensions, we provide explicit formulae for discriminating between divisors of zero and invertible multivectors, and for the computation of the inverse of a general invertible multivector. For algebras over vector spaces of dimension six or higher, we use isomorphisms between algebras, and between multivectors and matrix representations with multivector elements in Clifford algebras of lower dimension. Towards this end we provide explicit details of how to compute several forms of isomorphism that are essential to invert multivectors in arbitrarily chosen algebras. We also discuss briefly the computation of the inverses of matrices of multivectors by adapting an existing textbook algorithm for matrices to the multivector setting, using the previous results to compute the required inverses of individual multivectors.},
+    journal = {Applied Mathematics and Computation},
+    author = {Hitzer, Eckhard and Sangwine, Stephen},
+    year = {2017},
+    month = {Oct},
+    pages = {375–389}
+}