Towards a common Manifold Framework / Package

[antoine-levitt]

CC @kellertuer, @mateuszbaran, @pkofod

This is to continue the discussion from https://discourse.julialang.org/t/ann-manopt-jl/24906/17

[antoine-levitt]

Regarding types: I considered doing that (eg having a different type for points on the manifold and tangent vectors), but opted against it as over-engineering, so I have everything as just plain arrays. My motivation at the time was inserting manifolds into an existing package for optimization, and so the simplest things were, the happier I was. I still think this is a good idea: types are good, but complicate things (eg should they be AbstractArrays or not, what underlying array type should they use, etc.) So my prefered API would be an extension of the current one, with things like

retract(M::Manifold, x)
exp(M::Manifold, x, v)
parallel_transport(M::Manifold,x, y, v)

etc, where the x, y, v would simply be arrays.

Does this look good or is it too limitative for your use cases?

[antoine-levitt]

(I moved the discussion here for convenience; doesn't mean that we should keep the name or the code)

Another use case that would have to fit inside such a package: implicit manifolds defined by g(x)=0, using a nonlinear solver for retractions etc.

[kellertuer]

I do dispatch on all three types and can check dimensions and such based on types. I did an implementation in Matlab with matrices and felt much much much better with proper types, also to abstract from a matrix or vector representation, since I might want to do shape manifolds and such, where the inner representation can be a graph.

For implicit manifolds (in my mind) the manifold type would hold the information (g(x)=0) and the points and vectors are again representations of solutions of the solvers (independent on how you would like to solve that).

I know that's quite an abstract way to look at manifolds, points and tangent vectors (I like to also include cotangents).

To the moving of the discussion – I'm fine with that, I wanted to avoid exactly keeping the name of something, but with that cleared – I am fine and happy to lead the discussion (independent of where we do that).

Edit: Of course one could do converts from vectors to points on the sphere and vice versa, so with that you could still call the functions as you like and I'd still have my proper types

[mateuszbaran]

About types, I think a common interface wouldn't really impose that many restrictions on them. It's mostly about specifying a few abstract types and method stubs that formalize the interface (which can be implemented in any way). Manopt.jl and FunManifolds.jl already introduce them in a similar way. They could be implemented in any way. I also agree that having proper types makes working with points and tangent vectors much easier.

There is also the matter of argument order. Julia style guide suggests that mutated input should go very early in the order, so retract!(m::Manifold, x) does not conform. Do we break the rule here or not?

@kellertuer Could you give an example of what do you mean by convert from vector to point? Also, what is the benefit of specifying manifold in, for example, exp or log where other arguments have specific types? For example here.

BTW, let's also link the other thread from discourse: https://discourse.julialang.org/t/towards-a-common-manifold-framework-package/24910

[mateuszbaran]

One more point: efficient implementation of product manifolds (and power manifold, functional manifolds and some other) is very tricky. ManifoldProjections.jl uses array reshapes, Manopt.jl uses simple Arrays. Both approaches result in suboptimal performance. In FunManifolds.jl I use AbstractArray-wrapped tuples to be able to get isbits types of points on product manifolds (which in turn I then need to have efficient functional manifolds as Arrays of isbits types). The trouble with isbits types, though, is that they are immutable and to work around that I had to invent that condbc macro.

This is mainly to show what things should fit in the general framework, my implementation is definitely too specialized for a general package.

[antoine-levitt]

My concern with having points/vectors being types is that it forces you to do everything in that framework, adds a lot of complexity and is pretty inconvenient for playing around/debugging: you can't do println(norm(x-y)), or plot(x) or whatever; in Manopt.jl's case, you have to have a lot of getValue & co. The benefit is mainly subtractive: forbidding the user to do mathematically forbidden operations: having types for point and vectors don't really get you much more functionality. As a user, I usually find this kind of limitations extremely annoying. I feel like Julia encourages its users to avoid this sort of structure-heavy design, and that's one of the reasons I like it so much.

Possibly a good solution would be to have a package with a low-level API, in the style retract(M::Manifold,x::AbstractArray) etc. Then, it would be easy to design a higher-level API with types that would call into this lower-level API.

ManifoldProjections.jl uses array reshapes, Manopt.jl uses simple Arrays. Both approaches result in suboptimal performance

Really? Reshapes and views used to suck a lot pre-1.0 but seem fine now. At least I've had good experiences on pretty large optimization problems.

I was never satisfied with the design for power/product manifold, but I needed something quick and dirty for a project.

There is also the matter of argument order. Julia style guide suggests that mutated input should go very early in the order, so retract!(m::Manifold, x) does not conform. Do we break the rule here or not?

That's true; I kinda like the convention of putting dispatch operations first (so that every manifold-related function starts with a manifold), but don't have strong feelings about it.

[kellertuer]

Well the norm yould be just a real number, so printing that is easy. For Matrix manifold points I have overloaded the operators, for tangent vectors as well (on the most abstract layer). From my side the decision was also in favor of (what's called in Manopt.jl) TVectorE that also keeps its basepoint stored and can hence verify that operations (yes also those are overloaded) keep TVectors valid. I can understand your annoyance – but I would prefer to propose, that convert's might suffice, i.e. convert from an abastract array to an SnPoint or such.

In the end, converts or low level APIs are a solution to the same annoance you feel, just from two different side. And yes I can understand that annoyance.

Concerning the performance – for now Manopt.jl has also goals to educate and illustrate. Sure, performance is an issue still, but is not yet my main focus. I am also in favour of using views. For my sizes (512^2 point array of symmetric positive definite 3x3 matrices for example) the speed is enough (and at least 2x faster than Matlab/mex where I started out with). Finally, I also focused on documenting all formulae (instead of performance boosting) but of course being faster is always nice :)

I am also in favour of putting dispatch operations first (manifold always first, on tangent operators the base point second etc) – not because there's a style guide, but because I am used to think that way on manifolds. Despite being used to that order (and that it might take some time to get used to something different) I would be ok with changing that if a majority is for that.

Last but not least, I was wondering (and no offense meant) what we are aiming for with a comnmon framework. For best of cases that we can easily use one anothers manifolds. Then, the interfaces should be really fixed by conventions, but I think

• whether we do simple functions as base layer or conversions
• whether we aim for performance (in tricky-troubly isbits ;)) or validations (extended Vectors) and such
• whether we use views or not

