Advice on implementing a new group...

[domo-mr-roboto]

Hi all - thanks for the nice work on ManifoldsBase.jl and Manifolds.jl. I'm in the process of trying to determine how best to implement a new group - specifically, the "group of extended poses," which is in fact the semidirect product of SE(3) and $\mathbb{R}^3$. I began by defining the semidirect product of SE(3) and $\mathbb{R}^3$ as a SemidirectProductGroup, but this seems to lead to quite complex code that's not easy to parse (by me).

Notably, the group can be embedded in GL(6) and that's how I will use it, essentially exclusively. I'd like to develop a lightweight, matrix-based implementation - unsure whether to define a new group using Manifolds.jl (also unsure if there's a tutorial on implementing a new group?) or to try to use ManifoldsBase.jl only. Any advice would be appreciated; the best possible speed is important, and given the fixed size of the matrices involved, I intend to use StaticArrays.jl.

Also, a rather silly q (as a newbie) - I'm uncertain why the representation size is given as the manifold dimension + 1? What am I missing? Calling representation_size() on SO(2) gives (2, 2), so I'm somewhat confused.

J.

[kellertuer]

For now we indeed do not yet have a tutorial on how to implement a (Lie) group – but hopefully the ManifoldsBase.jl tutorial on a manifold is a good starting point – e.g. on how our functions are usually meant to be implemented. this includes the important difference between the allocating and in-place (fast!) variants we usually have in mind.
Then the code for one of the groups a good continuation.
Sure, StaticArrays should not be a problem, neither points nor tangent vectors are usually typed (more than necessary) which is exactly to be able to use other arrays formats if you would like to.

If you feel it is easier, then do not go the SemidirectProductGroup but implement your own group in the embedding to GL(6) as you write.

For the representation size – that is the array size used to represent points and tangent vectors for a certain manifold – which is unrelated to the dimension.
Example: For the sphere, your observation is correct, the dimension of the sphere is n and then its points are represented as vectors (arrays) of length n+1 (unit vectors to be precise).
SO(2) is represented as 2x2 rotation matrices, hence its representation size is (2,2)
The symmetric positive definite matrices of size $n\times n$ are a manifold of dimension $\frac{n(n+1)}{2}$, so here representation_size(M) would be (n,n) while manifold_dimension(M) returns n*(n+1)/2.

Feel free to ask further questions along the way and if you like, a PR for new manifolds and/or groups is always welcome.

[domo-mr-roboto]

Super thanks! As I become more familiar, I might be able to write up a tutorial. That makes good sense re: the sizing.

[kellertuer]

A tutorial on how to implement a group would be very cool indeed.

[mateuszbaran]

Hi! We don't currently have a tutorial for Lie groups but I think the best way to go is for example copy the implementation of HeisenbergGroup and modify methods to fit your group. SemidirectProductGroup is fairly advanced and not fully developed solution so I think it's fine to ignore it.

Representing that group by 6x6 matrices is one way to go but it might not be totally efficient. You can also consider using ArrayPartition of three static arrays, and derive optimized implementation using the embedding in GL. That's what I sometimes do when implementing methods for SE(n). You can also take a look how the ArrayPartition is used for special euclidean group.

Manifolds.jl interface generally can work with both SArray and MArray, though we don't always have the necessary specialized non-mutating methods for SArray -- feel free to report if you need some. We add them when there is demand.

[mateuszbaran]

By the way, I'm currently slowly moving towards robotics and IMU integration, so there will be more support for it in Manifolds.jl and the ecosystem (see RoME.jl). I think the standard Lie group approach doesn't map to the right mathematical concepts, which appear to be vector bundles instead. I can explain more if you are interested.

[domo-mr-roboto]

Great thanks! Good point about ArrayPartition - I'll investigate. And absolutely, would be happy to hear more about IMU preintegration (something else I'm interested in).

[mateuszbaran]

There has been a recent work and discussion here that you can read for a start: JuliaRobotics/RoME.jl#701 . IMUDeltaGroup is very similar to the group of extended poses.

[domo-mr-roboto]

Ah great - very helpful! (and sorry for late reply) Also I realized I mistyped above in the original post (minor detail) - it should be an inclusion in GL(5).

[domo-mr-roboto]

Following up on the above after far too long! I've continued to work with Julia and Manifolds.jl in my spare time (never sufficient). Hoping it's alright to slightly diverge from the original topic.

I'm still working with various Lie groups and looking at the LieGroups.jl package - notably, I'm in engineering and we're used to working with matrix Lie groups. I wanted to ask about one of the design decisions for Manifolds.jl and LieGroups.jl, namely, passing a 'group' instance along with all calls to exp, log, etc.

Typically, in our engineering work, it's most useful/easiest to deal with Lie group objects, so instead of having to call:

