Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

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

Merged
merged 5 commits into from
Jun 15, 2023

Conversation

Affie
Copy link
Contributor

@Affie Affie commented Jun 12, 2023

Hi, I wanted to test using Levenberg Marquardt in JuliaRobotics/IncrementalInference.jl#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
Copy link
Member

kellertuer commented Jun 12, 2023

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 kellertuer requested a review from mateuszbaran June 12, 2023 19:42
@kellertuer
Copy link
Member

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

Affie and others added 2 commits June 13, 2023 09:49
@mateuszbaran
Copy link
Member

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
Copy link

codecov bot commented Jun 13, 2023

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% <ø> (ø)

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

@Affie
Copy link
Contributor Author

Affie commented Jun 13, 2023

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
Copy link
Contributor Author

Affie commented Jun 13, 2023

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
Copy link
Member

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
Copy link
Member

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

@Affie
Copy link
Contributor Author

Affie commented Jun 13, 2023

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
Copy link
Member

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
Copy link
Contributor Author

Affie commented Jun 15, 2023

...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
Copy link
Member

...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
  • initial_residual_values = zeros(num_components) is a vector of numbers
  • initial_residual_values = zeros(eltype(p), num_components) is a vector of numbers but matching the type the entries of p have; which might be complicated for array partition or other non-just-array-types.
  1. 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
Copy link
Member

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.

Copy link
Member

@kellertuer kellertuer left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

THat's of course also fine with me.

@kellertuer kellertuer merged commit 9adae44 into JuliaManifolds:master Jun 15, 2023
@Affie Affie deleted the rlm_sparse_jac branch June 15, 2023 17:53
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants