implement left and right complements for ⋅,∨ products

[chakravala]

this is a work-in-progress for implementing the fully-generalized interior and regressive product #2

some help is required to fix the interaction between the metric and the interior product operations

Complement

The complement operations were implemented as two methods

complementleft(x)
complementright(x)

with the extra definition const ⋆ = complementright for the Hodge star linear map

determining the basis value for the complement operation is a fairly simple bitwise operation

complement(N::Int,B::UInt) = (~B)&(one(UInt)<<N-1)

while the parity of the left and right operations are determined by the indices and metric

parityright(V::Int,B,G,N) = isodd(V+B+Int((G+1)*G/2))
parityleft(V::Int,B,G,N) = (isodd(G) && iseven(N)) ⊻ parityright(V,B,G,N)
function parityright(V::VectorSpace,B,G=count_ones(B))
    b=indices(B)
    parityright(sum(V[:][b]),sum(b),G,ndims(V)))
 end
function parityleft(V::VectorSpace,B,G=count_ones(B))
    b=indices(B)
    parityleft(sum(V[:][b]),sum(b),G,ndims(V)))
 end

which I believe is the definition based on Browne's Grassmann algebra, where the parity is determined by the sum of the indices plus triangular number of the grade, in addition to the metric components

Regressive and Interior products

This leads to the definition of the Grassmann regressive product in terms of the complements

∨(a,b) = complementright(complementleft(a)∧complementleft(b))

which is based on the left and right complements as discussed in AntiWedge or Regressive Product

dot(a,b) = complementright(complementleft(a)∧b)

and then the interior or inner product is defined similarly, with a modification

Issues and results

for the positive definite VectorSpace metric, this approach does seem to work successfully

julia> using Grassmann; basis"4"
(⟨++++⟩, v, v₁, v₂, v₃, v₄, v₁₂, v₁₃, v₁₄, v₂₃, v₂₄, v₃₄, v₁₂₃, v₁₂₄, v₁₃₄, v₂₃₄, v₁₂₃₄)

julia> v1+v2+2v4+v3
1v₁ + 1v₂ + 1v₃ + 2v₄

julia> ans⋅ans
7

julia> v1*v1,v1⋅v1,v1∧v1
(v, 1v, 0v)

julia> v2*v2,v2⋅v2,v2∧v2
(v, v, 0v)

so all of the products work as expected when dealing with Euclidean + only spaces

however, when the Minkowski spacetime metric is used, some of the issues with the formulation show

julia> basis"-+++"
(⟨-+++⟩, v, v₁, v₂, v₃, v₄, v₁₂, v₁₃, v₁₄, v₂₃, v₂₄, v₃₄, v₁₂₃, v₁₂₄, v₁₃₄, v₂₃₄, v₁₂₃₄)

julia> v1*v1,v1⋅v1,v1∧v1
(-1v, 1v, 0v)

julia> v2*v2,v2⋅v2,v2∧v2
(v, -1v, 0v)

as you can see, in this case the parity of the inner product is reversed from what it should be

this is a result due to the definition of the interior product in terms of the right and left complement

Request for help

