Manifolds.jl icon indicating copy to clipboard operation
Manifolds.jl copied to clipboard

tutorial on how to SE2

Open dehann opened this issue 3 years ago • 53 comments

Good day,

Please see a new tutorial for SpecialEuclidean(2). There are things I still want to add, but this should be a good baseline point to clean up as needed and merge. I will add more PRs at a later date.

I'd suggest reading the entire tutorial before staring to write up comments (although I would happily take any feedback!). I suspect this tutorial is coming in from a very different angle than much of the existing documentation. So please look at it on the whole first, and feel free to tear it apart thereafter.

In it's current form, this text and structure is meant for someone who picked up Lie group operations as a mechanism for rigid transforms.

Second last note is that PR #355 should probably be merged first, since this one builds upon that. I'm likely to change some of the variable names in this tutorial to be more in line with those in #355, but also to be as closely recognizable to the existing Manifolds.jl code base.

Lastly, I might have found a bug during this and will open an issue on exp for SE(2) separately to clarify. If so, then I have another draft paragraph to add to the end of this tutorial.

cc @Affie who is likely working with this type of code right now.

Thanks for all the help!

Best, Dehann

dehann avatar May 22 '21 21:05 dehann

Codecov Report

Merging #366 (fa3eed3) into master (99849e7) will increase coverage by 7.20%. The diff coverage is n/a.

Impacted file tree graph

@@            Coverage Diff             @@
##           master     #366      +/-   ##
==========================================
+ Coverage   90.71%   97.92%   +7.20%     
==========================================
  Files          17       76      +59     
  Lines        1012     5915    +4903     
==========================================
+ Hits          918     5792    +4874     
- Misses         94      123      +29     
Impacted Files Coverage Δ
src/ode.jl 100.00% <0.00%> (ø)
src/Manifolds.jl 100.00% <0.00%> (ø)
src/forward_diff.jl 100.00% <0.00%> (ø)
src/SizedAbstractArray.jl
src/Sphere.jl
src/Euclidean.jl
src/PowerManifold.jl
src/Metric.jl
src/Rotations.jl
src/ProjectedDistribution.jl
... and 79 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 99849e7...fa3eed3. Read the comment docs.

codecov[bot] avatar May 22 '21 21:05 codecov[bot]

EDIT: should just mention, I'm documenting converstations related to this tutorial page and thereby understand what are the most important things to add on this tutorial page...

A conversation with @david-m-rosen , @Affie , @pvazteixeira showed we need to be accurate in this documentation on:

  • what is the data representation being used (realization of SE(d) to represent points),
  • what specific metric is being used,
  • what is the convention for the exponential map being used (particular be careful about Lie group exp vs Riemannian exp that corresponds to the metric listed above),

ref

  • https://juliamanifolds.github.io/Manifolds.jl/latest/manifolds/group.html#Manifolds.group_log-Tuple{AbstractGroupManifold,Vararg{Any,N}%20where%20N}
  • https://www.sciencedirect.com/science/article/pii/S0001870876800023

Keeping as breadcrum. I also want to explore / document (pseudo code):

# is p2 is right? ie do i understand exp correct?

p0 = SE2(  [0, 0, 0]  )
p1 = SE2(  [1, 0, pi/2]  )
exp_SE2(p, X) = ....
p2 = exp( p1, [1, 0, pi/2] )

log_SE2(p0, p2) == [1, 1, pi]

exp vs exp_SE2

(This code and test still needs to be refined)

dehann avatar May 27 '21 15:05 dehann

  • what is the data representation being used (realization of SE(d) to represent points),

The usual way we have here is: SE(d) (the manifold type SpecialEuclidean(d) or manifold types in general) are the abstract idea of the manifold, such that manifold_dimension can be implemented. They come with a default idea of points p,q (usually untyped in functions, having AbstractArrays in mind) and a default idea to represent tangent vectors X,Y (on group usually by default represented by elements from the Lie algebra, represented fitting to p,q, also untyped).