should for best of cases not matter. IF we find such a common interface, that would be really nice, each implementation of a manifold focusing on a different goal, but still being interchangable. Do you think we can find such an interface? Which functions should it cover?

Otherwise I would also favoru to have a collection of manifolds in such a Package, whether that's a good idea, i don't know, maybe I just like to collect things.

[mateuszbaran]

Last but not least, I was wondering (and no offense meant) what we are aiming for with a comnmon framework. For best of cases that we can easily use one anothers manifolds. Then, the interfaces should be really fixed by conventions, but I think

Yes, that's what I was thinking: the possibility to mix manifolds and solvers from different packages. There are many good suggestions here so let me make a first draft: https://gist.github.com/mateuszbaran/6ac0115866562d09ea5e856216568a73

Notice that there are no restrictions on types of points and tangent vectors, as @kellertuer suggests, the interface should be transparent to such details.

I'm eagerly awaiting comments. I'm not particularly attached to names I've used there.

Defining non-mutating operations using mutating operations is very tempting but doing it in a way that is AD-friendly is challenging. For example, in log, either one or both points can be represented by duals, and if any is, the resulting vector should be represented by duals as well. In the mutating version taking care of this is delegated to the caller.

I think that it would be totally OK if not all manifolds implemented all these functions but it'd be useful to explicitly state which functions particular solvers rely upon.

Possibly a good solution would be to have a package with a low-level API, in the style retract(M::Manifold,x::AbstractArray) etc. Then, it would be easy to design a higher-level API with types that would call into this lower-level API.

The API I've sketched is pretty much level-independent actually, so maybe that's enough? Whether xs or vs are AbstractArrays or something else would be up to a particular implementation (existing in other packages like Manopt.jl or FunManifolds.jl).

Otherwise I would also favoru to have a collection of manifolds in such a Package, whether that's a good idea, i don't know, maybe I just like to collect things.

The problem with this is that such a collection would in some way endorse these implementations and I don't think we could agree what is the right way to implement, for example, the product manifold. Even myself, I see that both my approach and the Manopt.jl's have different goals and should just coexist.

