Based on my brief look at this, I don’t believe what you’re trying to do is going to work.
Let’s just start with only thinking about how quaternions would be represented as a sub-matrix within your multivector matrix. To start off, you’d need to have a way to translate from multivectors to quaternions. Let’s do that first. In Grassmann.jl this is with:
julia> using Grassmann, AbstractTensors; basis"3"
(⟨×××⟩, v, v₁, v₂, v₃, v₁₂, v₁₃, v₂₃, v₁₂₃)
julia> 𝕚,𝕛,𝕜
(-1v₂₃, 1v₁₃, -1v₁₂)
julia> hyperplanes(ℝ3)
3-element Vector{Single{⟨×××⟩, 2, B, Int64} where B}:
-1v₂₃
1v₁₃
-1v₁₂
julia> 𝕚*𝕛 == 𝕜, 𝕛*𝕜 == 𝕚, 𝕜*𝕚 == 𝕛, 𝕚*𝕛*𝕜 == -1
(true, true, true, true)
The generalization of this is to all grades is to use the transformation -I\omega:
julia> basis = Ref(-I).*Values(v,v1,v2,v3,v12,v13,v23,v123)
(-1v₁₂₃, -1v₂₃, 1v₁₃, -1v₁₂, 1v₃, -1v₂, 1v₁, 1v)
Now we have an actual multivector basis containing the quaternion basis, unlike what you’ve been doing in your posts (although the ordering is different, it doesn’t matter since matrix algebra is invariant modulo swapping the order of the basis). Now let’s define the multivector “involutions” required for the matrix definition using similar concepts as written in the paper you referenced. First, we need a few constructors:
fun(a,b) = a==b ? b : ~b
create(x,b) = Chain(fun.(Ref(x),b))
(Single{V,G,B,T} where {G,B})(::Zero{V}) where {V,T} = zero(T)*One{V}()
Then, evaluating this:
julia> A,B,C,D,E,F,G,H = create.(basis,Ref(basis))
8-element Values{8, Chain{⟨××××××××⟩, 1, Single{⟨×××⟩, G, B, Int64} where {G, B}, 8}} with indices SOneTo(8):
-1v₁₂₃*v₁ + 1v₂₃*v₂ + -1v₁₃*v₃ + 1v₁₂*v₄ + 1v₃*v₅ + -1v₂*v₆ + 1v₁*v₇ + 1v*v₈
1v₁₂₃*v₁ + -1v₂₃*v₂ + -1v₁₃*v₃ + 1v₁₂*v₄ + 1v₃*v₅ + -1v₂*v₆ + 1v₁*v₇ + 1v*v₈
1v₁₂₃*v₁ + 1v₂₃*v₂ + 1v₁₃*v₃ + 1v₁₂*v₄ + 1v₃*v₅ + -1v₂*v₆ + 1v₁*v₇ + 1v*v₈
1v₁₂₃*v₁ + 1v₂₃*v₂ + -1v₁₃*v₃ + -1v₁₂*v₄ + 1v₃*v₅ + -1v₂*v₆ + 1v₁*v₇ + 1v*v₈
1v₁₂₃*v₁ + 1v₂₃*v₂ + -1v₁₃*v₃ + 1v₁₂*v₄ + 1v₃*v₅ + -1v₂*v₆ + 1v₁*v₇ + 1v*v₈
1v₁₂₃*v₁ + 1v₂₃*v₂ + -1v₁₃*v₃ + 1v₁₂*v₄ + 1v₃*v₅ + -1v₂*v₆ + 1v₁*v₇ + 1v*v₈
1v₁₂₃*v₁ + 1v₂₃*v₂ + -1v₁₃*v₃ + 1v₁₂*v₄ + 1v₃*v₅ + -1v₂*v₆ + 1v₁*v₇ + 1v*v₈
1v₁₂₃*v₁ + 1v₂₃*v₂ + -1v₁₃*v₃ + 1v₁₂*v₄ + 1v₃*v₅ + -1v₂*v₆ + 1v₁*v₇ + 1v*v₈
Finally, we want to take the determinant of this in 8 dimensions, since the matrix will be 8x8 in order to determine the possibility of it leading to an invertible matrix definition:
julia> A∧B∧C∧D∧E∧F∧G∧H
0v*v₁₂₃₄₅₆₇₈
We get a 0 determinant, which means that if we were to attempt defining a matrix for multivectors in a similar way to the one for quaternions you metnioned, then we don’t have the possibility of getting an invertible matrix, unlike as for quaternions:
julia> qbasis = Ref(-I).*Values(v1,v2,v3,v123)
4-element Values{4, Single{⟨×××⟩, G, B, Int64} where {G, B}} with indices SOneTo(4):
-1v₂₃
1v₁₃
-1v₁₂
1v
julia> W,X,Y,Z = create.(qbasis,Ref(qbasis))
4-element Values{4, Chain{⟨××××⟩, 1, Single{⟨×××⟩, G, B, Int64} where {G, B}, 4}} with indices SOneTo(4):
-1v₂₃*v₁ + -1v₁₃*v₂ + 1v₁₂*v₃ + 1v*v₄
1v₂₃*v₁ + 1v₁₃*v₂ + 1v₁₂*v₃ + 1v*v₄
1v₂₃*v₁ + -1v₁₃*v₂ + -1v₁₂*v₃ + 1v*v₄
1v₂₃*v₁ + -1v₁₃*v₂ + 1v₁₂*v₃ + 1v*v₄
julia> W∧X∧Y∧Z
-4v*v₁₂₃₄
As we can see, my calculation comes out compatible with the references, the determinant is not zero which means the resulting matrix would be invertible for quaternions, but for multivectors the same calculation leads to a 0 determinant.
In conclusion, I have demonstrated that what you’re trying to do is not really possible. If you attempt to create a matrix for multivectors that is similarly constructed as the one for quaternions you referenced, then this is not possible based on my interpretations. This technique is possible for quaternions, but does not appear appropriate for multivectors.