For different representations, use subtypes of AbstractManifoldPoint and TVector , respectively. See https://juliamanifolds.github.io/Manifolds.jl/latest/manifolds/hyperbolic.html for an example. And sure, if the default is not well documented, that should be improved.

  • what specific metric is being used,

Similar as for points/vectors the default has a metric in mind, otherwise use new AbstractMetric subtypes and the Metricmanifold{M,G} to distinguish, see https://juliamanifolds.github.io/Manifolds.jl/latest/manifolds/symmetricpositivedefinite.html#Default-metric:-the-linear-affine-metric fo an example.

  • what is the convention for the exponential map being used (particular be careful about Lie group exp vs Riemannian exp that corresponds to the metric listed above),

exp/log mean the Riemannian https://juliamanifolds.github.io/Manifolds.jl/latest/interface.html#The-exponential-and-the-logarithmic-map,-and-geodesics for the Lie group exponential see group_exp https://juliamanifolds.github.io/Manifolds.jl/latest/manifolds/group.html#Manifolds.group_exp-Tuple{AbstractGroupManifold,Vararg{Any,N}%20where%20N}

All these functions (exp/log/group_exp/group_log/retract/...) have as their first argument always the manifold, hence the MetricManifold can also here be used for dispatching on different metrics. Having this in ming, I do not understand your exp_S2? Also p0,p1 would (for the default representation) just be arrays.

For both – different representations and different metric, we would be very happy to extend our code base to more of these; if you are missing a representation or a metric, I‘m happy to help with the structure and even documentation, if you provide the code (that is content within the functions then). :)

kellertuer avatar May 27 '21 16:05 kellertuer

Thanks!

I do not understand your exp_S2? Also p0,p1 would (for the default representation) just be arrays.

I just added a piece of pseudo code to remember and illustrate, one mistake I made was the difference between exp and exp_SE2 -- I was trying this in code and were not quite getting the results I expect it yet. Punchline for me is that a bit more documentation is needed on this rigid transform tutorial to show the difference between which metric and which exponential map is used. Even though I was aware, I still stepped in the trap :-P

After the discussion we had today, the interesting thing about SpecialEuclidean(2) is the subtle difference between:

  • ProductManifold T(n) x SO(n) where exp and log are closer to Riemannian metrics, along with what does group_exp mean in this compared; VERSUS
  • The SemidirectProductManifold i.e. SpecialEuclidean(n) = T(n) x| SO(n) and how to use exp and group_exp. I know these functions are documented individually, and I'd ilke to showcase some easy examples using them.

dehann avatar May 27 '21 18:05 dehann

I'm still not following the discussion from @dehann's example and semi-direct product manifolds.

I wrote the pseudo code in julia as I understand it:

julia> using Manifolds

julia> using StaticArrays

julia> M = SpecialEuclidean(2)
SpecialEuclidean(2)

julia> p0 = ProductRepr(SA[0.0, 0.0], SA[1.0 0.0; 0 1.0])
ProductRepr with 2 submanifold components:
 Component 1 =
  2-element SVector{2, Float64} with indices SOneTo(2):
   0.0
   0.0
 Component 2 =
  2×2 SMatrix{2, 2, Float64, 4} with indices SOneTo(2)×SOneTo(2):
   1.0  0.0
   0.0  1.0

julia> X1 = hat(M, p0, SA[1., 0, pi/2])
ProductRepr with 2 submanifold components:
 Component 1 =
  2-element MVector{2, Float64} with indices SOneTo(2):
   1.0
   0.0
 Component 2 =
  2×2 MMatrix{2, 2, Float64, 4} with indices SOneTo(2)×SOneTo(2):
   0.0     -1.5708
   1.5708   0.0

julia> p1 = exp(M, p0, X1)
ProductRepr with 2 submanifold components:
 Component 1 =
  2-element MVector{2, Float64} with indices SOneTo(2):
   1.0
   0.0
 Component 2 =
  2×2 MMatrix{2, 2, Float64, 4} with indices SOneTo(2)×SOneTo(2):
   6.12323e-17  -1.0
   1.0           6.12323e-17

julia> X2 = hat(M, p1, SA[1., 0, pi/2])
ProductRepr with 2 submanifold components:
 Component 1 =
  2-element MVector{2, Float64} with indices SOneTo(2):
   1.0
   0.0
 Component 2 =
  2×2 MMatrix{2, 2, Float64, 4} with indices SOneTo(2)×SOneTo(2):
   0.0     -1.5708
   1.5708   0.0

julia> p2 = exp(M, p1, X2)
ProductRepr with 2 submanifold components:
 Component 1 =
  2-element MVector{2, Float64} with indices SOneTo(2):
   2.0
   0.0
 Component 2 =
  2×2 MMatrix{2, 2, Float64, 4} with indices SOneTo(2)×SOneTo(2):
   -1.0          -1.22465e-16
    1.22465e-16  -1.0

What I understand from the above code is it is just calling exp on all the parts/components.


What will the manifolds.jl equivalent method be for the following composision using the tangent space?

julia> p0 = SA[1.0 0.0 0.0;
               0.0 1.0 0.0;
               0.0 0.0 1.0];

julia>  RT = SA[0.0  -1.0  1.0;
                1.0   0.0  0.0;
                0.0   0.0  1.0];

julia> p1 = RT * p0
3×3 SMatrix{3, 3, Float64, 9} with indices SOneTo(3)×SOneTo(3):
 0.0  -1.0  1.0
 1.0   0.0  0.0
 0.0   0.0  1.0

julia> p2 = RT * p1
3×3 SMatrix{3, 3, Float64, 9} with indices SOneTo(3)×SOneTo(3):
 -1.0   0.0  1.0
  0.0  -1.0  1.0
  0.0   0.0  1.0

Affie avatar May 28 '21 08:05 Affie

What I understand from the above code is it is just calling exp on all the parts/components. Exactly, since the group manifold is just a decorator: exp/log are not reimplemented for this decorator, so they “pass through” to the underlying product manifold where the (implicit default) metric is the product metric, which decouples.

What will the manifolds.jl equivalent method be for the following composision using the tangent space?

is p0 the same p0 but as one array and then in homogeneous coordinates? What's RT? a tangent vector at p0? An element from the Lie algebra? Another point on the manifold? Maybe you want to apply the group action (matrix multiplication here)?

For the last question with yes, it's

M = SpecialEuclidean(2)
p = ProductRepr([0.0,0.0], [1.0 0.0;0.0 1.0])
q = ProductRepr([1.0,0.0], [0.0 -1.0;1.0 0.0])

then

julia> r = compose(M,p,X)
ProductRepr with 2 submanifold components:
 Component 1 =
  2-element Vector{Float64}:
   1.0
   0.0
 Component 2 =
  2×2 Matrix{Float64}:
   0.0  -1.0
   1.0   0.0

which is what you have, too, see

julia> affine_matrix(M,r)
3×3 Matrix{Float64}:
 0.0  -1.0  1.0
 1.0   0.0  0.0
 0.0   0.0  1.0

(for sure works the same with static arrays, I was just lazy).

kellertuer avatar May 28 '21 09:05 kellertuer

Thanks for all your patience and help while I'm learning more about manifolds and the Manifolds.jl implimentation.

Affie avatar May 28 '21 13:05 Affie

Let me assure you that it is surely not so easy with all these different representations to properly distinguish them (hence all my questions), one has to be very precise there. But great that my explanation helped :) Especially on Lie groups there is quite a set of operations to get right (group operation, exp/log, group_exp/group_log, ...).

kellertuer avatar May 28 '21 13:05 kellertuer

I've seen the plus operator (⊕) used in robotics. This example is how I understand it can be used with Manifolds.jl (we should maybe add something like this in this tutorial?).