X = log(SO3, g)

it would be more convenient to call:

X = log(g)

where g is an element of SO(3) and a subtype of AbstractArray (or similar - note that I haven't fully thought this through). The approach is similar to what Rotations.jl does with quaternions (i.e., you don't need to pass an argument to identify a quaternion, you just pass an instance of the quaternion type).

Perhaps there are use cases that I'm not familiar with or design choices that motivate the approach? I'd like to start building on top of LieGroups.jl, but want to make sure I understand its design first.

Thanks!

[domo-mr-roboto]

Thanks to you both for the detailed and thoughtful responses! This makes good sense - particularly the need to produce a large number of new types. In some sense, I'm coming from a 'biased' use case, where only a small number of Lie groups (maybe four or five at most) are used. A couple of follow-up questions on two points (unintended pun):

1. For rand(), I certainly like the idea of calling rand(SE3) to get a random element - I wonder though does this in some ways conflict with the way the base rand() function is defined? That is, we call rand(3, 3) to get a 3 x 3 array, and not rand(Array{Float64, 2}, 3, 3) - the base rand() does not require a 'type' argument. A call to get a random integer within some range, rand(1:10) say, also infers the required output from the single input (i.e., dispatch based on UnitRange{Int64}). Indeed though, I can't see a clean approach to do this in some other way for groups (could do rand(SE3IdentityElement()), that is, just pass any element of the group, but that's quite ugly).

2. In terms of memory, is there overhead for an array that's all of the same type? That is, if I have an array of 1,000,000 points on a Lie group (as in the example), only the type of each point and the data itself need to be stored, which would be the same as storing any other type and data? Although there might be padding needed if the type isn't a primitive type (depending on memory layout, I'd imagine). I just want to make sure I'm understanding the overhead issue.

My main thought related to the group/object being passed as the first argument to associated calls is that when building estimators (for, e.g., tracking the position of a mobile robot or similar), I end up writing code that might look something like:

H = inv(Jᵣ(G, err))*Ad(G, inv(G, F))

(inverse of right Jacobian wrt to error term, times adjoint of element of group G, where the element is the inverse of the element F...) I find myself sometimes wanting to write:

H = inv(Jᵣ(err))*Ad(inv(F))

Or, in an object-oriented language, something like:

H = Jr(err).inv()*F.inv().adjoint()

But I do certainly understand the desire for the clean separation of the manifold/group logic from the representation (etc.) and the need to support a much wider range of operations (connections, etc.).

Actually (all the way back in late 2023), the group I was working on implementing in Julia (outside of Manifolds.jl at the time) was the Galilean group, which is actually exactly the IMUDeltaGroup mentioned above, it turns out.

I have a working implementation of the Galilean group, but built entirely with matrices (very naively), mainly because I originally need to hack something together quickly'ish. Maybe my next step can be to try convert this to use the LieGroups.jl framework And thanks again for the effort to put these very useful packages together!

[kellertuer]

For the first – not really or fully. Or: It depends ;)

