Bivector and rotor splits for arbitrary GAs

clifford/_layout.py

@@ -17,7 +17,6 @@
 )
 from . import _numba_utils
 from .io import read_ga_file
-from . import _settings
 from ._multivector import MultiVector
 from ._layout_helpers import (
     BasisBladeOrder, BasisVectorIds, canonical_reordering_sign_euclidean
@@ -685,8 +684,6 @@ def inv_func(self):
         @_numba_utils.njit
         def leftLaInvJIT(value):
             intermed = _numba_val_get_left_gmt_matrix(value, k_list, l_list, m_list, mult_table_vals, n_dims)
-            if abs(np.linalg.det(intermed)) < _settings._eps:
-                raise ValueError("multivector has no left-inverse")
             sol = np.linalg.solve(intermed, identity.astype(intermed.dtype))
             return sol
 

clifford/invariant_decomposition.py

@@ -0,0 +1,184 @@
+"""
+.. currentmodule:: clifford.invariant_decomposition
+
+=====================================================
+invariant_decomposition (:mod:`clifford.invariant_decomposition`)
+=====================================================
+
+.. versionadded:: 1.4.0
+
+
+This file implements the invariant decomposition (aka bivector split) of bivectors into
+mutually commuting orthogonal simple bivectors, based on the method of :cite:`roelfs2021thesis`, chapter 6.
+
+The invariant decomposition enables closed form exponentials and logarithms, and the factorization of
+rotors into simple rotors.
+
+Example usage::
+
+>>> B = 1*e12 + 2*e34
+>>> bivector_split(B)
+[1^e12, 2^e34]
+
+Implemented functions
+---------------------
+
+.. autofunction:: bivector_split
+.. autofunction:: rotor_split
+.. autofunction:: exp
+.. autofunction:: log
+
+
+Helper functions
+----------------
+
+.. autofunction:: _bivector_split
+.. autofunction:: single_split
+
+"""
+import math
+from functools import reduce
+
+import numpy as np
+
+from ._settings import _eps
+from . import _numba_utils
+
+
+@_numba_utils.njit(cache=True)
+def _single_split_even_values(Wm_array, li, r):
+    ND = np.zeros((2, Wm_array[0, :].shape[0]), dtype=np.complex_)
+    for i in range(0, Wm_array.shape[0]//2+1):
+        ND[0, :] += Wm_array[2*i, :] * li**(r - i)
+    for i in range(0, Wm_array.shape[0]//2):
+        ND[1, :] += Wm_array[2*i+1, :] * li**(r - i - 1)
+    return ND
+
+
+@_numba_utils.njit(cache=True)
+def _single_split_odd_values(Wm_array, li, r):
+    ND = np.zeros((2, Wm_array[0, :].shape[0]), dtype=np.complex_)
+    for i in range(0, Wm_array.shape[0]//2):
+        ND[0, :] += Wm_array[2 * i + 1, :] * li ** (r - i)
+        ND[1, :] += Wm_array[2*i, :] * li**(r - i)
+    return ND
+
+
+def single_split_even(Wm, li, r):
+    """Helper function to compute a single split for a given set of W_m and
+    eigenvalue lambda_i, when the total number of terms in the split is even.
+    """
+    Wm_array = np.array([W.value for W in Wm])
+    ND = _single_split_even_values(Wm_array, li, r)
+    N = Wm[0].layout.MultiVector(ND[0, :])
+    D = Wm[0].layout.MultiVector(ND[1, :])
+    return N*D.leftLaInv()
+
+
+def single_split_odd(Wm, li, r):
+    """Helper function to compute a single split for a given set of W_m and
+    eigenvalue lambda_i, when the total number of terms in the split is odd.
+    """
+    Wm_array = np.array([W.value for W in Wm])
+    ND = _single_split_odd_values(Wm_array, li, r)
+    N = Wm[0].layout.MultiVector(ND[0, :])
+    D = Wm[0].layout.MultiVector(ND[1, :])
+    return N*D.leftLaInv()
+
+
+def _bivector_split(Wm, return_all=True):
+    """Internal helper function to perform the decomposition, given a set of Wm.
+    Parameters
+    ----------
+    return_all : bool, optional
+        If `True`, returns all the :math:`b_i`.
+        If `False`, return all :math:`b_i` except for the one with the smallest magnitude.
+    """
+    # The effective value of k is determined by the largest non-zero W.
+    # remove the highest grade zeros to prevent meaningless lambda_i = 0 values.
+    for W in Wm[::-1]:
+        if np.linalg.norm(W.value) > _eps:
+            break
+        else:
+            Wm = Wm[:-1]
+
+    k = (len(Wm) - 1)
+    r = k // 2
+    Wm_sq = np.array([(W ** 2).value[0] * (-1) ** (k - m) for m, W in enumerate(Wm)])
+    ls = np.roots(Wm_sq)
+
+    Bs = []
+    # Sort to have the value closest to zero last.
+    ls_sorted = sorted(ls, key=lambda li: -li)
+    # Exclude the smallest value if asked.
+    for li in (ls_sorted if return_all else ls_sorted[:-1]):
+        if k % 2 == 0:
+            Bs.append(single_split_even(Wm, li, r))
+        else:
+            Bs.append(single_split_odd(Wm, li, r))
+    return (Bs, ls_sorted)
+
+
+def bivector_split(B, k=None, roots=False):
+    """Bivector split of the bivector B based on the method of :cite:`roelfs2021thesis`, chapter 6.
+
+
+    Parameters
+    ----------
+    roots : bool, optional
+        If `True`, return the values of the :math:`\lambda_i` in addition to the :math:`b_i`.
+        If `False`, return only the :math:`b_i`.
+    """
+    dim = B.layout.dims
+    if k is None:
+        k = dim // 2
+
+    Wm = [(B**m)(2*m) / math.factorial(m) for m in range(0, k + 1)]
+    Bs, ls = _bivector_split(Wm, return_all=False)
+    Bs = Bs + [B - sum(Bs)]
+    return (Bs, ls) if roots else Bs
+
+
+def rotor_split(R, k=None, roots=False):
+    dim = R.layout.dims
+    if k is None:
+        k = dim // 2
+
+    Wm = [R(2 * m) for m in range(0, k + 1)]
+    Ts, ls = _bivector_split(Wm, return_all=False)
+
+    Rs = [(1 + ti) for ti in Ts]
+    Rs = [Ri.normal() if np.isreal((Ri*~Ri).value[0]) else Ri / np.sqrt((Ri*~Ri).value[0]) for Ri in Rs]
+    P = reduce(lambda tot, x: tot*x, Rs, 1.0 + 0.0*R)
+    Rs = Rs + [R*P.leftLaInv()]
+    return (Rs, ls) if roots else Rs
+
+
+def exp(B):
+    Bs, ls = bivector_split(B, roots=True)
+    R = B.layout.scalar
+    for Bi, li in zip(Bs, ls):
+        if np.isreal(li) and li < 0:
+            beta_i = np.sqrt(-li)
+            R *= np.cos(beta_i) + (np.sin(beta_i) / beta_i) * Bi
+        elif np.isreal(li) and np.abs(li) < _eps:
+            R *= 1 + Bi
+        else:
+            beta_i = np.sqrt(li)
+            R *= np.cosh(beta_i) + (np.sinh(beta_i) / beta_i) * Bi
+    return R
+
+
+def log(R):
+    Rs, ls = rotor_split(R, roots=True)
+    logR = R.layout.MultiVector()
+    for Ri, li in zip(Rs, ls):
+        if li < 0:
+            norm = np.sqrt(- (Ri(2) ** 2).value[0])
+            logR += np.arccos(Ri.value[0]) * Ri(2) / norm
+        elif np.abs(li) < _eps:
+            logR += Ri(2)
+        else:
+            norm = np.sqrt((Ri(2)**2).value[0])
+            logR += np.arccosh(Ri.value[0]) * Ri(2) / norm
+    return logR

clifford/test/__init__.py

@@ -2,13 +2,20 @@
 import pytest
 import numpy as np
 
+from clifford._numba_utils import DISABLE_JIT
+
 
 @pytest.fixture
 def rng():
     default_test_seed = 1  # the default seed to start pseudo-random tests
     return np.random.default_rng(default_test_seed)
 
 
+too_slow_without_jit = pytest.mark.skipif(
+    DISABLE_JIT, reason="test is too slow without JIT"
+)
+
+
 def run_all_tests(*args):
     """ Invoke pytest, forwarding options to pytest.main """
     pytest.main([os.path.dirname(__file__)] + list(args))

clifford/test/test_invariant_decomposition.py

@@ -0,0 +1,176 @@
+from functools import reduce, lru_cache
+import itertools
+import operator
+
+import pytest
+import numpy as np
+
+from clifford import Layout, BasisVectorIds
+from clifford.invariant_decomposition import bivector_split, rotor_split, exp, log
+
+import clifford as cf
+from . import rng  # noqa: F401
+from . import too_slow_without_jit
+
+
+# Test some known splits in various algebras. The target bivector is `B`,
+# the expected split is `Bs`, the expected eigenvalues li = Bi**2 are
+# given in `ls`. Lastly, the expected logarithm is `logR`.
+# We use `lru_cache` here so that these layouts can be reused between tests,
+# but are not constructed at all if not needed.
+@lru_cache
+def sta_split():
+    alg = Layout([1, 1, 1, -1], ids=BasisVectorIds(['x', 'y', 'z', 't']))
+    ex, ey, ez, et = alg.basis_vectors_lst
+    return {'B': 2 * ex * ey + 4 * ez * et,
+            'Bs': [4 * ez * et, 2 * ex * ey],
+            'ls': [16.0, -4.0],
+            'logR': 2 * ex * ey + 4 * ez * et}
+
+@lru_cache
+def pga3d_split():
+    alg = Layout([1, 1, 1, 0], ids=BasisVectorIds(['x', 'y', 'z', 'w']))
+    ex, ey, ez, ew = alg.basis_vectors_lst
+    return {'B': 2 * ex * ey + 4 * ez * ew,
+            'Bs': [4 * ez * ew, 2 * ex * ey],
+            'ls': [0.0, -4.0],
+            'logR': 2 * ex * ey + 4 * ez * ew}
+
+@lru_cache
+def r22_split():
+    alg = Layout([1, 1, -1, -1])
+    e1, e2, e3, e4 = alg.basis_vectors_lst
+    return {'B': 0.5 * (e1*e2 + e1*e4 - e2*e3 - e3*e4),
+            'Bs': [0.25 * ((1-1j)*e1*e2 + (1+1j)*e1*e4 + (-1-1j)*e2*e3 + (-1+1j)*e3*e4),
+                   0.25 * ((1+1j)*e1*e2 + (1-1j)*e1*e4 + (-1+1j)*e2*e3 + (-1-1j)*e3*e4)],
+            'ls': [0.5j, -0.5j],
+            'logR': 0.5 * (e1*e2 + e1*e4 - e2*e3 - e3*e4)}
+
+@lru_cache
+def r6_split():
+    alg = Layout([1, 1, 1, 1, 1, 1])
+    e1, e2, e3, e4, e5, e6 = alg.basis_vectors_lst
+    return {'B': 2*e1*e2 + 5*e3*e4 + 7*e5*e6,
+            'Bs': [2*e1*e2, 5*e3*e4, 7*e5*e6],
+            'ls': [-4.0, -25.0, -49.0],
+            # The log is by no means unique, and this example illustrates it.
+            # With the conventions of this implementation this is the answer, which is a
+            # total of 4 pi away from the input B.
+            'logR': (2 - np.pi)*e1*e2 + (5 - np.pi)*e3*e4 + (7 - 2*np.pi)*e5*e6}
+
+@lru_cache
+def r4_split():
+    alg = Layout([1, 1, 1, 1])
+    e1, e2, e3, e4 = alg.basis_vectors_lst
+    delta = 1
+    return {'B': 2*e1*e2 + (2+delta)*e3*e4,
+            'Bs': [(2+delta)*e3*e4, 2*e1*e2],
+            'ls': [-(2+delta)**2, -4.0],
+            'logR': (2-np.pi)*e1*e2 + (2+delta-np.pi)*e3*e4}
+
+
+class TestInvariantDecomposition:
+    """ Test the invariant decomposition of bivectors and rotors, and the resulting exp and log functions.
+    """
+    @pytest.mark.parametrize('known_split', [
+        pytest.param(sta_split, id='sta'),
+        pytest.param(pga3d_split, id='pga'),
+        pytest.param(r22_split, id='r22'),
+        pytest.param(r6_split, id='r6'),
+        pytest.param(r4_split, id='r4', marks=[
+            pytest.mark.xfail(reason='unknown')]),
+    ])
+    def test_known_splits(self, known_split):
+        known_split = known_split()
+
+        B = known_split['B']
+        Bs, ls = bivector_split(B, roots=True)
+        # Test the bivector split
+        for calculated, known in zip(Bs, known_split['Bs']):
+            np.testing.assert_allclose(calculated.value, known.value, atol=1E-5, rtol=1E-5)
+        np.testing.assert_allclose(ls, known_split['ls'])
+
+        # Test the exp function agrees with the taylor expansion used by `.exp()`
+        R = exp(B)
+        Rraw = reduce(operator.mul, [Bi.exp() for Bi in Bs])
+        np.testing.assert_allclose(Rraw.value, R.value, atol=1E-5, rtol=1E-5)
+
+        for Bi in Bs:
+            # Test if the simple bivectors are exponentiated correctly.
+            np.testing.assert_allclose(Bi.exp().value, exp(Bi).value, atol=1E-5, rtol=1E-5)
+
+    @pytest.mark.parametrize('known_split', [
+        pytest.param(sta_split, id='sta'),
+        pytest.param(pga3d_split, id='pga'),
+        pytest.param(r22_split, id='r22'),
+        pytest.param(r6_split, id='r6'),
+        pytest.param(r4_split, id='r4'),
+    ])
+    def test_known_rotor_splits(self, known_split):
+        known_split = known_split()
+        B = known_split['B']
+
+        R = exp(B)
+
+        # Split R into simple rotors.
+        Rs = rotor_split(R)
+        # Test simpleness of the Ri
+        for Ri in Rs:
+            np.testing.assert_allclose(Ri.value, Ri(0, 2).value, atol=1E-5, rtol=1E-5)
+        # Test commutativity
+        for Ri, Rj in itertools.combinations(Rs, 2):
+            np.testing.assert_allclose((Ri*Rj).value, (Rj*Ri).value, atol=1E-5, rtol=1E-5)
+
+        # Reconstruct R from the rotor_split.
+        Rre = reduce(operator.mul, Rs, B.layout.scalar)
+        np.testing.assert_allclose(R.value, Rre.value, atol=1E-5, rtol=1E-5)
+
+        logR = log(R)
+        np.testing.assert_allclose(logR.value, known_split['logR'].value, atol=1E-5, rtol=1E-5)
+
+    @pytest.mark.parametrize('r', range(2))
+    @pytest.mark.parametrize('p, q', [
+        pytest.param(p, total_dims - p,
+            marks=[pytest.mark.veryslow, too_slow_without_jit] if total_dims >= 6 else [])
+        for total_dims in [0, 1, 2, 3, 4, 5, 6, 7, 8]
+        for p in range(total_dims + 1)
+    ])
+    def test_unknown_splits(self, p, q, r, rng):  # noqa: F811
+        Ntests = 10
+        layout, blades = cf.Cl(p, q, r)
+        for i in range(Ntests):
+            B = layout.randomMV(rng=rng, grades=[2])
+            Bs, ls = bivector_split(B, roots=True)
+            for Bi, li in zip(Bs, ls):
+                # To be simple, you must square to a scalar.
+                Bisq = Bi**2
+                np.testing.assert_allclose(Bisq.value, Bisq(0).value,
+                                           rtol=1E-5, atol=1E-5)
+                np.testing.assert_allclose(Bisq.value[0], li,
+                                           rtol=1E-5, atol=1E-5)
+
+            # Assert that the bivectors sum to the original
+            np.testing.assert_allclose(sum(Bs).value, B.value,
+                                       rtol=1E-5, atol=1E-5)
+
+            # Assert that the bivectors are commutative
+            for x, y in itertools.combinations(Bs, 2):
+                np.testing.assert_allclose((x*y).value, (y*x).value,
+                                           rtol=1E-5, atol=1E-5)
+
+            R = exp(B)
+
+            # Make the absolute tolerance of the rotor tests dependent on total dimension
+            default_atol = 1E-5
+            if p + q + r > 7:
+                default_atol = 1E-2
+
+            # Assert that the rotor split multiplies into the original.
+            Rs, ls = rotor_split(R, roots=True)
+            Rre = reduce(operator.mul, Rs, 1)
+            np.testing.assert_allclose(R.value, Rre.value,
+                                       rtol=1E-5, atol=default_atol)
+            # Assert that the exp(log(R)) is the original rotor R.
+            logR = log(R)
+            np.testing.assert_allclose(R.value, exp(logR).value,
+                                       rtol=1E-5, atol=default_atol)