Bivector and rotor splits for arbitrary GAs

[tBuLi]

This PR aims to add the bivector split of chapter 6 of my thesis to clifford. When finished, this PR will add the decomposition of arbitrary bivectors into mutually committing orthogonal simple bivectors, and the decomposition of rotors into mutually commuting orthogonal simple rotors. Additionally, this enables us to implement closed form implementations for the exp/log function of simple bivectors/rotors respectively.

ToDo

• Invariant decomposition of bivectors into simple bivectors
• Closed form exponential of bivectors
• Invariant decomposition of rotors into simple rotors
• Logarithm of rotors
• Exception handling: the algorithm will fail when the eigenvalues are degenerate. We need to discuss how to handle these cases, although this cannot occur in any the popular small algebras of 3DCGA, 3DPGA, STA, and it might be okay to leave it up to the user.
• Optimize tests: a lot of copy-pasting went into the current tests, pytest probably has a solution for this.
• Documentation.

Curious to hear your opinions on the work done so far.

[hugohadfield]

This looks great, hopefully the inclusion of the general decomposition of bivectors means we can remove the various algebra specific bivector decompositions in other places in the code :) I'll have a dig for them later so we can do some spring cleaning

[tBuLi]

@hugohadfield that sounds great! Cleaning up the code base is what GA is for ;).

I have made one potentially contentious change to _layout.py: I've removed the checking if a MV is invertible by removing the following lines:

if abs(np.linalg.det(intermed)) < _settings._eps:
    raise ValueError("multivector has no left-inverse")