please let me know if any of the geometric algebra enthusiasts out there (such as @waldyrious, @enkimute, @arsenovic, @hugohadfield, @joanlasenby could help clarify or figure out what the solution is to making this implementation work correctly in terms of the left and right complements

any comments, corrections, suggestions, or references are welcome

[codecov-io]

Codecov Report

Merging #3 into master will increase coverage by 1%.
The diff coverage is 26.53%.

@@           Coverage Diff           @@
##           master       #3   +/-   ##
=======================================
+ Coverage   24.92%   25.92%   +1%     
=======================================
  Files           7        7           
  Lines         967     1022   +55     
=======================================
+ Hits          241      265   +24     
- Misses        726      757   +31

Impacted Files	Coverage Δ
src/multivectors.jl	22.85% <ø> (ø)	⬆️
src/forms.jl	16.66% <0%> (ø)	⬆️
src/utilities.jl	43.42% <0%> (-1.79%)	⬇️
src/symbolic.jl	16.66% <13.33%> (-6.87%)	⬇️
src/algebra.jl	16.85% <28.9%> (+3.83%)	⬆️
src/generators.jl	81.81% <0%> (+0.23%)	⬆️

---

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 8ad0169...5bb7fde. Read the comment docs.

[enkimute]

I'd have to look into Lengyels definition more closely (although I would like to note that picking overbar and underbar for 'left' and 'right' complement are not the most intuitive choices ;) )

Dorst has some details on various inner products in a separate article here :

https://www.researchgate.net/publication/2842332_The_Inner_Products_of_Geometric_Algebra

Gunn also has an alternate definition of duality, and thus the regressive product, in his work on PGA using a degenerate metric. (He uses Poincare duality).

https://www.researchgate.net/publication/265672134_Geometry_Kinematics_and_Rigid_Body_Mechanics_in_Cayley-Klein_Geometries

Ganja uses Poincare duality when dealing with the degenerate metric, and the left contraction inner product from Dorst. (although a symmetric inner product is also available)

when I have some time I'll check the compatibility between Gunn's and Lengyels definition.

[chakravala]

Thanks for your comments @enkimute, after looking at Gunn's work on Poincare duality, I did some experimentation. So if I understand correctly, the point is that you are taking the dual of the complement.

In Grassmann we already have the dual involution ' available, and in this pull request you will find both of the complementleft and complementright operations, allowing you to experiment with their compositions in various different signature algebras.

julia> basis"++++"
(⟨++++⟩, v, v₁, v₂, v₃, v₄, v₁₂, v₁₃, v₁₄, v₂₃, v₂₄, v₃₄, v₁₂₃, v₁₂₄, v₁₃₄, v₂₃₄, v₁₂₃₄)

julia> Grassmann.complementleft(v1)' ∧ v1'
1w¹²³⁴

julia> (Grassmann.complementleft(v1) ∧ v1)'
1w¹²³⁴

As you can see, the dual ' operation is an outermorphism, which distributes with complements.

Let's define a norm using the left and right complement and also the dual operations

julia> norm(x) = Grassmann.complementright((Grassmann.complementleft(x)∧x)')'
norm (generic function with 1 method)

julia> norm(v1),norm(v2)
(1v, v)

But now let's change the basis to the Minkowski example from above

julia> basis"-+++"
(⟨-+++⟩, v, v₁, v₂, v₃, v₄, v₁₂, v₁₃, v₁₄, v₂₃, v₂₄, v₃₄, v₁₂₃, v₁₂₄, v₁₃₄, v₂₃₄, v₁₂₃₄)

julia> norm(v1),norm(v2)
(1v, -1v)

We still have the same exact issue as before, apparently.

Could you try to express your definitions of inner product using these operations ' and ⋆?

Also, have you used the definition of Hodge duality before?

[chakravala]

In the end, I resolved the issue by sticking to the method presented by Browne in Grassmann algebra,

L = grade(a) + grade(b)
(-1)^(L*(L-ndims(V))) * complementright(complementright(a)∧complementright(b)))

his formula for the regressive product is only in terms of complementright and a grade factor

The expression of a regressive product in terms of an exterior product and a double application of the complement operation poses an interesting conundrum. On the left hand side, there is a product which is defined in a completely non-metric way, and is thus independent of any metric imposed on the space. On the right hand side we have a double application of the complement operation, each of which requries a metric. This leads us to the notion that the complement operation applied twice in this way is independent of any metric used to define it. A second application of the complement operation effectively 'cancels out' the metric elements introduced during the first application.

sticking to this method, I found I could make a very efficient interior product implementation, where I essentially skip over the required intermediate allocations and determine the parity in one step:

function interior_calc(N,M,S,A,B)
    γ = complement(N,B)
    p,C,t = regressive_calc(N,M,S,A,γ)
    return t ? p⊻parityright(S,B) : p, C, t
end

which in turn depends on the parity calculation obtained from the regressive product

function regressive_calc(N,M,S,A,B)
    α,β = complement(N,A),complement(N,B)
    if (count_ones(α&β)==0) && !(M ∈ (1,3,5,7,9,11) && isodd(α) && isodd(β))
        C = α ⊻ β
        L = count_ones(A)+count_ones(B)
        pa,pb,pc = parityright(S,A),parityright(S,B),parityright(S,C)
        return isodd(L*(L-N))⊻pa⊻pb⊻parity(N,S,α,β)⊻pc, complement(N,C), true
    else
        return false, zero(Bits), false
    end
end

however, since this is a rather cumbersome calculation requiring many micoseconds to compute a single boolean value for the parity.. I have opted to cache the value efficiently into the compiled code

julia> using Grassmann; basis"-+++"
(⟨-+++⟩, v, v₁, v₂, v₃, v₄, v₁₂, v₁₃, v₁₄, v₂₃, v₂₄, v₃₄, v₁₂₃, v₁₂₄, v₁₃₄, v₂₃₄, v₁₂₃₄)

julia> @btime h = 2v1+v3
  37.794 ns (3 allocations: 80 bytes)
2v₁ + 0v₂ + 1v₃ + 0v₄

julia> @btime h⋅h
  105.948 ns (2 allocations: 160 bytes)
-3

and as a result I am getting good sub-microsecond results for high-performnace computing

it would still be interesting to find out more about Lengyel's approach, but this one seems solid now

by the way, @enkimute your GanjaCheatSheets helped me a lot from your blog :)

[enkimute]

That is a very interesting quote - I will read up on that, will have to check the way he defined the complement and figure out if this may provide a method that works for degenerate metrics in the same unified fashion..

Cool update ! and happy the cheat-sheets came in handy 👍