Ellipsoid manifold

[blegat]

Suppose I want to optimize over the set {x | x' * A' * A * x = 1}, e.g., the boundary of an ellipsoid. Is this manifold implemented somewhere ?

[mateuszbaran]

AFAIK it's not implemented yet. Probably people just reparameterize the objective as y=Ax and optimize on a sphere? I guess ellipsoid could be nice for convenience. Do you have in mind only the nice case of square, full-rank A or more general?

[kellertuer]

I would also probably use the reparametrization y=Ax, I think geodesics on ellipsoids are a bit nasty formula wise, and if A is not of full rank, it might be really not that nice to do, since then we basically have two sides of a flat ellipsoidal shape in a plane (in the 3D case with one eigenvalue 0)

[blegat]

If I was minimizing f(y) where y in A S where S is the sphere then I could just min f(A x) where x in S. Here, it's more like f(x) where A x in S. Say x is dim n, S is dim m - 1 and A is n x m and rank(A) is m. The set is more like a semi-ellipsoid embedding in R^n which looks like a cylinder that infinite in n - m directions. I could do the SVD on A and then reformulate but the SVD might be more expensive than the manifold optimization on the semi-ellipsoid (in my case, the objective is `norm(B x)so it's some kind of generalized Rayleigh quotient and the proof that the SOCP found by manifold optimization will be global minimum seems to generalize easily to the semi-ellipsoid case). The formulas for the projections onto the manifold and the projections onto the tangent space weren't too hard and I would make gradient descent work for simple examples with aConstantStepsize(0.1)and gradient descent. Then I tried without setting the stepsize but then I need to define the injectivity radius. I triedπ/2withProjectionRetraction` but it didn't seem to converge but I haven't investigated too much.

[kellertuer]

Why do you need the infectivity radius? We only use that as a default in a few places, if you do not provide (e.g. the constant in ConstantStepsize or the initial step size in Armijo, I think).

Constant Step size and gradient descent is usually not the best of ideas, since there is no convergence guarantee. We still use that as a default, but only since Armijo requires more things available.

If you do not set the (constant) stepwise, in ordre to get a good guess it uses the infectivity radius (we can not just roll a dice here). And if instead of exp we use a retraction, the convergence is (though only a bit) less likely.

[blegat]

Why do you need the infectivity radius?

Yes, it's in the Armijo line search. By the way, while doing this I thought it would be nice to add a Manopt tutorial showing the minimal interface needed in order to optimize on it. I basically did a trial and error following the MethodError's

[kellertuer]

to get started, there is this tutorial https://manoptjl.org/stable/tutorials/Optimize!/ can you be a bit more precise what you are missing? I am happy to add more details there or more tutorials.

For the error in Armijo – I checked and actually we do not need it for the initial step but for the safeguard (stop_when_stepsize_exceeds) which uses max_stepsize see https://manoptjl.org/stable/plans/stepsize/#Manopt.ArmijoLinesearch (and yes there a link is broken, will fix that soon-ish).

So let me know which method errors you run in and what tutorial you would have wished for. I am happy to add that.

[mateuszbaran]

@kellertuer there is currently no tutorial on implementing new manifolds for particular solvers, and basically the only way to figure out if manifold X can be used with solver Y is "does it throw errors when you try?"

[kellertuer]

Hm, but I am unsure what to do about this. It would mean we come up with some place somewhere (no clue where) that says:

“If you want to use solver Y on your now manifold, you at least need the following functions?”

And for 95% this is retraction, inverse retraction, norm. And the rest depends on single step sizes or stopping criteria.

So one could write that for every “ingredient” of a solver, but I am not sure where to put such a thing.

[mateuszbaran]

Maybe let's put it in documentation page of each solver?

[kellertuer]

We could do a “Technical Details” section for each solver basically, but that would not help with Armijo linesearch, which seems to have been the reason here.

[blegat]

What is missing is a tutorial doing struct CustomManifold <: ManifoldsBase.AbstractManifold and showing which functions need to be implemented for, e.g., gradient descent to work. If it depends on the solver, I guess just the minimal interface to start with. If you look at the existing implementations of manifolds in this package, it's also mixed it things that are there to speed up but it's best to first present the minimal interface to make it work. The stacktrace is

ERROR: MethodError: no method matching injectivity_radius(::Macaulay.Ellipsoid{Float64, LinearAlgebra.Adjoint{Float64, Matrix{Float64}}})

Closest candidates are:
  injectivity_radius(::AbstractManifold, ::AbstractRetractionMethod)
   @ ManifoldsBase ~/.julia/packages/ManifoldsBase/vKGuo/src/ManifoldsBase.jl:507
  injectivity_radius(::AbstractManifold, ::Any)
   @ ManifoldsBase ~/.julia/packages/ManifoldsBase/vKGuo/src/ManifoldsBase.jl:503
  injectivity_radius(::AbstractManifold, ::Any, ::AbstractRetractionMethod)
   @ ManifoldsBase ~/.julia/packages/ManifoldsBase/vKGuo/src/ManifoldsBase.jl:504
  ...

Stacktrace:
 [1] max_stepsize(M::Macaulay.Ellipsoid{Float64, LinearAlgebra.Adjoint{Float64, Matrix{Float64}}})
   @ Manopt ~/.julia/packages/Manopt/A4a09/src/plans/stepsize.jl:39
 [2] default_stepsize(M::Macaulay.Ellipsoid{Float64, LinearAlgebra.Adjoint{Float64, Matrix{Float64}}}, ::Type{Manopt.GradientDescentState}; retraction_method::ProjectionRetraction)
   @ Manopt ~/.julia/packages/Manopt/A4a09/src/solvers/gradient_descent.jl:89

