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).