search

JuliaManifolds / Manopt.jl

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

Riemannian Levenberg-Marquardt initial residual and jacobian as kwargs #268

Closed Affie closed 1 year ago

Affie commented 1 year ago

Hi, I wanted to test using Levenberg Marquardt in https://github.com/JuliaRobotics/IncrementalInference.jl/pull/1724 but came across 2 issues that I fixed in this PR.

  1. initial_residual_values = similar(p, get_objective(nlso).num_components) Didn't work for me on SE(n) with p an ArrayPartion of (SVector,SMatrix)

  2. initial_jacF = similar(p, get_objective(nlso).num_components, manifold_dimension(M) I wanted to use a sparse Jacobian so I had to supply the initial jacobian myself.

I'm open to suggestions of other ways to do it, but thought I'd rather open a PR than an issue.

PS. I almost have Riemannian Levenberg-Marquardt working as one of our solvers, but I'm not getting the correct results yet, and still need to figure out some performance issues. But thanks for another great addition to julia manifolds.

kellertuer commented 1 year ago

To me this looks good. However Mateusz is the expert here, he wrote the algorithm.

  1. He also designed the similar function, I would probably have written that as a [copy(M, p) for _ in 1:get_objective(nlso).num_components] That said, I would like to know the example that failed and/or the error message? It seems our default here should be improved so it also works “out of the box“ for SE(n)

  2. sure moving both to keywords is fine, just the naming as I wrote could maybe be improved.

Ah, and thanks for your kind words :) While I do care a lot about Manifolds.jl – is Manopt.jl my start of all this and still the package I care about most. So I am very happy when it is useful!

kellertuer commented 1 year ago

Ah, and I noticed: Of course new keyword arguments have to be documented in the doctoring in the beginning. We order them in alphabetical order in line 31-37

mateuszbaran commented 1 year ago

Nice, I think this is a good change.

Regarding sparsity, I would need to check if you actually properly exploit it. Just using a block array may not be enough, and I haven't tested this implementation with sparse Jacobians. You can take a look here for some comments about sparse variants of LM: https://arxiv.org/pdf/2206.05188.pdf . Do you have an example I could run?

codecov[bot] commented 1 year ago

Codecov Report

Merging #268 (3bfb3f0) into master (ec293ad) will not change coverage. The diff coverage is n/a.

@@           Coverage Diff           @@
##           master     #268   +/-   ##
=======================================
  Coverage   99.70%   99.70%           
=======================================
  Files          73       73           
  Lines        6465     6465           
=======================================
  Hits         6446     6446           
  Misses         19       19
Impacted Files Coverage Δ
src/solvers/LevenbergMarquardt.jl 100.00% <ø> (ø)

:mega: We’re building smart automated test selection to slash your CI/CD build times. Learn more

Affie commented 1 year ago

Thanks for the article, I'm sure I'm not properly exploiting sparsity yet, but at least it did improve performance quite a bit. It's a bit hard to make an MWE for you to look at as it integrates with IncrementalInference and RoME. I'll try and make something similar but in the meantime, here is my attempt so far: https://github.com/JuliaRobotics/IncrementalInference.jl/pull/1724/files#diff-ce2c498913efe786709e2e5d9cdb6a01adebe4a344429cfb02f8dcc4d77ea60d

EDIT: To run a simple example on that branch (requires this pr though)

#]add IncrementalInference#23Q2/enh/manopt
#]add RoME#master
using RoME, IncrementalInference
fg = generateGraph_Hexagonal(; graphinit=false, landmark=false)
v,result = IIF.solve_RLM_sparse(fg)

This generates and solves a small factor graph of a robot driving in a hexagonal, which results in: image

Affie commented 1 year ago

To me this looks good. However Mateusz is the expert here, he wrote the algorithm.

  1. He also designed the similar function, I would probably have written that as a [copy(M, p) for _ in 1:get_objective(nlso).num_components] That said, I would like to know the example that failed and/or the error message? It seems our default here should be improved so it also works “out of the box“ for SE(n)

Yes, I think the default is the issue in my case. I have a vector of points on SE(n) and get this error without the change.

ERROR: MethodError: Cannot `convert` an object of type Float64 to an object of type ArrayPartition{Float64, Tuple{StaticArraysCore.SVector{2, Float64}, StaticArraysCore.SMatrix{2, 2, Float64, 4}}}