julia> using Manifolds

julia> using StaticArrays

julia> M = SpecialEuclidean(2);

# Identity element of SE2
julia> I_SE2 = ProductRepr(SA[0.0, 0.0], SA[1.0 0.0; 0 1.0]);

julia> p0 = ProductRepr(SA[0.0, 0.0], SA[1.0 0.0; 0 1.0]);

julia> X1 = hat(M, p0, SA[1., 0, pi/2]);

julia> p1 = exp(M, p0, X1);

julia> X2 = hat(M, p1, SA[1., 0, -pi/2]);

# right ⊕ : p2 = p1 ⊕ X2
julia> p2 = compose(M, p1, exp(M, I_SE2, X2))
ProductRepr with 2 submanifold components:
 Component 1 =
  2-element MVector{2, Float64} with indices SOneTo(2):
   1.0
   1.0
 Component 2 =
  2×2 MMatrix{2, 2, Float64, 4} with indices SOneTo(2)×SOneTo(2):
   1.0  0.0
   0.0  1.0

# left ⊕  :  p2 = X2 ⊕ p1
julia> p2 = compose(M, exp(M, I_SE2, X2), p1)
ProductRepr with 2 submanifold components:
 Component 1 =
  2-element MVector{2, Float64} with indices SOneTo(2):
    1.0
   -1.0
 Component 2 =
  2×2 MMatrix{2, 2, Float64, 4} with indices SOneTo(2)×SOneTo(2):
   1.0  0.0
   0.0  1.0

# exp map on product manifold components
julia> p2 = exp(M, p1, X2)
ProductRepr with 2 submanifold components:
 Component 1 =
  2-element MVector{2, Float64} with indices SOneTo(2):
   2.0
   0.0
 Component 2 =
  2×2 MMatrix{2, 2, Float64, 4} with indices SOneTo(2)×SOneTo(2):
   1.0  0.0
   0.0  1.0

What got me: I was expecting p2 = exp(M, p1, X2) to look like p2 = p1 ⊕ X2

I also still need to figure out how group_exp works and fits into this (think I'm missing some fundamentals).

Please let me know what you think.

Affie avatar Jun 04 '21 10:06 Affie

I've seen the plus operator (⊕) used in robotics.

The question is: How is it defined? Having seen it – does not give (me) the context. So for me it is unclear “what got you”, since I am missing the context

For exp keep in mind it is the component wise exponential map, so for the matrix component its the exponential map on the rotation matrices, for the R^2 component (of the translation Group thereon) it is just +. This is both exactly, what the exp does there.

kellertuer avatar Jun 04 '21 12:06 kellertuer

Sorry I could have been clearer, it is defined in the code comment:

p2 = X2 ⊕ p1 as p2 = compose(M, exp(M, I_SE2, X2), p1)

I'm just making a statement of my pitfalls coming from engineering.

For exp keep in mind it is the component wise exponential map, so for the matrix component its the exponential map on the rotation matrices, for the R^2 component (of the translation Group thereon) it is just +. This is both exactly, what the exp does there.

Yes, this is how I now understand it.

Affie avatar Jun 04 '21 14:06 Affie

HI, perhaps I can ask a related question -- I'm trying to understand how to use exp/retract to "walk" through a series of points on a product manifold (in the general sense). I have vector Xa which is defined in the tangent space around p1, and Xb which is defined in the tangent space of p2:

manifoldsWormWalk

neither group_exp or exp for SpecialEuclidean(2) seem to do the operation shown in the figure above...

  • I think the confusion comes from which exp is being used when (and if there might be an implementation gap).
  • Can I use exp/retract to do these local operations, or must I first use exp defined against group identity elements and use then a compose operation?

i.e. retract with vectors local to each point VS compose using group elements is the only answer?

dehann avatar Jun 04 '21 20:06 dehann

You are again a little imprecise here. So if you have p1 as a point on the manifold and Xa is a tangent vector at p1 then sure exp)M,p1,Xa) is a point on the manifold, as is retract(M,p1,Xa,r) if you have the retraction r defined (note that exp is a special case of a retraction). For the definition of a retraction take a look at http://sma.epfl.ch/~nboumal/book/index.html (p. 46, Sec 3.6 or precisely Def 3.41). And exp fulfils some ode additionally to that (Sec 10.2, p. 242, precisely Def 10.13).

