ksil / LFPSQP.jl

Julia implementation of Locally Feasibly Projected Sequential Quadratic Programming
MIT License
24 stars 2 forks source link

Error after first step #2

Closed rakshith95 closed 5 months ago

rakshith95 commented 1 year ago

Hello,

I'm trying to run an equality-constrained problem, but I get the following error:


   step |          f     ||c||      |Δf|    ||Δx||  |   S iter      res  |   M   iter  (pcg)  |        α  flag
--------------------------------------------------------------------------------------------------------------
      0 |  0.000e+00   2.4e+04                      |                    |                    |               
ERROR: MethodError: -(::ReverseDiff.TrackedArray{Float64, Float64, 1, Vector{Float64}, Vector{Float64}}, ::StaticArraysCore.SVector{2, Float64}) is ambiguous. Candidates:
  -(a::AbstractArray, b::StaticArraysCore.StaticArray) in StaticArrays at /home/rakshith/.julia/packages/StaticArrays/x7lS0/src/linalg.jl:17
  -(x::ReverseDiff.TrackedArray{V, D}, y::AbstractVector) where {V, D} in ReverseDiff at /home/rakshith/.julia/packages/ReverseDiff/YkVxM/src/derivatives/linalg/arithmetic.jl:107
  -(x::ReverseDiff.TrackedArray{V, D}, y::AbstractArray) where {V, D} in ReverseDiff at /home/rakshith/.julia/packages/ReverseDiff/YkVxM/src/derivatives/linalg/arithmetic.jl:107
Possible fix, define
  -(::ReverseDiff.TrackedArray{V, D}, ::StaticArraysCore.StaticArray{S, T, 1} where {T, S<:Tuple}) where {V, D}
