diff --git a/clifford/_layout.py b/clifford/_layout.py index 682c5478..3e6a7dba 100644 --- a/clifford/_layout.py +++ b/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 diff --git a/clifford/invariant_decomposition.py b/clifford/invariant_decomposition.py new file mode 100644 index 00000000..ac6eed82 --- /dev/null +++ b/clifford/invariant_decomposition.py @@ -0,0 +1,185 @@ +""" +.. currentmodule:: clifford.invariant_decomposition + +===================================================== +invariant_decomposition (:mod:`clifford.invariant_decomposition`) +===================================================== + +.. versionadded:: 1.5.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:: + + >>> from clifford.g4 import * + >>> B = 1*e12 + 2*e34 + >>> bivector_split(B) + [((1+0j)^e12), ((2+0j)^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): + r"""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 diff --git a/clifford/test/__init__.py b/clifford/test/__init__.py index 18ffbb02..393c1927 100644 --- a/clifford/test/__init__.py +++ b/clifford/test/__init__.py @@ -2,6 +2,8 @@ import pytest import numpy as np +from clifford._numba_utils import DISABLE_JIT + @pytest.fixture def rng(): @@ -9,6 +11,11 @@ def rng(): 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)) diff --git a/clifford/test/test_invariant_decomposition.py b/clifford/test/test_invariant_decomposition.py new file mode 100644 index 00000000..a3904b8c --- /dev/null +++ b/clifford/test/test_invariant_decomposition.py @@ -0,0 +1,180 @@ +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(maxsize=None) +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(maxsize=None) +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(maxsize=None) +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(maxsize=None) +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(maxsize=None) +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) diff --git a/clifford/test/test_multivector_inverse.py b/clifford/test/test_multivector_inverse.py index 6f1a8584..38a829ab 100644 --- a/clifford/test/test_multivector_inverse.py +++ b/clifford/test/test_multivector_inverse.py @@ -3,13 +3,7 @@ import clifford as cf from . import rng # noqa: F401 - -from clifford._numba_utils import DISABLE_JIT - - -too_slow_without_jit = pytest.mark.skipif( - DISABLE_JIT, reason="test is too slow without JIT" -) +from . import too_slow_without_jit class TestClosedForm: diff --git a/docs/api/index.rst b/docs/api/index.rst index d6509af4..c06eeb57 100644 --- a/docs/api/index.rst +++ b/docs/api/index.rst @@ -11,4 +11,5 @@ API operator transformations taylor_expansions + invariant_decomposition numba diff --git a/docs/api/invariant_decomposition.rst b/docs/api/invariant_decomposition.rst new file mode 100644 index 00000000..d14c1cf5 --- /dev/null +++ b/docs/api/invariant_decomposition.rst @@ -0,0 +1 @@ +.. automodule:: clifford.invariant_decomposition diff --git a/docs/refs.bib b/docs/refs.bib index 0eedd75f..86e63054 100644 --- a/docs/refs.bib +++ b/docs/refs.bib @@ -135,3 +135,14 @@ @book{hestenes-space-time year = {2015}, keywords = {Mathematical physics, Space and time} } + +@misc{roelfs2021thesis, + abstract = {The best theory currently available to describe the formation of hadronic particles like protons and neutrons is Quantum Chromodynamics (QCD). It predicts the existence of quarks and gluons, although neither of these particles have been observed directly. The theory predicts this is due to a phenomenon called confinement, which forbids quarks and gluons to travel unaccompanied. However, they are believed to be real physical entities, as the predictions of this model are in excellent agreement with a wide range of experiments. But perhaps one of the most exciting predictions of this theory has not yet been observed: glueballs. Glueballs are neutral particles consisting purely of gluons. Some glueball candidates have been identified in particle accelerators, but physicist cannot be sure yet whether these are really the glueballs predicted by QCD. This is due to the great difficulty of separating them from other particles which have similar masses. Therefore more accurate predictions of the glueball mass and other spectral properties are needed. Another type of particle which is unobservable, is the Faddeev-Popov ghost. These particles are introduced to gauge fix the theory: a mathematical necessity to allow the calculation of scattering amplitudes. As the Faddeev-Popov ghosts are obtained by rewriting the mathematical machinery of gauge fixing in the language of particle physics, they are not just unobservable; they are unphysical. The primary goal of this thesis was to use lattice QCD to calculate the spectral properties of gluons, glueballs, and Faddeev-Popov ghosts. This was done by finding the Källén-Lehmann spectral representation from the propagators of these particles. Additionally, a method for calculating the Faddeev-Popov ghosts in linear covariant gauges on the lattice was proposed. As QCD is a local gauge theory over the non-abelian group SU(3), a new set of mathematical tools to study SU(3) was also developed. These tools find their origin in geometric algebra, but have also been translated to the more traditional matrix representation of SU(3). The results are presented using both the geometric algebra and the matrix representation. Lastly, the exploration of geometric algebra lead to the discovery of a novel technique for the numerical computation of derivatives of holomorphic functions, using complex quaternions. This chapter stands somewhat apart from the others, although the complex quaternions are isomorphic to SU(2) x U(1), which is the gauge group of the electroweak force. Therefore this paper does briefly describe a different way of understanding SU(2) x U(1) compared to the usual matrix formalism, although it then quickly digresses into the derivation of holomorphic functions.}, + year = {2021}, + title = {Spectroscopic and Geometric Algebra Methods for Lattice Gauge Theory}, + language = {eng}, + author = {Roelfs, Martin}, + url = {https://lirias.kuleuven.be/retrieve/620480}, + isbn = {978-94-6416-639-2} +} +