Closest candidates are:
  convert(::Type{V}, ::ManifoldsBase.ValidationMPoint{V}) where V<:AbstractArray
   @ ManifoldsBase ~/.julia/packages/ManifoldsBase/74WVY/src/ValidationManifold.jl:84
  convert(::Type{T}, ::Manifolds.HyperboloidPoint, ::Manifolds.HyperboloidTVector) where T<:(AbstractVector)
   @ Manifolds ~/.julia/packages/Manifolds/9j170/src/manifolds/HyperbolicHyperboloid.jl:59
  convert(::Type{T}, ::Factorization) where T<:AbstractArray
   @ LinearAlgebra ~/.julia/juliaup/julia-1.9.1+0.x64.linux.gnu/share/julia/stdlib/v1.9/LinearAlgebra/src/factorization.jl:59
  ...

Stacktrace:
  [1] setindex!(A::Vector{ArrayPartition{Float64, Tuple{StaticArraysCore.SVector{2, Float64}, StaticArraysCore.SMatrix{2, 2, Float64, 4}}}}, x::Float64, i1::Int64)
    @ Base ./array.jl:969
  [2] macro expansion
    @ ~/.julia/packages/StaticArrays/O6dgq/src/broadcast.jl:160 [inlined]
  [3] _broadcast!
    @ ~/.julia/packages/StaticArrays/O6dgq/src/broadcast.jl:148 [inlined]
  [4] _copyto!
    @ ~/.julia/packages/StaticArrays/O6dgq/src/broadcast.jl:72 [inlined]
  [5] copyto!
    @ ~/.julia/packages/StaticArrays/O6dgq/src/broadcast.jl:65 [inlined]
  [6] materialize!
    @ ./broadcast.jl:884 [inlined]
  [7] materialize!(dest::Vector{ArrayPartition{Float64, Tuple{StaticArraysCore.SVector{2, Float64}, StaticArraysCore.SMatrix{2, 2, Float64, 4}}}}, bc::Base.Broadcast.Broadcasted{StaticArraysCore.StaticArrayStyle{1}, Nothing, typeof(identity), Tuple{StaticArraysCore.SVector{21, Float64}}})
    @ Base.Broadcast ./broadcast.jl:881
  [8] (::IncrementalInference.CostF_RLM!{ArrayPartition{Float64, Tuple{StaticArraysCore.SVector{2, Float64}, StaticArraysCore.SMatrix{2, 2, Float64, 4}}}})(M::IncrementalInference.NPowerManifold{ℝ, Manifolds.GroupManifold{ℝ, Manifolds.ProductManifold{ℝ, Tuple{TranslationGroup{Tuple{2}, ℝ}, SpecialOrthogonal{2}}}, Manifolds.SemidirectProductOperation{Manifolds.RotationAction{TranslationGroup{Tuple{2}, ℝ}, SpecialOrthogonal{2}, Manifolds.LeftAction}}}}, x::Vector{ArrayPartition{Float64, Tuple{StaticArraysCore.SVector{2, Float64}, StaticArraysCore.SMatrix{2, 2, Float64, 4}}}}, p::Vector{ArrayPartition{Float64, Tuple{StaticArraysCore.SVector{2, Float64}, StaticArraysCore.SMatrix{2, 2, Float64, 4}}}})
    @ IncrementalInference ~/.julia/packages/IncrementalInference/h3MMp/src/ParametricManoptDev.jl:88
  [9] initialize_solver!
    @ ~/.julia/packages/Manopt/bYnCJ/src/solvers/LevenbergMarquardt.jl:204 [inlined]
 [10] initialize_solver!(amp::Manopt.DefaultManoptProblem{IncrementalInference.NPowerManifold{ℝ, Manifolds.GroupManifold{ℝ, Manifolds.ProductManifold{ℝ, Tuple{TranslationGroup{Tuple{2}, ℝ}, SpecialOrthogonal{2}}}, Manifolds.SemidirectProductOperation{Manifolds.RotationAction{TranslationGroup{Tuple{2}, ℝ}, SpecialOrthogonal{2}, Manifolds.LeftAction}}}}, Manopt.NonlinearLeastSquaresObjective{Manopt.InplaceEvaluation, IncrementalInference.CostF_RLM!{ArrayPartition{Float64, Tuple{StaticArraysCore.SVector{2, Float64}, StaticArraysCore.SMatrix{2, 2, Float64, 4}}}}, IncrementalInference.JacF_RLM!{IncrementalInference.CostF_RLM!{ArrayPartition{Float64, Tuple{StaticArraysCore.SVector{2, Float64}, StaticArraysCore.SMatrix{2, 2, Float64, 4}}}}, Vector{ArrayPartition{Float64, Tuple{StaticArraysCore.SVector{2, Float64}, StaticArraysCore.SMatrix{2, 2, Float64, 4}}}}}, ManifoldsBase.DefaultOrthonormalBasis{ℝ, ManifoldsBase.TangentSpaceType}}}, dss::Manopt.DebugSolverState{Manopt.LevenbergMarquardtState{Vector{ArrayPartition{Float64, Tuple{StaticArraysCore.SVector{2, Float64}, StaticArraysCore.SMatrix{2, 2, Float64, 4}}}}, Manopt.StopWhenAny{Tuple{Manopt.StopAfterIteration, Manopt.StopWhenGradientNormLess, Manopt.StopWhenStepsizeLess}}, ManifoldsBase.ExponentialRetraction, Vector{ArrayPartition{Float64, Tuple{StaticArraysCore.SVector{2, Float64}, StaticArraysCore.SMatrix{2, 2, Float64, 4}}}}, Matrix{ArrayPartition{Float64, Tuple{StaticArraysCore.SVector{2, Float64}, StaticArraysCore.SMatrix{2, 2, Float64, 4}}}}, Vector{ArrayPartition{Float64, Tuple{StaticArraysCore.SVector{2, Float64}, StaticArraysCore.SMatrix{2, 2, Float64, 4}}}}, Float64}})
    @ Manopt ~/.julia/packages/Manopt/bYnCJ/src/solvers/debug_solver.jl:8
 [11] solve!
    @ ~/.julia/packages/Manopt/bYnCJ/src/solvers/solver.jl:118 [inlined]
 [12] LevenbergMarquardt!(M::IncrementalInference.NPowerManifold{ℝ, Manifolds.GroupManifold{ℝ, Manifolds.ProductManifold{ℝ, Tuple{TranslationGroup{Tuple{2}, ℝ}, SpecialOrthogonal{2}}}, Manifolds.SemidirectProductOperation{Manifolds.RotationAction{TranslationGroup{Tuple{2}, ℝ}, SpecialOrthogonal{2}, Manifolds.LeftAction}}}}, nlso::Manopt.NonlinearLeastSquaresObjective{Manopt.InplaceEvaluation, IncrementalInference.CostF_RLM!{ArrayPartition{Float64, Tuple{StaticArraysCore.SVector{2, Float64}, StaticArraysCore.SMatrix{2, 2, Float64, 4}}}}, IncrementalInference.JacF_RLM!{IncrementalInference.CostF_RLM!{ArrayPartition{Float64, Tuple{StaticArraysCore.SVector{2, Float64}, StaticArraysCore.SMatrix{2, 2, Float64, 4}}}}, Vector{ArrayPartition{Float64, Tuple{StaticArraysCore.SVector{2, Float64}, StaticArraysCore.SMatrix{2, 2, Float64, 4}}}}}, ManifoldsBase.DefaultOrthonormalBasis{ℝ, ManifoldsBase.TangentSpaceType}}, p::Vector{ArrayPartition{Float64, Tuple{StaticArraysCore.SVector{2, Float64}, StaticArraysCore.SMatrix{2, 2, Float64, 4}}}}; retraction_method::ManifoldsBase.ExponentialRetraction, stopping_criterion::Manopt.StopWhenAny{Tuple{Manopt.StopAfterIteration, Manopt.StopWhenGradientNormLess, Manopt.StopWhenStepsizeLess}}, debug::Vector{Manopt.DebugWarnIfCostIncreases}, expect_zero_residual::Bool, kwargs::Base.Pairs{Symbol, Any, Tuple{Symbol, Symbol, Symbol}, NamedTuple{(:evaluation, :initial_residual_values, :initial_jacobian_f), Tuple{Manopt.InplaceEvaluation, Vector{Float64}, SparseArrays.SparseMatrixCSC{Float64, Int64}}}})
    @ Manopt ~/.julia/packages/Manopt/bYnCJ/src/solvers/LevenbergMarquardt.jl:185
 [13] #LevenbergMarquardt#647
    @ ~/.julia/packages/Manopt/bYnCJ/src/solvers/LevenbergMarquardt.jl:114 [inlined]
 [14] LevenbergMarquardt
    @ ~/.julia/packages/Manopt/bYnCJ/src/solvers/LevenbergMarquardt.jl:110 [inlined]
 [15] #LevenbergMarquardt#646
    @ ~/.julia/packages/Manopt/bYnCJ/src/solvers/LevenbergMarquardt.jl:108 [inlined]