Stacktrace:
  [1] f(x::ReverseDiff.TrackedArray{Float64, Float64, 1, Vector{Float64}, Vector{Float64}})
    @ Main ~/ws/VisualGeometryToolkit.jl/src/affine/optimize_aff_SV.jl:60
  [2] ReverseDiff.GradientTape(f::typeof(f), input::SubArray{Float64, 1, Vector{Float64}, Tuple{UnitRange{Int64}}, true}, cfg::ReverseDiff.GradientConfig{ReverseDiff.TrackedArray{Float64, Float64, 1, Vector{Float64}, Vector{Float64}}})
    @ ReverseDiff ~/.julia/packages/ReverseDiff/YkVxM/src/api/tape.jl:199
  [3] gradient!(result::SubArray{Float64, 1, Vector{Float64}, Tuple{UnitRange{Int64}}, true}, f::Function, input::SubArray{Float64, 1, Vector{Float64}, Tuple{UnitRange{Int64}}, true}, cfg::ReverseDiff.GradientConfig{ReverseDiff.TrackedArray{Float64, Float64, 1, Vector{Float64}, Vector{Float64}}})
    @ ReverseDiff ~/.julia/packages/ReverseDiff/YkVxM/src/api/gradients.jl:41
  [4] grad!
    @ ~/.julia/packages/LFPSQP/JBVQh/src/autodiff_generators.jl:8 [inlined]
  [5] optimize(f::typeof(f), grad!::LFPSQP.var"#grad!#7"{typeof(f), ReverseDiff.GradientConfig{ReverseDiff.TrackedArray{Float64, Float64, 1, Vector{Float64}, Vector{Float64}}}}, c!::typeof(c!), jac!::LFPSQP.var"#jac!#9"{typeof(c!), ReverseDiff.JacobianConfig{ReverseDiff.TrackedArray{Float64, Float64, 1, Vector{Float64}, Vector{Float64}}, Vector{ReverseDiff.TrackedReal{Float64, Float64, Nothing}}}}, hess_lag_vec!::LFPSQP.var"#hess_lag_vec!#12"{LFPSQP.var"#grad!#7"{typeof(f), ReverseDiff.GradientConfig{ReverseDiff.TrackedArray{ForwardDiff.Dual{nothing, Float64, 1}, ForwardDiff.Dual{nothing, Float64, 1}, 1, Vector{ForwardDiff.Dual{nothing, Float64, 1}}, Vector{ForwardDiff.Dual{nothing, Float64, 1}}}}}, LFPSQP.var"#jac!#9"{typeof(c!), ReverseDiff.JacobianConfig{ReverseDiff.TrackedArray{ForwardDiff.Dual{nothing, Float64, 1}, ForwardDiff.Dual{nothing, Float64, 1}, 1, Vector{ForwardDiff.Dual{nothing, Float64, 1}}, Vector{ForwardDiff.Dual{nothing, Float64, 1}}}, Vector{ReverseDiff.TrackedReal{ForwardDiff.Dual{nothing, Float64, 1}, ForwardDiff.Dual{nothing, Float64, 1}, Nothing}}}}, Int64, Vector{ForwardDiff.Dual{nothing, Float64, 1}}, Matrix{ForwardDiff.Dual{nothing, Float64, 1}}, Vector{ForwardDiff.Dual{nothing, Float64, 1}}, Vector{ForwardDiff.Dual{nothing, Float64, 1}}, Int64}, x0::Vector{Float64}, xl::Nothing, xu::Nothing, m::Int64, param::LFPSQPParams)
    @ LFPSQP ~/.julia/packages/LFPSQP/JBVQh/src/optimize.jl:259
  [6] optimize(f::Function, c!::Function, x0::Vector{Float64}, xl::Nothing, xu::Nothing, m::Int64, param::LFPSQPParams)
    @ LFPSQP ~/.julia/packages/LFPSQP/JBVQh/src/optimize.jl:103
  [7] optimize (repeats 2 times)
    @ ~/.julia/packages/LFPSQP/JBVQh/src/optimize.jl:108 [inlined]
  [8] aff_optimize(OE1::OrientedEllipse{Float64}, OE2::OrientedEllipse{Float64}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ Main ~/ws/VisualGeometryToolkit.jl/src/affine/optimize_aff_SV.jl:91
  [9] aff_optimize(OE1::OrientedEllipse{Float64}, OE2::OrientedEllipse{Float64})
    @ Main ~/ws/VisualGeometryToolkit.jl/src/affine/optimize_aff_SV.jl:69
 [10] top-level scope
    @ ~/ws/VisualGeometryToolkit.jl/src/affine/optimize_aff_SV.jl:104

Do you have any idea why this would happen?

ksil commented 1 year ago

Hi Rakshith, thanks for the comment! My strong guess is that your function is somehow incompatible with ReverseDiff, which is a common Julia automatic differentiation library that I also use in LFPSQP here. Without seeing the function, I would also guess that you may have some vector / storage you're using inside that function that isn't typed appropriately (check out this link: https://juliadiff.org/ReverseDiff.jl/limits/). If you can fix the function and test it successfully with ReverseDiff, I would recommend that. Otherwise, if you can obtain the gradient/Hessian of the function in another manner, you can feed that into LFPSQP. Hope this helps!

rakshith95 commented 1 year ago

Hello @ksil , thanks for the hint. I call the LinearAlgebra.norm() function in my objective, so that could be causing the problem. I will try it without the function call, and let you know if it works.

I was also wondering on what happens if I have a constrained optimization problem, and my initial guess is infeasible. Does the solver fail, or does it start from some other initial point? (zeros / closest feasible point / something else? )

Thank you for your work!

ksil commented 1 year ago

Theoretically, you should always start LFPSQP with a feasible point. That being said, because every step involves a solve, in practice, I have found that it works if the starting point is decently close to the constraint manifold. So, your mileage may vary. If you have issues, you can always perform one robust function solve in order to get an initial point for LFPSQP that is feasible.

rakshith95 commented 1 year ago

If you have issues, you can always perform one robust function solve in order to get an initial point for LFPSQP that is feasible.

This is sort of what I want to do with LFPSQP. I'm trying to get a feasible point closest ( by Frobenius norm) to my initial data. What do you mean by

one robust function solve ?

ksil commented 1 year ago

For example, if you're trying to minimize f(x) over constraints c(x) = 0, but you don't have an x that satisfies c(x) = 0, you can first solve for c(x) = 0 with some nonlinear solver or even with LFPSQP by minimizing || c(x) ||^2. Then, you could start the original optimization routine with the feasible point that was found.