[kellertuer]

As I wrote above, that is very solver-dependent. For the default gradient-descent-with-constant-stepsize you only need a retraction (and default_retraction_method) and for the stopping criterion maybe the norm of tangent vectors (norm(M, p, X)).

And what you are missing is either an upper bound on the maximal step size (for example anything beyond π does not make much sense as a step size on the sphere) or the injectivitiy radius. But this is independent of gradient descent – any solver with a step computation would require that as soon as you want to use Armijo there.

That is why I am not yet so sure what to write in such a tutorial, since it is really super dependent on the solver and its single ingredients (beyond the tutorial https://manoptjl.org/stable/tutorials/Optimize!/ to get started with optimization in general and https://juliamanifolds.github.io/ManifoldsBase.jl/stable/tutorials/implement-a-manifold/ on how to get started with manifolds in general).

If you look at the existing implementations of manifolds in this package, it's also mixed it things that are there to speed up but it's best to first present the minimal interface to make it work.

Again https://juliamanifolds.github.io/ManifoldsBase.jl/stable/tutorials/implement-a-manifold/ is how to get started with implementing a manifold.

“A simple interface to make it work” does not exist for solvers

• some (like CG) need a vector transport as well
• some need an inverse retraction as well
• some need also themselves the infectivity_radius (trust region for example)

so independent of which tutorial I would write, you could always find a next solcver for which the tutorial would not be enough 😕 I am not even sure where to put such a tutorial. It does not belong into Manopt.jl, since it is not about implementing a solver, it also would not so much belong here, since I would discuss again details you need to implement here to be used in Manopt.jl. So any tips what you would prefer would help me a lot. I am just not sure what to write nor where to put that.

[blegat]

I think Manopt would make more sense since it's at the intersection between ManifoldsBase and Manopt. I think a tutorial showing it for one solver is fine. Once the user has something working for one solver, it's easier to try to get it working for another solver. But getting it to work for the first solver could be discouraging for the user since he does not know how far he is and whether is at 10% of the MethodErrors or at 90% ^^ Also, since it's a tutorial and not a manual, it's fine if it just shows it for one solver. A tutorial is never meant to do "everything", it's just one sample from the space of possible things you can do ^^

[kellertuer]

I will think about it, since I already linked the ManifoldsBase tutorial, it would maybe be even a good idea to start with the same manifold idea from there (Sphere with a radius) and extend that example to explain which things are needed.

Might take a bit until I find time to write this – we are in the middle of a semester here in Norway, so I am busy with teaching as well – but I will think about that then and try to find a reasonable way of presentation.

[kellertuer]

I gave the tutorial idea a few thoughts yesterday evening, and my main problem is – I think – that it will be really just for one solver then, because most solvers have a few additional functions they would need from the API ManifoldsBase.jl defines – and which ones they do need is quite individual – besides a retraction and an inverse retraction, those two every solver probably needs. I still have to think how to turn something for this into a good tutorial – for now I do not have a good idea; something that works for exactly one solver is not a good tutorial how to approach that for the others then.

[mateuszbaran]

What about adding a function to Manopt that takes solver as argument and lists what needs to be implemented? Something like required_manifold_features(quasi_Newton) that returns a string like:

* retract!(M::MyManifold, q, p, X, t::Real)
* default_vector_transport_method(M::MyManifold)
* vector_transport_to!(M::MyManifold, Y, p, X, q, ::MyMethod) # MyMethod is returned by default_vector_transport_method(M::MyManifold)
* inner(M::MyManifold, p, X, Y)

[mateuszbaran]

It could potentially have more details but this would be the general idea.

[kellertuer]

That sounds reasonable, but maybe it is more accessible when doing a section “Technical Details” at the end of every solver documentation page?

[mateuszbaran]

That was my idea from two weeks ago that you didn't particularly like: https://github.com/JuliaManifolds/Manifolds.jl/issues/658#issuecomment-1759333881 :smile:

[kellertuer]

Yes, but your new comment would also not help with Armijo ;) Since we could only refer in the solver to “this solver does line search” – and maybe do a technical details on the line search documentation? Listing that for every line search might be a bit crowded, but a way to go.

[mateuszbaran]

A method like that could also call a similar method for the default line search and concatenate the strings. And point out that there are different line searches with different requirements. Anyway, I think this isn't something that can be centralized and a tutorial should just point to places where one needs to look for solver-specific information.

[kellertuer]

Sure there would be a general remark like “This solver uses a line search – check page X for further details”. That would be fine with me. The main decision is whether to write that as a function (and what should that function return?) or sections in the docs. But together with the docs – sure one can write then a tutorial.

[mateuszbaran]

I think docs would be easier to prepare than a function. A function could automatically aggregate information but that may not be necessary when a tutorial would explain how to do it.

[kellertuer]

Both would require about the same amount of manual collection of information I think and then the documentation section is the one more accessible to users. I think I will go for the docs section then – when I find time :)

[kellertuer]

I did the “Technical Details” section on every solver and a tutorial to implement a manifold (going along what a solver needs) – so I will close this for now, if you feel we should work on an ellipsoid manifold or we can help there in any way, feel free to reopen this :)

Of course new manifolds are always welcome as a PR.