mateuszbaran commented 1 year ago

Thanks for the example. I've taken a look at performance. There seem to be no major issues that can be solved by small changes. It looks like it could be made 2x-3x faster by fixing small issues. For example calcfacs = CalcFactorManopt.(facs, Ref(varIntLabel)) is not ideal, there seems to be a type instability inside so I suggest changing it to map.

Do the results look correct in this example though? Or, what is wrong there? I think correctness should be solved first.

kellertuer commented 1 year ago

...but could we improve our default that it at least does not error when initialising LM?

Affie commented 1 year ago

Thanks for the example. I've taken a look at performance. There seem to be no major issues that can be solved by small changes. It looks like it could be made 2x-3x faster by fixing small issues. For example calcfacs = CalcFactorManopt.(facs, Ref(varIntLabel)) is not ideal, there seems to be a type instability inside so I suggest changing it to map.

Thanks, I still haven't figured out type stability vs compile time quite yet. Also, still have to figure out how to improve type stability on the CalcFactorManopt struct. I did see it was dynamically dispatching in places. The other performance issue we are currently working on is changing over to support StaticArrays (mostly done) we still have to get rid of all the allocations in the cost function.

Do the results look correct in this example though? Or, what is wrong there? I think correctness should be solved first.

The hex result is correct. I do not get the correct results in a larger graph on SE(3) yet. I'm still learning Manopt and suspect it might not have converged yet as it does look close.