Maybe such basic package shouldn't have any implementations of manifolds? As I see it, some reference (not even exported) for testing solvers would have some usefulness but maybe that's all.

Really? Reshapes and views used to suck a lot pre-1.0 but seem fine now. At least I've had good experiences on pretty large optimization problems.

That depends on the application and what is considered fine. Using heavily optimised linear algebra from StaticArrays can easily make an order of magnitude of difference. For example here: https://gist.github.com/mateuszbaran/9fa1c7ffcd3f9445c80503e80bd0dad8 is a benchmark that compares simple dot product using both reshapes and a vector of SArrays. I've even used UnsafeArrays to get the version with reshapes faster. Still, it's about 8-10 times slower than vector of SArrays. Using sized reshape is still 2x slower than vector of SArrays but maybe it's a better idea overall? I'll have to think about it at some point.

[antoine-levitt]

Cool! I commented on https://gist.github.com/mateuszbaran/9fa1c7ffcd3f9445c80503e80bd0dad8

Otherwise I would also favoru to have a collection of manifolds in such a Package, whether that's a good idea, i don't know, maybe I just like to collect things.

Me too :-) we could call this package ManifoldMuseum.jl and have it contain a number of manifolds with implementations of the basic operations on them.

[kellertuer]

Thanks Mateusz for the ideas. Actually Manopt.jl has retractions (but not yet fallbacks providing errors see here, but for example Stiefel has retractions (that I will rename to retract – I like that!).
My solvers have an option retr that defaults to exp where you can plug in retractions then. But let'S only start with solver interfaces after we've fixed the Manifold interface.

My problem with transparency (well or let's say something that is not yet clear to me) – if we introduce a type for points and vectors but still allow just arrays instead of these types – does that lead to unwanted side effects?

I love the idea of a ManifoldMuseum.jl– that sounds neat!

[antoine-levitt]

So I played around with these ideas a bit, and there are a few things that bother me:

• In principle the library should not depend on plain old Float64 vectors but should be able to work on more complicated structures: multidimensional, AD, GPU arrays, distributed, whatever. This does not play well with things like random_point(M) and friends that create things out of thin air...

• How much information should a manifold hold? In my implementation I've basically taken the approach that the manifold is an empty shell and just a way to dispatch on the correct manifold. But if you want things like random_point and dimension this can't really work...

• Julia doesn't have a great way to think about supporting both mutating and non-mutating variants (eg static arrays). I say let's forget about immutable arrays.

[kellertuer]

1. For the first point one could store the element-type in Manifold?
2. For me the manifold holds all information that is needed to for example do the things you call out-of-thin-air, for example Sphere(n) has a field dimension where n is stored. For more general cases, when you keep in mind that a discrete approximation to a(n infinite dimensional) shape manifold is a finite (graph) manifold – M also stores that graph. In short: M stores everything that is inherent to the manifold itself. My ones for example also store a small string with their name, power an array of its dimensions but for now not more.

For the third point I first have to understand the differences more – I think I am working on mutable ones but my manifold points are immutable? I think we should either decide for one or have a proper type hierarchy to dispatch to one or the other case.

[antoine-levitt]

For the first point one could store the element-type in Manifold?

That's not sufficient: eg points could be 2D arrays, or CuArrays, and then it'd be nice to keep everything of the same type.

My ones for example also store a small string with their name

That's typically the sort of things I'd do with dispatch (or not do at all)

For the third point I first have to understand the differences more – I think I am working on mutable ones but my manifold points are immutable? I think we should either decide for one or have a proper type hierarchy to dispatch to one or the other case.

No, immutables is for things like static arrays (https://github.com/JuliaArrays/StaticArrays.jl) where you can't do things like x .= 0 for instance.

[kellertuer]

Ah, then I am also for leaving out immutable arrays for now.

Concerning the type – hm, we could (ah you will not like that) introduce a Parametric Type?

[mateuszbaran]

• Good point. Storing element type in the manifold is an option but it will almost certainly fail in AD (I've tried a similar thing and ultimately gave up). Still, random_point can be useful in many scenarios.
• Dimension of the manifold can be stored as a type parameter, this IMHO makes sense.
• I'm OK with not caring about support for immutable arrays.

For the third point I first have to understand the differences more – I think I am working on mutable ones but my manifold points are immutable? I think we should either decide for one or have a proper type hierarchy to dispatch to one or the other case.

On top of what @antoine-levitt said, in my opinion most code should be agnostic of array type (that is, using lots of <:AbstractArray).

[kellertuer]

I use random_point quite often

Yes, I am in favour of <: AbstractArray coming from Matlab I was too lazy to do that for my types but I will refactor that in the future.

[antoine-levitt]

OK, here's a bare-bones version with Sphere, based on https://gist.github.com/mateuszbaran/6ac0115866562d09ea5e856216568a73 : https://github.com/JuliaNLSolvers/ManifoldMuseum. I'll add you two (and anyone interested) as committers, feel free to break everything

[kellertuer]

I don't understand your exp/log

• for me exp depends on a point x and a tangent v (no y)
• for me log depends on two points x,y (no v),

and for log I am also not sure what! should do – replace the value in x? That would be type instable for the approach with types I'd like to stick to. And – what if I really want to have y=exp(M,x,v)`? We could also continue over there further. Also, what to do with different types of noise (I am also doing image processing, there different types of noise are very important).

[antoine-levitt]

That's using the convention that mutated arguments go first (which are now second because we put manifold types first). exp!(y,x,v) computes exp(x,v) and puts the result in y. With the fallback exp(M::Manifold, x, v) = exp!(M,copy(x),x,v) you implement exp! and get exp for free.

Also, what to do with different types of noise (I am also doing image processing, there different types of noise are very important).

Where does that come in? In random_point? I'd do as for the choices of retractions and put that info in the manifold structure.

[kellertuer]

Thanks for the explanation, on exp

I am not that much a fan of putting that into M, but I see that that's a logical way to do (and I have no better conjecture instead). So in M it is.

[mateuszbaran]

OK, here's a bare-bones version with Sphere, based on https://gist.github.com/mateuszbaran/6ac0115866562d09ea5e856216568a73 : https://github.com/JuliaNLSolvers/ManifoldMuseum. I'll add you two (and anyone interested) as committers, feel free to break everything

Thanks, I'll soon start working on some tests and examples.

Also, what to do with different types of noise (I am also doing image processing, there different types of noise are very important).

Where does that come in? In random_point? I'd do as for the choices of retractions and put that info in the manifold structure.

OK, random_point doesn't look great to me now. We should probably add new methods to rand with our own VariateForm (like in Distributions.jl): ManifoldPointVariate, TangentVectorVariate or something.

[sethaxen]

I'm excited that this is happening. I'm pretty new to both Julia and topology/geometry, but I've been thinking about a general interface for manifolds and would be happy to contribute. My interest is statistics on Riemannian manifolds (and their products), such as matrix Lie groups.

A few additional desiderata:

1. The interface shouldn't force a specific set of coordinates for a manifold but instead permitted using different coordinate sets. I confess I have no idea how to do this.
2. An interface for Lie groups would be great and not too hard to add.
3. The interface should separate the concepts of manifold, coordinate, metric, and group. I know there tend to be commonly used metrics/coordinates for certain manifolds/groups, but for a general interface to be truly useful it would be nice to separate these.
4. It would also be great if the interface for product/power manifolds also handled some of their properties correctly (e.g. default product metric, product tangent bundle, etc), while permitting them to be overriden by the user or for special cases (e.g. IIRC the Lie Algebra of SE(n) is the tangent space of SE(n) at the identity but is not the same thing as the tangent space of its manifold Rn x SO(n) at the same point, due to the semidirect product)
5. As much should be handled at compile time to enable dispatch as possible. It would be ideal if the overhead from the manifolds interface was negligible.
6. Some basic framework for maps (e.g. pullbacks) between manifolds would be great.

[kellertuer]

Thanks for your input. I think to summarise the current state

1. the interface will be just a manifold type. coordinates could then be handled introducing a subtype. All manifolds I had in mind up to know (the sphere, SPDs,...) have coordinate free formulae, but you're right that the interface will not force certain coordinates
2. Yes, we already have an idea how to do that (similar to thing you gain for product you will gain certain things for Lie)
3. I am not yet sure how to do that, despite Traits for groups, we'll see.
4. Yes, that's covered – I did something similar here already. The hope is that you get all three “for free” if you define your own manifold
5. yes
6. not yet sure how to do that on the general level, but that's a good point to keep in mind.

[kellertuer]

And I think we should start to track these things in separate threads, since we already have a repository, see https://github.com/JuliaNLSolvers/ManifoldMuseum

[sethaxen]

Awesome! Will do!

[dkarrasch]

This is truly exciting progress here. I would like to bring a couple of aspects to your attention, just to see if there are chances that this can be handled within the targeted framework.

I'm interested in ODE integration, say on the sphere (ocean surface). As you can imagine, ocean surface velocity fields are provided only on the ocean surface. Viewing the sphere as embedded in R^3 and equipping it with manifold projections or whatever is really not helpful here. So one would want to work in (geographic) coordinates, that have the manifold constraint already built in. As long as you only want to integrate the ODE in coordinates there is no need to even start with a manifold framework code, you simply treat everything as if working in R^2 and that's it. However, in my field, we are often interested in stretching by the induced flow, which is essentially measured by pullbacks of metric tensors (or their inverses, interpreted as diffusion tensor fields). Funnily, diffusion tensor fields are, of course, not given in longitude-latitude coordinates, but in meters. Long story short, I wonder if handling coordinates/units and changes thereof, tensors and pullbacks is within the scope of this effort, or if that should be put on top of this framework elsewhere. See also Ferrite-FEM/Tensors.jl#107.

[antoine-levitt]

I don't have much experience with this but I think metric tensors can be handled in the current framework: just do a flat manifold with a varying inner product. We might want to add a function to get the metric tensor (which would fallback to I). That's one design objectives on my end: make the interface easy with lots of fallbacks, so that it's only a couple of lines to define a custom manifold that does what's of interest to your particular application.

Units are always a pain, I don't think we want to bother with Unitful in this package (support should be automatic as long as we don't specialize on Float64s).

Also, there are tons of things one could do with manifolds, but I think it'd be more productive to focus on functionality that is helpful to build non-trivial algorithms (ie optimisation, ODE solving) on top. But admittedly that might be too narrow-minded of me.

[dkarrasch]

I didn't necessarily mean to handle units explicitly. It's just that I looks like [cos(\theta)^2 0; 0 1] in lon-lat-coordinates. If you "image" I has units meter, than the other representation has different units. I don't work with explicit units either.

just do a flat manifold with a varying inner product

Yes, if there was a way to get these objects "automatically" knowing the coordinates of the current point and the corresponding coordinate system, that would be great. I guess my "feature wish" is to have some generic structure for coordinate systems, so that one would define different charts and transitions between them, and get transformation rules/pullbacks for free. This would somehow mean that any coordinate representation (vector/matrix/tensor) is associated to a chart. Not sure if that's adding too much overhead...

[sethaxen]

I guess my "feature wish" is to have some generic structure for coordinate systems, so that one would define different charts and transitions between them, and get transformation rules/pullbacks for free. This would somehow mean that any coordinate representation (vector/matrix/tensor) is associated to a chart.

I've been trying to think of a way to do this in my own work that doesn't add lots of overhead. Supposing we have an atlas for a manifold. To get a point from a set of coordinates, we need to define which chart should be applied to those coordinates to take them to a point. But if we perturb the coordinates by some amount, they may now be outside of the region covered by that chart. So before using the coordinates, we need to choose a chart, and then we need to check if the coordinates fit in that chart. It's a mess for efficiency.

[kellertuer]

For the atlas I would try to model the charts as (U,𝝋) i.e. domain and chart as nodes in a graph, where the edges at least have the intersection function (then of U,V) available and the chart transition. This way one could also look at all adjacent nodes of (U,𝝋) looking for a “neighboring” chart with a domain V that contains your perturbed point.

[sethaxen]

That's really clever. I'm still concerned about its efficiency though. I guess if we are stuck using an atlas with multiple charts, it might be the best we can do.

There's also the problem that we can't computationally represent a specific point on a manifold without selecting a set of coordinates and a chart with which to describe it. This isn't a problem for manifolds that have a canonical set of coordinates (e.g. matrices in matrix Lie groups).

Maybe the discussion of points, coordinates, and charts should move to its own issue.