and instead will let numpy decide. Maybe there is a good reason this line was added that I'm not aware of, but in my experiments with bivector splits I've regularly walked into this condition whereas I knew the requested inversion was possible. Also decreasing _eps further doesn't always work, whereas removing this check all together works like a charm. So I'm curious to hear what the reason was for adding this check. (If interested in a MWE, the r6 rotor_split in the tests is uninvertible according to this check, but it isn't.)

[lgtm-com]

This pull request introduces 1 alert when merging f4e56da into 31651d0 - view on LGTM.com

new alerts:

• 1 for Unused import

[tBuLi]

@hugohadfield, when it comes to cleaning up the code base, should that be part of this PR or is it better to leave this one focused on adding the general solution and do the clean-up separately?

[hugohadfield]

Good question, I think the file that is likely to be affected most by this merge is:
https://github.com/pygae/clifford/blob/master/clifford/tools/g3c/rotor_parameterisation.py
and the bivector split is specifically here:

clifford/clifford/tools/g3c/rotor_parameterisation.py

Line 52 in 31651d0

def decompose_bivector(F):

There is also some functionality from @arsenovic here:

clifford/clifford/tools/__init__.py

Line 462 in 31651d0

def log_rotor(V):

And the exponential is chosen currently as the taylor series by default here:

clifford/clifford/_multivector.py

Line 101 in 31651d0

def exp(self) -> 'MultiVector':

[hugohadfield]

It looks to me like the linear algebra method for inverting multivectors is not happy with the complex numbers and we might have to merge this: #373 first in order to get an inverse that is happy

[hugohadfield]

Looking at the talking points on here, specifically:

Exception handling: the algorithm will fail when the eigenvalues are degenerate. We need to discuss how to handle these cases, although this cannot occur in any the popular small algebras of 3DCGA, 3DPGA, STA, and it might be okay to leave it up to the user.

How does this failure manifest itself? Like, does the algorithm just give the wrong result or do we get a divide by zero or something like that?

Can we construct a test case that doesn't work? If so we should include some test cases to ensure that we handle these failures in a reasonable way :)

I'm also tempted to build in something to test that exp and log definitely round trip for the principal logarithm of rotors in loads of different spaces

[tBuLi]

Due to numerical errors, there are two possibilities. In theory the degenerate cases lead to a "no inverse exists" error, although in practice I've had instances where they are different enough for the code to just continue anyway, especially if the eigenvalues pick up a small imaginary part like lets say, lambda_1 = 1.0 + 1e-8j and lambda_2 = 1.0 - 1e-8j. But despite not giving an error, the result will not be a 2-blade. The best would be to compute all abs((Bi**2)(4)) and verify that they are close to zero, but this is a bit expensive. Perhaps a warning when the eigenvalues start to become close will suffice?

I'll think about and implement a test case later.

I'm also tempted to build in something to test that exp and log definitely round trip for the principal logarithm of rotors in loads of different spaces

I agree with that, I think it would be a great idea to make some tests to verify that

R = exp(B)
assert R == exp(log(R))

if I understood you correctly.

[tBuLi]

So if I create the known_split

def r4_split():
    alg = Layout([1, 1, 1, 1])
    e1, e2, e3, e4 = alg.basis_vectors_lst
    return {'B': 2*e1*e2 + 2*e3*e4,
            'Bs': [2*e3*e4, 2*e1*e2],
            'ls': [-4.0, -4.0],
            'logR': 2*e1*e2 + 2*e3*e4}

then test_known_splits errors with the traceback

clifford/test/test_invariant_decomposition.py:72 (TestInvariantDecomposition.test_known_splits[known_split0])
self = <clifford.test.test_invariant_decomposition.TestInvariantDecomposition object at 0x16e35fb20>
known_split = {'B': (2^e12) + (2^e34), 'Bs': [(2^e34), (2^e12)], 'logR': (2^e12) + (2^e34), 'ls': [-4.0, -4.0]}

    @pytest.mark.parametrize('known_split', [r4_split()])
    def test_known_splits(self, known_split):
        B = known_split['B']
>       Bs, ls = bivector_split(B, roots=True)

test/test_invariant_decomposition.py:76: 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 
invariant_decomposition.py:101: in bivector_split
    Bs, ls = _bivector_split(Wm, return_all=False)
invariant_decomposition.py:82: in _bivector_split
    Bs.append(single_split(Wm, li))
invariant_decomposition.py:53: in single_split
    return N*D.inv()
_multivector.py:795: in inv
    return self._pick_inv(fallback=True)
_multivector.py:767: in _pick_inv
    return self.hitzer_inverse()
_multivector.py:741: in hitzer_inverse
    return self.layout._hitzer_inverse(self)
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 

>   raise ValueError('Multivector has no inverse')
E   ValueError: Multivector has no inverse

_layout.py:634: ValueError

However, if I split B = 2*e1*e2 + (2+1e-7)*e3*e4 instead, I get the split
[-(0.06108^e12) + (2.06108^e34), (2.06108^e12) - (0.06108^e34)],
which both square to -4.2518 - (0.2518^e1234) and thus aren't simple. So numerically speaking, the problems with the inverse being not defined already start a bit earlier. Perhaps it will be enough to give a warning when the eigenvalues start to become degenerate?

[tBuLi]

I've now added tests for exp(log(R)) == R in higher dimensions as well. This slows down the tests significantly but was well worth it because it smoked out some bugs in the implementation. Most importantly, it made me realise that you have to sort from boost to rotations, i.e. from positive to negative eigenvalues, before performing the split, to prevent boosts from ever getting a negative scalar part. However, this means that zero could be somewhere in the middle, and thus I've had to implement the even and odd splits properly like in my thesis, whereas before I used a single form which is valid for non-zero eigenvalues and then computed the remaining split with B - sum(Bs).

In the spaces 3-5-1, 4-4-1, 5-3-0, and 5-3-1 I now get an error that the log doesn't work up to the desired precision, but that appears to be purely a numerical precision problem, its still correct-ish since np.linalg.norm(R.value - exp(log(R)).value) ~= 1e-7. Perhaps we should make the demanded precision in the test dimension dependent?

@hugohadfield, is it possible to only iterate spaces with p >= q? That saves some computation time and the even subalgebras of p,q and q,p spaces appear to be isomorphic anyway.

[hugohadfield]

I think adding the dimension dependent accuracy check is a good idea, it would also be a good idea to profile the tests a little and see where the code is spending most of its time, unfortunately I suspect it is in algebra creation. We should also think about maybe numba JITing the various functions in this PR, that way users can include them inside jitted functions for speed etc.

[codecov]

Codecov Report

Merging #398 (e206c24) into master (e63f856) will decrease coverage by 38.77%.
The diff coverage is 37.36%.

@@             Coverage Diff             @@
##           master     #398       +/-   ##
===========================================
- Coverage   71.37%   32.59%   -38.78%     
===========================================
  Files          73       75        +2     
  Lines        9459     9635      +176     
  Branches     1211     1244       +33     
===========================================
- Hits         6751     3141     -3610     
- Misses       2538     6381     +3843     
+ Partials      170      113       -57     

Impacted Files	Coverage Δ
clifford/_layout.py	56.13% <ø> (-40.83%)	⬇️
clifford/test/test_invariant_decomposition.py	31.81% <31.81%> (ø)
clifford/invariant_decomposition.py	40.65% <40.65%> (ø)
clifford/test/__init__.py	72.72% <100.00%> (-16.17%)	⬇️
clifford/test/test_multivector_inverse.py	42.30% <100.00%> (-57.70%)	⬇️
clifford/test/test_dpga.py	9.52% <0.00%> (-90.48%)	⬇️
clifford/test/test_degenerate.py	10.25% <0.00%> (-89.75%)	⬇️
clifford/test/test_dg3c.py	15.92% <0.00%> (-84.08%)	⬇️
clifford/_parser.py	16.39% <0.00%> (-83.61%)	⬇️
clifford/test/test_clifford.py	15.67% <0.00%> (-83.40%)	⬇️
... and 51 more

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update e63f856...e206c24. Read the comment docs.

[eric-wieser]

@hugohadfield, I assume this jitted code is your doing. Is there a reason you define the values function vs just using the multivector jit support?

[eric-wieser]

Maybe #412 will resolve this.

[arsenovic]

whats the status on this? any chance we could get it pushed in some form so i can use it?
i mainly want to use in 4d.

[arsenovic]

just checking in again. this is a really great feature!