JuliaManifolds / Manopt.jl

🏔️Manopt. jl – Optimization on Manifolds in Julia
http://manoptjl.org
Other
321 stars 40 forks source link

Nelder-Mead seems to be much worse than in Optim.jl #179

Closed mateuszbaran closed 1 year ago

mateuszbaran commented 1 year ago

See the following code, optiomization of the Rosenbrock function using Optim.jl and Manopt.jl:

using Optim, Manopt
using Manifolds

f(x) = (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
f_manopt(::Euclidean, x) = (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2

optimize(f, [0.0, 0.0], Optim.NelderMead())

M = Euclidean(2)
opts = Manopt.NelderMead(M, f_manopt; return_options=true)

when executed I get:

julia> using Optim, Manopt

julia> using Manifolds

julia> f(x) = (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
f (generic function with 1 method)

julia> f_manopt(::Euclidean, x) = (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
f_manopt (generic function with 1 method)

julia> optimize(f, [0.0, 0.0], Optim.NelderMead())
 * Status: success

 * Candidate solution
    Final objective value:     3.525527e-09

 * Found with
    Algorithm:     Nelder-Mead

 * Convergence measures
    √(Σ(yᵢ-ȳ)²)/n ≤ 1.0e-08

 * Work counters
    Seconds run:   0  (vs limit Inf)
    Iterations:    60
    f(x) calls:    117

julia> M = Euclidean(2)
Euclidean(2; field = ℝ)

julia> opts = Manopt.NelderMead(M, f_manopt; return_options=true)
NelderMeadOptions{Vector{Float64}, StopAfterIteration, Float64, Float64, Float64, Float64, ExponentialRetraction, LogarithmicInverseRetraction}([[0.12499192873081919, 0.009135702816670005], [0.12499192873081919, 0.009135702816670005], [0.12499192873081919, 0.009135702816670005]], StopAfterIteration(200000, "The algorithm reached its maximal number of iterations (200000).\n"), 1.0, 2.0, 0.5, 0.5, [0.12499192873081919, 0.009135702816670005], [0.7698476042280331, 0.7698476042280331, 0.7698476042280331], ExponentialRetraction(), LogarithmicInverseRetraction())

julia> f(opts.x)
0.7698476042280331

So Manopt calls objective much more times than Optim and still returns a much worse solution. Could that be a bug in Manopt?

kellertuer commented 1 year ago

Hm Nelder Mead. That is a strange algorithm and it might be hard to get the right default values. If I read that correctly we use the 4 default values from the paper of Nelder&Mead, while the default in Optim is more advanced and changes during the iterations (that is not possible or would require some implementation).

So I would only speak of a bug if with completely the same parameters the behaviour would be tremendously different.

On the other hand … sure there can be a bug, I have not used the algorithm much, but all tests that I wrote show that the function value decreases.

By the way, if you just need opts.x you can set return_options=false and the solver will directly return x.

mateuszbaran commented 1 year ago

I've checked it, and the implementation seems incorrect. Moreover, the initial simplex can be constructed much better if we have an initial point and a basis at that point. I will make a PR soon with improvements.

mateuszbaran commented 1 year ago

My only issue is how to make an API for specifying a point and a basis given that by default NelderMead already has a parameter for population. Do you have a preference?

kellertuer commented 1 year ago

The only disadvantage with such an PR now is, that NelderMead is the one algorithms that is already completely (with running tests) done in the new framework, so a PR on that would lead to me merging for two hours that into the new one. So while these ideas are all nice it is unfortunate that these ideas hit the same time as I am trying to cope with your first idea and polishing the code overall in the Rework PR.

Besides that, we can just type the popultation, i.e. AbstractArray is the current population and Tuple{Point, Basis} Is the other – since I think the size of the population is always n+1?

mateuszbaran commented 1 year ago

OK, I can wait with that PR.

Yes, the size of population is always n+1. How about, instead of a Tuple, making a new type NelderMeadSimplex that can be initialized from a point and a basis? Something like NelderMeadSimplex(M, x0, DefaultOrthonormalBasis()). This way we could add more customization there.

kellertuer commented 1 year ago

...and the current variant is a constructor NelderMeadSimplex(M, X)? Sure sounds reasonable.

mateuszbaran commented 1 year ago

The current would be just NelderMeadSimplex(M) (no X).

kellertuer commented 1 year ago

Ah, I thought one could make the X an optional keyword, which gets filled with random points (like the current default) and if the user specifies X, then it is the identity? Hm sure but in chat case he could specify the population also directly. Then yours sounds more reasonable.

mateuszbaran commented 1 year ago

I will figure out the details of directly specifying a population later, maybe it will be nicer to also handle it by NelderMeadSimplex.