kellertuer commented 1 year ago

The hex result is correct. I do not get the correct results in a larger graph on SE(3) yet. I'm still learning Manopt and suspect it might not have converged yet as it does look close.

Oh I was just referring to the initialisation, not to your solver run. You can easily look at some debug if that helps (for example whether your gradient norm is small). Correct I referred to as “the default does not error on SE(3)” ;)

Affie commented 1 year ago

...but could we improve our default that it at least does not error when initialising LM?

@kellertuer, I don't know what defualt initialization will work here, I think it needs the element type, I just used initial_residual_values = zeros(num_components) but that won't work for ForwardDiff. I couldn't find an eltype like function that works in all situations. for example initial_residual_values = zeros(eltype(p), num_components)

[copy(M, p) for _ in 1:get_objective(nlso).num_components]

Doesn't work as it's not on the manifold but Rn. If I'm not misunderstanding RLM completely.

What else is needed for this PR, is it just better default values for initial_residual_values and initial_jacobian_f?

kellertuer commented 1 year ago

...but could we improve our default that it at least does not error when initialising LM?

@kellertuer, I don't know what defualt initialization will work here, I think it needs the element type, I just used initial_residual_values = zeros(num_components) but that won't work for ForwardDiff. I couldn't find an eltype like function that works in all situations. for example initial_residual_values = zeros(eltype(p), num_components)

[copy(M, p) for _ in 1:get_objective(nlso).num_components]

Doesn't work as it's not on the manifold but Rn. If I'm not misunderstanding RLM completely.

That's why I used M – the manifold one is on and not Rn, but it was also just a guess, since I am still trying to understand which values we need.

So here is why I am confused.

1) for initial_residual_values the name indicates it might just be numbers? But since p is a point on the manifold your two variants

2) The second variable I do not understand yet either. initial_jacobian_f for me indicates it might be a matrix? but the default is similar(p, get_objective(nlso).num_components, manifold_dimension(M) so it might even be coefficients with respect to a matrix (due to the manifold dimension and not some representation size used). If that is the case this could be resolved with solving the same question as for 1.

I would have thought 2) might also be a collection of points or tangent vectors, hence my copy idea, but maybe I misread it and both are really just a vector and a matrix. Then we would need a good way to find the eltype (let's call it E and then we could just use zero(E, [...size information from above...]) for both

So @mateuszbaran what could we best do here to get the right number type even fof ArrayPartition and SVDMPoint points on manifolds?

mateuszbaran commented 1 year ago

So @mateuszbaran what could we best do here to get the right number type even fof ArrayPartition and SVDMPoint points on manifolds?

I would propose merging this PR as-is and make an issue about improving default for initial residuals for later.