Can you explain why you think exp is not doing such an operation? It really does!

What you have to be careful about is whether (for group manifolds/Lie groups) your Xa is a tangent vector or a Lie group element. My text above only holds if it is a tangent vector as you stated! If you have an element from the Lie group you first have to use the differential of the group operation, see for example https://juliamanifolds.github.io/Manifolds.jl/latest/manifolds/group.html#Manifolds.translate_diff-Tuple{AbstractGroupManifold,Vararg{Any,N}%20where%20N}

For the group exponential, that is a little different because it has as input an element from the Lie algebra (and not from the tangent space) see https://juliamanifolds.github.io/Manifolds.jl/latest/manifolds/group.html#Manifolds.group_exp-Tuple{AbstractGroupManifold,Vararg{Any,N}%20where%20N}

kellertuer avatar Jun 04 '21 20:06 kellertuer

Can you explain why you think exp is not doing such an operation? It really does!

@kellertuer , I think we are close to breaking through on this. Here is a minimum example showing my confusion on exp about some random user point p.

M = SpecialEuclidean(2);
p1 = ProductRepr(SA[0.0, 0.0], SA[1.0 0.0; 0 1.0]);
Xa = hat(M, p1, SA[1., 0, pi/2]);
p2 = exp(M, p1, Xa);
Xb = hat(M, p2, SA[1., 0, -pi/2]);
p3_test = exp(M, p2, Xb)
# \theta=0.0 (works okay in this example)
@assert isapprox( [1 0; 0 1], p3_test.parts[2])

# but why t=[2,0]?
p3_test.parts[1]' |> println
# should be [1,1]?
@assert isapprox( [1.0,1.0], p3_test.parts[1])

(simplified code and similar to the picture just above)

So if you have p1 as a point on the manifold and Xa is a tangent vector at p1 then sure exp)M,p1,Xa) is a point on the manifold

I would expect this to be point p2, which we can then use with a new vector Xb to "walk" to p3, no? That's what this code snippet is trying to do.

All I can think of is we are using the wrong exp for SpecialEuclidean(2) somehow, maybe but I'm not sure... I checked group_exp(M, X) but that is not defined about a user point (I think it is assuming the group identity element for it's coordinate frame). Perhaps there is still an implementation gap around group_exp, or I'm just missing a major piece of fundamental info. Either way, I think something about this should be added to this tutorial page.

dehann avatar Jun 06 '21 14:06 dehann

So if you have p1 as a point on the manifold and Xa is a tangent vector at p1 then sure exp)M,p1,Xa) is a point on the manifold

I would expect this to be point p2, which we can then use with a new vector Xb to "walk" to p3, no? That's what this code snippet is trying to do.

Let me phrase that carefully: Sure exp(M,p1,Xa) is some point on the manifold, and sure if at this some point you move somewhere further, then you end up at some third point.

Keep in mind that in the first (R2 component) exp is really just addition! So with that, p1 has [0,0], adding Xa ([1,0]). and Xb(again [1,0]) of course yields [2,0], because that space is - for exp – really just the Euclidean R2.

Why I was asking about more details is exactly: What do you expect, where does this expectation come from? All these details I am missing, so I can just guess and not really help. So with just code without any explanation, I don‘t know why you think it should be [1,1]because all I can see in the Code – it is perfectly fine and computes the correct things for the input you provide and the functions you use.

So again: Can you explain why you think exp is not doing such an operation? It really does!

References, you own computations by hand, something like that is what I would call more detail.