1. There is rand(Float64, (3,3) where you can specify the “element type” of the array to be Float64. and rand(Float64) does just return a float. So sure passing a manifold and getting a point thereon is a bit of a stretch, but rand would need something like the size information (are we on the 2-sphere or the 3-sphere?) anyways, which is beyond something like just a float type.
   There is Identity(SE3) which is an indicator (but memoryless) for the identity on SE3 but calling rand on that and expecting a point on SE3 really seems a bit stange.
   I would consider the idea here to be closer to the Float64 idea, i.e. specifying the one element you want to have, with the small caveat that we provide the manifold the point is on and not the concrete point type. The reason is that the point type would just be arrays. So maybe easier on manifold
   rand(Sphere(2)) returns an arbitrary unit vector in R3 – but as a default array.
   and sure we follow the interface of rand, that you alwaus are allowed to pass an rng as first parameter.
   One thing we could do is to also accept numbers and then return an element on the corresponding power manifold. But since that one also has at least 2 different representations (cf. plan vs nested) that is also not so easy.

2. In the simple scenarios you have in mind, the point types might be quite lightweight, for best of cases without overhead.
   But if the manifold carries some extra information – see e.g. the MetricManifold that allows to use different metrics on a manifold – then we have to copy that information which metric to use to every point, since there is no “central” place to store the metric. For worst of cases, the metric might consist of storing a d-by-d(d the manifold dimension) SPD matrix.
   Sure if you do that careful, you would just reference this matrix in every point and not copy it – but you would have to be veeeeery careful with that not to store that metric of potentiall d^2 numbers too often.
   Also, if I hand you, say, a matrix where every entry is meant to be a point on the sphere, but it came from another application/pckage – so it is a matrix A where every entry A[i,j] is a unit vector itself (so something like Array{Array{Float64,1},2} or so).
   In the current implementation you just “throw” that into the power manifold call and you are fine, or with single enties you just walk into the sphere and you are fine.
   If all of those have to be ypoed, yo first have to all SpherePoint(a) for a in A or even a PowerManifoldPoint around that? Gets a bit technical. Sure points can be typed in our ecosystem to distinguish different representations, but the most common / eastiest one is usually untyped.

I always confuse your right/left/center Jacobians.
For me – or maybe even matheatically there is only one, namely the Jacobian as the collection (or matrix to be precise) of derivatives of a multivariate vector-valued function, but yes our adoint (https://juliamanifolds.github.io/LieGroups.jl/stable/interface/group/#Base.adjoint-Tuple{AbstractLieGroup,%20Any,%20Any}) has the Lie group upfront, and inv (https://juliamanifolds.github.io/LieGroups.jl/stable/interface/group/#Base.inv-Tuple{AbstractLieGroup,%20Any}) as well, where (again) we might have different views on what an adoint is, ours might be the general case, where you have to set g=e to get yours.

An advantage of this approach is, that we have a inv!(G, h, g) where we compute the inverse of g in the memory of h avoiding unnecessary allocations. So also inv!(G,g,g) is possible, to “overwrite” g with its inevrse in-place.
Not only is avoiding allocations one thing that makes Julia fast – on large Liegroups (again in imaging imagine a million pixel of a rotation each as a “superLieGroup”) it might be crucial not not copy all of that.

Overall, I personally see more advantages in splitting that and having the interface we currently offer. Even after now quite a few years working with that.
Your mileage might vary, but then you simply do not fit the idea we follow here. I am convinced that any package with certain modelling decisions might mean that others do not agree with those decisions.
Then maybe start your own package.
If you want exactly the style you write (H = inv(Jᵣ(err))*Ad(inv(F))) – then you have to do your own CoolLieGroups.jl package. But as you wrote, setting up those packages is a bit of an effort – and we discussed the overall design ideas for quite a while (and had a start with types points indeed).

For the Gallilean Group – that was brought up when we started with LieGroups, cf JuliaManifolds/LieGroups.jl#11, but since I for example do most of this in my free time and my university job keeps me busy as well – some things might only happen in the more distant future. But if you want to contribute something like that I would not only be very happy but also be willing to help along the way, review the PRs and help with the docs.

[mateuszbaran]

1. For rand(), I certainly like the idea of calling rand(SE3) to get a random element - I wonder though does this in some ways conflict with the way the base rand() function is defined?

It does not because we require manifolds to subtype our type (AbstractManifold), so it's just a bit of stretching the standard semantics of rand. We were quite careful to not introduce any possible ambiguities. Note that Julia also has sampling from an Array for example (rand([1, 3])), so the first argument of rand can already be quite diverse in Julia.

My main thought related to the group/object being passed as the first argument to associated calls is that when building estimators (for, e.g., tracking the position of a mobile robot or similar), I end up writing code that might look something like:

H = inv(Jᵣ(G, err))*Ad(G, inv(G, F))

I think you might be interested in our new library: https://github.com/JuliaManifolds/GeometricKalman.jl (and the associated preprint linked there). It provides a more thought-through description of Kalman filtering than the common Lie group approach. In particular, I advocate against using the terms "left Jacobian" and "right Jacobian". There is really no symmetry there, it's just one standard Jacobian of two different functions. Moreover, those Jacobians quite often have structure that can be used to speed up computations, or we might want to consider doing in-place computation of the Jacobian, so optimized code wouldn't look short and neat anyway.

[domo-mr-roboto]

Great, thanks to you both! (and sorry for my delayed response)

I'm convinced re: the use of the manifold / group argument 😄. And quite right about the Jacobian naming - I'm just using the naming convention that seems to have become prevalent in the estimation literature.

Oh! I didn't realize that the Galilean group had already been implemented in RoME.jl - that's great. I'm happy to work on porting this over to LieGroups.jl if it's useful. Also, is a tutorial still needed?

GeometricKalman.jl looks very interesting as well - I'll check that out too!

[kellertuer]

No worries, this is never urgent!

Actually, @Affie said on the last Juliamanifolds community call, he might give it a try as well, so maybe check with him as well.

Yes for now there is the interface pages for the Lie group (https://juliamanifolds.github.io/LieGroups.jl/stable/interface/group/) and algebra (https://juliamanifolds.github.io/LieGroups.jl/stable/interface/algebra/), but a tutorial “How to implement a Lie group” similar to “How to define a manifold” https://juliamanifolds.github.io/ManifoldsBase.jl/stable/tutorials/implement-a-manifold/ would be a great thing to have, I think.

Actually any tutorial illustrating certain uses of Lie groups is great :)

For the time delay, no worries, this is all not urgent (from our side at least).