Bundle adjustment model or CaNNOLeS optimizer not working?

[ojwoodford]

Hi there. All these packages look interesting, and I wanted to try them out on bundle adjustment. I was expecting the following code to work:

using BundleAdjustmentModels, CaNNOLeS
model = BundleAdjustmentModel("problem-16-22106")
stats = cannoles(model)

but I get the following error:

ERROR: MethodError: no method matching hess_structure_residual!(::BundleAdjustmentModel{Float64, Vector{Float64}}, ::Vector{Int64}, ::Vector{Int64})
Stacktrace:
 [1] hess_structure_residual(nls::BundleAdjustmentModel{Float64, Vector{Float64}})
   @ NLPModels ~/.julia/packages/NLPModels/ZITAG/src/nls/api.jl:286
 [2] solve!(solver::CaNNOLeSSolver{Float64, Vector{Float64}}, nls::BundleAdjustmentModel{Float64, Vector{Float64}}, stats::SolverCore.GenericExecutionStats{Float64, Vector{Float64}, Vector{Float64}, Any}; x::Vector{Float64}, λ::Vector{Float64}, method::Symbol, linsolve::Symbol, max_f::Int64, max_time::Float64, max_inner::Int64, ϵtol::Float64, verbose::Int64, check_small_residual::Bool, always_accept_extrapolation::Bool, δdec::Float64)
   @ CaNNOLeS ~/.julia/packages/CaNNOLeS/EgvjT/src/CaNNOLeS.jl:175
 [3] solve!(solver::CaNNOLeSSolver{Float64, Vector{Float64}}, nls::BundleAdjustmentModel{Float64, Vector{Float64}}, stats::SolverCore.GenericExecutionStats{Float64, Vector{Float64}, Vector{Float64}, Any})
   @ CaNNOLeS ~/.julia/packages/CaNNOLeS/EgvjT/src/CaNNOLeS.jl:114
 [4] solve!(solver::CaNNOLeSSolver{Float64, Vector{Float64}}, model::BundleAdjustmentModel{Float64, Vector{Float64}}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ SolverCore ~/.julia/packages/SolverCore/kzlVQ/src/solver.jl:40
 [5] solve!(solver::CaNNOLeSSolver{Float64, Vector{Float64}}, model::BundleAdjustmentModel{Float64, Vector{Float64}})
   @ SolverCore ~/.julia/packages/SolverCore/kzlVQ/src/solver.jl:38
 [6] cannoles(nls::BundleAdjustmentModel{Float64, Vector{Float64}}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ CaNNOLeS ~/.julia/packages/CaNNOLeS/EgvjT/src/CaNNOLeS.jl:111
 [7] cannoles(nls::BundleAdjustmentModel{Float64, Vector{Float64}})
   @ CaNNOLeS ~/.julia/packages/CaNNOLeS/EgvjT/src/CaNNOLeS.jl:103

Am I doing something wrong? Also, can CaNNOLeS optimize robustified least squares costs, or just pure least squares costs?

[abelsiqueira]

Hi @ojwoodford, thanks for reporting this. The issue here is that BundleAdjustmentModels doesn't have second-order information on the residual function, and CaNNOLeS tries to use that, by default.
One solution (would be) to use method = :Newton_noFHess in the CaNNOLeS call, but I just noticed that it doesn't work.
I created JuliaSmoothOptimizers/CaNNOLeS.jl#110 to track it and I hope to fix it soon.

I don't know the robustified least squares problem. What changes in the optimization problem?

[dpo]

Yes I guess it is possible to implement a thin modeling layer between the raw model and the solver if we provide the derivatives of the $\rho$ functions.

[ojwoodford]

Ceres is written based on the assumption that each group of residuals only depend on a subset of variables. This is for practical implementation purposes. However, this is not a necessary assumption to make. All residuals could easily depend on all variables. This has no effect on robustification.

One necessary assumption to make for robustification is that each rho function (or robustifier) applies to a subset of residuals (a subset can be 1 or more residuals), and that these subsets are mutually exclusive.

I am actually in the process of writing my own NLS optimizer in Julia here, which closely follows the Ceres way of doing things. It already has support for robust costs. It would be straightforward for me to add the solver described in your paper to this library. However, at present I'm not sure if the algorithm can be adapted for use with robustified costs, or not. It is not obvious that it can be. So my question is currently more theoretical than practical.

[AntoninKns]

Sorry for taking so long to react. I am currently not adding new features but I am upgrading documentation and solving the latest issues in order to help people to add new features and to make a proper release.
Concerning the idea of adding robustified NLS it could be a nice idea indeed, but I cannot do it right now.
That being said, if you have any questions about the code I would be glad to answer them.
Here is a link to the initial conception of the code and a detailed version of the jacobian matrix.

[ojwoodford]

The main change that needs to be made is in the solver. It should not assume that cost = 0.5 * r(x)' * r(x), and instead should call a separate cost function, that each model should supply. The model should supply 3 functions:

• f(x): The robustified cost (scalar, or vector to be summed)
• r(x): The weighted residual vector, length N
• J(x): The weighted Jacobian matrix, NxM, where M is the number of model variables

Then any model can easily robustified.

[tmigot]

Hi @ojwoodford ! Currently, most of our NLS API assumes that the cost function is 0.5 * r(x) * r(x), but using the function obj instead is possible. I am adding it to our modeling packages and will keep you updated when it is implemented in the solvers.

[ojwoodford]

Hi again. I can't find any example in the documention of a BundleAdjustmentModel actually being optimized. It would help to add this to the BundleAdjustmentModels tutorial. I'm again having trouble getting a BundleAdjustmentModel to optimize, now using any optimizer. See the following:

using BundleAdjustmentModel, CaNNOLeS, JSOSolvers
model = BundleAdjustmentModel("problem-16-22106")
stats = tron(model) 

ERROR: MethodError: no method matching jtprod_residual!(::BundleAdjustmentModel{Float64, Vector{…}}, ::Vector{Float64}, ::Vector{Float64}, ::Vector{Float64})
The function `jtprod_residual!` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  jtprod_residual!(::NLPModelsModifiers.FeasibilityFormNLS, ::AbstractVector, ::AbstractVector, ::AbstractVector)
   @ NLPModelsModifiers ~/.julia/packages/NLPModelsModifiers/zy8cI/src/feasibility-form-nls.jl:465
  jtprod_residual!(::NLPModelsModifiers.FeasibilityResidual, ::AbstractVector, ::AbstractVector, ::AbstractVector)
   @ NLPModelsModifiers ~/.julia/packages/NLPModelsModifiers/zy8cI/src/feasibility-residual.jl:119
  jtprod_residual!(::NLPModelsModifiers.SlackNLSModel, ::AbstractVector, ::AbstractVector, ::AbstractVector)
   @ NLPModelsModifiers ~/.julia/packages/NLPModelsModifiers/zy8cI/src/slack-model.jl:751
  ...

Stacktrace:
 [1] objgrad!(nls::BundleAdjustmentModel{…}, x::Vector{…}, g::Vector{…}, Fx::Vector{…}; recompute::Bool)
   @ NLPModels ~/.julia/packages/NLPModels/Cfwb3/src/nls/api.jl:556
 [2] solve!(solver::TrunkSolverNLS{…}, nlp::BundleAdjustmentModel{…}, stats::SolverCore.GenericExecutionStats{…}; callback::JSOSolvers.var"#79#80", x::Vector{…}, atol::Float64, rtol::Float64, Fatol::Float64, Frtol::Float64, max_eval::Int64, max_iter::Int64, max_time::Float64, verbose::Int64, subsolver_verbose::Int64)
   @ JSOSolvers ~/.julia/packages/JSOSolvers/tX6Sk/src/trunkls.jl:269
 [3] solve!(solver::TrunkSolverNLS{…}, model::BundleAdjustmentModel{…}; kwargs::@Kwargs{…})
   @ SolverCoreNLPModelsExt ~/.julia/packages/SolverCore/UI1dA/ext/SolverCoreNLPModelsExt.jl:60
 [4] trunk(::Val{…}, nlp::BundleAdjustmentModel{…}; x::Vector{…}, bk_max::Int64, monotone::Bool, nm_itmax::Int64, subsolver::Symbol, kwargs::@Kwargs{})
   @ JSOSolvers ~/.julia/packages/JSOSolvers/tX6Sk/src/trunkls.jl:220
 [5] trunk(::Val{:GaussNewton}, nlp::BundleAdjustmentModel{Float64, Vector{Float64}})
   @ JSOSolvers ~/.julia/packages/JSOSolvers/tX6Sk/src/trunkls.jl:203
 [6] trunk(nlp::BundleAdjustmentModel{Float64, Vector{Float64}}; variant::Symbol, kwargs::@Kwargs{})
   @ JSOSolvers ~/.julia/packages/JSOSolvers/tX6Sk/src/trunkls.jl:5
 [7] jso_func()
   @ Main ~/.julia/dev/scratch/compare_bundleadjusters.jl:7
 [8] top-level scope
   @ ~/.julia/dev/scratch/compare_bundleadjusters.jl:35
Some type information was truncated. Use `show(err)` to see complete types.

stats = trunk(model) 

ERROR: MethodError: no method matching jtprod_residual!(::BundleAdjustmentModel{Float64, Vector{…}}, ::Vector{Float64}, ::Vector{Float64}, ::Vector{Float64})
The function `jtprod_residual!` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  jtprod_residual!(::NLPModelsModifiers.FeasibilityFormNLS, ::AbstractVector, ::AbstractVector, ::AbstractVector)
   @ NLPModelsModifiers ~/.julia/packages/NLPModelsModifiers/zy8cI/src/feasibility-form-nls.jl:465
  jtprod_residual!(::NLPModelsModifiers.FeasibilityResidual, ::AbstractVector, ::AbstractVector, ::AbstractVector)
   @ NLPModelsModifiers ~/.julia/packages/NLPModelsModifiers/zy8cI/src/feasibility-residual.jl:119
  jtprod_residual!(::NLPModelsModifiers.SlackNLSModel, ::AbstractVector, ::AbstractVector, ::AbstractVector)
   @ NLPModelsModifiers ~/.julia/packages/NLPModelsModifiers/zy8cI/src/slack-model.jl:751
  ...

Stacktrace:
 [1] objgrad!(nls::BundleAdjustmentModel{…}, x::Vector{…}, g::Vector{…}, Fx::Vector{…}; recompute::Bool)
   @ NLPModels ~/.julia/packages/NLPModels/Cfwb3/src/nls/api.jl:556
 [2] solve!(solver::TrunkSolverNLS{…}, nlp::BundleAdjustmentModel{…}, stats::SolverCore.GenericExecutionStats{…}; callback::JSOSolvers.var"#79#80", x::Vector{…}, atol::Float64, rtol::Float64, Fatol::Float64, Frtol::Float64, max_eval::Int64, max_iter::Int64, max_time::Float64, verbose::Int64, subsolver_verbose::Int64)
   @ JSOSolvers ~/.julia/packages/JSOSolvers/tX6Sk/src/trunkls.jl:269
 [3] solve!(solver::TrunkSolverNLS{…}, model::BundleAdjustmentModel{…}; kwargs::@Kwargs{…})
   @ SolverCoreNLPModelsExt ~/.julia/packages/SolverCore/UI1dA/ext/SolverCoreNLPModelsExt.jl:60
 [4] trunk(::Val{…}, nlp::BundleAdjustmentModel{…}; x::Vector{…}, bk_max::Int64, monotone::Bool, nm_itmax::Int64, subsolver::Symbol, kwargs::@Kwargs{})
   @ JSOSolvers ~/.julia/packages/JSOSolvers/tX6Sk/src/trunkls.jl:220
 [5] trunk(::Val{:GaussNewton}, nlp::BundleAdjustmentModel{Float64, Vector{Float64}})
   @ JSOSolvers ~/.julia/packages/JSOSolvers/tX6Sk/src/trunkls.jl:203
 [6] trunk(nlp::BundleAdjustmentModel{Float64, Vector{Float64}}; variant::Symbol, kwargs::@Kwargs{})
   @ JSOSolvers ~/.julia/packages/JSOSolvers/tX6Sk/src/trunkls.jl:5
 [7] jso_func()
   @ Main ~/.julia/dev/scratch/compare_bundleadjusters.jl:7
 [8] top-level scope
   @ ~/.julia/dev/scratch/compare_bundleadjusters.jl:35
Some type information was truncated. Use `show(err)` to see complete types.

stats = cannoles(model, method=:Newton_noFHess; max_time=1000.0)

┌ Warning: linsolve ma57 not available. Using :ldlfactorizations instead
└ @ CaNNOLeS ~/.julia/packages/CaNNOLeS/XBbAl/src/CaNNOLeS.jl:318
ERROR: AssertionError: Dϕ < 0
Stacktrace:
 [1] line_search(x::Vector{…}, r::Vector{…}, λ::Vector{…}, dx::SubArray{…}, dr::SubArray{…}, dλ::Vector{…}, Fx::Vector{…}, cx::Vector{…}, Jx::SparseMatricesCOO.SparseMatrixCOO{…}, Jcx::SparseMatricesCOO.SparseMatrixCOO{…}, ϕ::CaNNOLeS.var"#ϕ#10", xt::Vector{…}, rt::Vector{…}, λt::Vector{…}, Ft::Vector{…}, ct::Vector{…}, F!::CaNNOLeS.var"#F!#8"{…}, c!::CaNNOLeS.var"#c!#9"{…}, ρ::Float64, δ::Float64, η::Float64, trial_computed::Bool, merit::Symbol, params::CaNNOLeS.ParamCaNNOLeS{…})
   @ CaNNOLeS ~/.julia/packages/CaNNOLeS/XBbAl/src/CaNNOLeS.jl:1085
 [2] solve!(solver::CaNNOLeSSolver{…}, nls::BundleAdjustmentModel{…}, stats::SolverCore.GenericExecutionStats{…}; callback::CaNNOLeS.var"#13#14", x::Vector{…}, λ::Vector{…}, use_initial_multiplier::Bool, max_iter::Int64, max_eval::Int64, max_time::Float64, max_inner::Int64, atol::Float64, rtol::Float64, Fatol::Float64, Frtol::Float64, verbose::Int64, always_accept_extrapolation::Bool, δdec::Float64)
   @ CaNNOLeS ~/.julia/packages/CaNNOLeS/XBbAl/src/CaNNOLeS.jl:672
 [3] solve!
   @ ~/.julia/packages/CaNNOLeS/XBbAl/src/CaNNOLeS.jl:418 [inlined]
 [4] #solve!#1
   @ ~/.julia/packages/SolverCore/UI1dA/ext/SolverCoreNLPModelsExt.jl:60 [inlined]
 [5] cannoles(nls::BundleAdjustmentModel{…}; linsolve::Symbol, method::Symbol, kwargs::@Kwargs{…})
   @ CaNNOLeS ~/.julia/packages/CaNNOLeS/XBbAl/src/CaNNOLeS.jl:415
 [6] cannoles
   @ ~/.julia/packages/CaNNOLeS/XBbAl/src/CaNNOLeS.jl:402 [inlined]
 [7] jso_func()
   @ Main ~/.julia/dev/scratch/compare_bundleadjusters.jl:7
 [8] top-level scope
   @ ~/.julia/dev/scratch/compare_bundleadjusters.jl:35
Some type information was truncated. Use `show(err)` to see complete types.

I also find cannoles takes way longer than expected. My own solver takes 3 seconds to optimize the problem, while cannoles seems to be taking over 30s per iteration!