All I can think of is we are using the wrong exp for SpecialEuclidean(2) somehow, maybe but I'm not sure... I checked group_exp(M, X) but that is not defined about a user point (I think it is assuming the group identity element for it's coordinate frame). Perhaps there is still an implementation gap around group_exp, or I'm just missing a major piece of fundamental info. Either way, I think something about this should be added to this tutorial page.

Watt is a “user point” here? Which coordinate frame do you have in mind?

kellertuer avatar Jun 06 '21 15:06 kellertuer

I think I know what's going on here. I guess what you want to do is to have an object on a plane with an orientation. Let is start at your p1:

p1 = ProductRepr(SA[0.0, 0.0], SA[1.0 0.0; 0 1.0]);

Let us also make the identity element of SE(2) for convenience:

id = ProductRepr(SA[0.0, 0.0], SA[1.0 0.0; 0 1.0])

Now, we want to apply two actions to this object: translation by 1 in the direction it's facing and rotation by pi/2:

a1 = exp(M, id, hat(M, id, SA[1.0, 0, pi/2]))

The other operation is similar but performs rotation by pi/2 in the other direction:

a2 = exp(M, id, hat(M, id, SA[1.0, 0, -pi/2]))

Now, to do the walking you seem to expect (I think?) we need right action of the group operation:

GOA = GroupOperationAction(M, RightAction())

and we can start walking by applying our operations:

julia> u = apply(GOA, a1, p1)
ProductRepr with 2 submanifold components:
 Component 1 =
  2-element MVector{2, Float64} with indices SOneTo(2):
   1.0
   0.0
 Component 2 =
  2×2 MMatrix{2, 2, Float64, 4} with indices SOneTo(2)×SOneTo(2):
   6.12323e-17  -1.0
   1.0           6.12323e-17

julia> u2 = apply(GOA, a2, u)
ProductRepr with 2 submanifold components:
 Component 1 =
  2-element MVector{2, Float64} with indices SOneTo(2):
   1.0
   1.0
 Component 2 =
  2×2 MMatrix{2, 2, Float64, 4} with indices SOneTo(2)×SOneTo(2):
   1.0  0.0
   0.0  1.0

Does that seem right?

mateuszbaran avatar Jun 06 '21 15:06 mateuszbaran

You probably need to be a bit more careful in distinguishing elements of SE(2) corresponding to location/orientation of an object and elements of SE(2) that correspond to operations (translations, rotations) of said object.

mateuszbaran avatar Jun 06 '21 15:06 mateuszbaran

All these details I am missing, so I can just guess

Understood. I think there are two sets of conventions here which makes it tricky to connect. What I'd like to do is use the Manifolds.jl way (which I'm taking as the right way) two write down this "walking from point to point on a general manifold". Think I'm using Manifolds.jl incorrectly.

I don‘t know why you think it should be [1,1] because all I can see in the Code – it is perfectly fine and computes the correct things for the input you provide and the functions you use.

Yes, I'm expecting translation from [1,1] from the code snippet ... This helps me pin down that exp is not doing what I expect, thanks.

want to do is to have an object on a plane with an orientation.

The reason I expect [1,1] is easiest to explain via cartoon avatar on a XY plane, and tracking a yaw angle for the avatar, hence SE(2). Starting from the coordinates for p1 at [x,y,yaw] = [0,0,0], walk to p2 using vector coordinates [1,0,pi/2] (i.e. walk forward and turn 90deg left, assuming z-axis points up). What I called Xa above. Now walk from p2 to p3 with new vector coordinates [1,0,-pi/2], what I called Xb above. Note, "forward" from p2 is actually along the y-axis in the "global coordinates". I'm taking "global coordinates" here as the origin of the XY plane where we started. In "global reference" I would then expect the coordinates of p3 to be [1,1,0].

366_walk

The question though is whether this can be done with exponential maps directly from vectors (attempted here), or whether it can only be done with actions on group elements?

The latest example (https://github.com/JuliaManifolds/Manifolds.jl/pull/366#issuecomment-855419533) is still actions using group elements. So, it looks unlikely to me that the idea of "walking" on a plane with orientation (which we represent with elements from SE(2)) can be done using vectors and the exponential map directly (as i tried to show here). It is looking more likely that we have to use compose or actions where we first use exp to get a group element around the identity and then as a second stage use compose / actions to effect the rigid-transform.

Note, the group action / compose version is already in the PR. Thanks, I know this has taken a lot of time. I will try write all this succinctly at the end of this How to SE(2) documentation page.

dehann avatar Jun 06 '21 16:06 dehann

The answer is already directly above (1) turn your Lie algebra elements (tangent vectors at the identity) into points on the manifold and apply the group action then. That was what Was asking and Mateusz guesses correctly. Your Xa, Xb are tangent vectors at the identity not at p1,p2.

kellertuer avatar Jun 06 '21 16:06 kellertuer

The question though is whether this can be done with exponential maps directly from vectors (attempted here), or whether it can only be done with actions on group elements?

The latest example (#366 (comment)) is still actions using group elements. So, it looks unlikely to me that the idea of "walking" on a plane with orientation (which we represent with elements from SE(2)) can be done using vectors and the exponential map directly. It is looking more likely that we have to use compose or actions where we first use exp to get a group element around the identity and then as a second stage use compose / actions to effect the rigid-transform.

After looking at your example, I think group operation actions is the right abstraction here. There are, although, other ways around it. I think I could come up with something relying on differentials of compose/apply, and then using FVector we could add overloads that would let you work on SE(2) with the compact translation-angle representation if you need that but what I wrote above is a much simpler solution. So you don't really have to use exp but this way is easier.

mateuszbaran avatar Jun 06 '21 16:06 mateuszbaran

Alright, great thanks. Let's then just use groups. I'll add to the end of this new documentation page a note to show that

  • Manifolds.exp(M,p,Xa) is not the same as compose(M, p, exp(M, I_SE2, Xa)) (linking to this #366 discussion), and
  • The second example of using actions, as per https://github.com/JuliaManifolds/Manifolds.jl/pull/366#issuecomment-855419533

My next steps are to make sure all the updates are in 355 and 366 so you can merge when happy. Thanks again for all the feedback!

dehann avatar Jun 06 '21 17:06 dehann

Great, looking forward to the update. I would prefer though, if it was self-contained and that it does not need links to comments/issues.

kellertuer avatar Jun 06 '21 17:06 kellertuer

I've looked into group operations in Lie algebras and it's actually something we don't have yet in Manifolds.jl: https://encyclopediaofmath.org/wiki/Adjoint_representation_of_a_Lie_group . We probably need to introduce Lie brackets as a separate operation, probably lie_bracket(G::GroupManifold, X, Y) with a default implementation like

function lie_bracket(G::SpecialEuclidean, X, Y)
    mX = screw_matrix(G, X)
    mY = screw_matrix(G, Y)
    bracket = mX*mY-mY*mX
    return from_screw_matrix(G, bracket)
end

Actually, maybe it'd be a good idea to include it in this tutorial?

@sethaxen Did you think before about Lie brackets for SpecialEuclidean?

mateuszbaran avatar Jun 07 '21 08:06 mateuszbaran

I just found this article and it's quite helpful: https://arxiv.org/pdf/1812.01537.pdf. It also gives the differences I tried to show in https://github.com/JuliaManifolds/Manifolds.jl/pull/366#issuecomment-854613353 clearer (p11): image

Affie avatar Jun 07 '21 18:06 Affie

I just found this article and it's quite helpful: https://arxiv.org/pdf/1812.01537.pdf. It also gives the differences I tried to show in #366 (comment) clearer (p11):

Thanks for this link, it's quite helpful, but note that it has a slightly different approach to Lie groups than Manifolds.jl. As far as I can see our SE(n) is what they call T(n)xSO(n). I don't know what exactly their ρ is though. Well, there are some equations from which it could be deduced but I don't see any justification. I tried looking in the book [11] they cite but they use a different notation, and page 35 talks about parametrizations (we now have support for parametrizations so we could add that one).

mateuszbaran avatar Jun 07 '21 19:06 mateuszbaran

Yeah, my current understanding is the difference between

  • T(n) x SO(n) a direct product manifold (aka "more Riemannian") and I think Manifolds.jl in general is coming from this angle vs.
  • T(n) x| SO(n) as SemiDirectProduct (aka Lie Algebra SE(n) as used in rigid transform write-ups), but also as per current docs.
  • (and new reference suggests a third option which I still need to digest).

The differences in convention makes it difficult to connect and communicate quickly, hence this "How to Rigid-Transforms" example page...

dehann avatar Jun 07 '21 19:06 dehann

Thanks for this link, it's quite helpful, but note that it has a slightly different approach to Lie groups than Manifolds.jl. As far as I can see our SE(n) is what they call T(n)xSO(n).

I think this difference is where the confusion came from (at least for me).

Going by the component wise way manifolds.jl works with Lie groups, I would say the <R^n, SO(n)> is the one implemented?

Affie avatar Jun 07 '21 19:06 Affie

  • T(n) x SO(n) a direct product manifold (aka "more Riemannian") and I think Manifolds.jl in general is coming from this angle vs.

This is N = ProductManifold(TranslationGroup(3), SpecialOrthogonal(3)) and exp/log on this is the third point

  • T(n) x| SO(n) as SemiDirectProduct (aka Lie Algebra SE(n) as used in rigid transform write-ups), but also as per current docs.

Yes, that‘s M = SpecialEuclidean(3). Note that its type is more precisely

GroupManifold{ℝ, ProductManifold{ℝ, Tuple{TranslationGroup{Tuple{3}, ℝ}, SpecialOrthogonal{3}}}, Manifolds.SemidirectProductOperation{RotationAction{TranslationGroup{Tuple{3}, ℝ}, SpecialOrthogonal{3}, LeftAction}}}

So it decorates the manifold N above with a Semidirect group action that consists of a couple of group actions.

  • (and new reference suggests a third option which I still need to digest).

The differences in convention makes it difficult to connect and communicate quickly, hence this "How to Rigid-Transforms" example page...

kellertuer avatar Jun 07 '21 20:06 kellertuer

  • T(n) x SO(n) a direct product manifold (aka "more Riemannian") and I think Manifolds.jl in general is coming from this angle vs. Direct vs semidirect product only refers to the group operation, and I'd say that neither is "more Riemannian" than the other. One could of course talk about metric invariance under group operation action but it's a different story.

Going by the component wise way manifolds.jl works with Lie groups, I would say the <R^n, SO(n)> is the one implemented?

Riemannian structure goes component-wise, yes, but the group structure does not, and we can add parametrizations that also mix components so it's actually complicated... these plus operations in circles or squares seem to mix group composition and parametrization, where pluses for (their) SE(n) and T(n)xSO(n) have the same group composition (the semidirect product one I think) but different parametrizations. The <R^n, SO(n)> one has a direct product composition.

So, in Manifolds.jl, to compute their plus you'd have to:

  1. go from matrix representation of M to ProductRepr representation (at least for now, can be extended later),
  2. go back from parametrization of tau to a ProductRepr group element using get_point,
  3. do apply of group operation action.

Or at least that's one way, one could also use parametrization of the Lie algebra and add exp after than instead of point (2).

mateuszbaran avatar Jun 07 '21 20:06 mateuszbaran

Or, still alternatively, you could define in Manifolds.jl those plus operations as actions of a yet another group on SE(n). There is just way too many ways to describe the same computation. Still, as far as I can see they key part of those pluses is our apply.

mateuszbaran avatar Jun 07 '21 20:06 mateuszbaran