Differences between NLsolve and Optim in solving system of equations

Specific Domains Optimization (Mathematical)
21

Thank you!

22

Does that solver have some type instability?

using NLSolvers,ForwardDiff

function F_rosenbrock!(Fx, x)
    Fx[1] = 1 - x[1]
    Fx[2] = 10(x[2]-x[1]^2)
    return Fx
end

function customnlsolve(f!,x0,method=TrustRegion(Newton(), Dogleg()),options=NEqOptions())
    len = length(x0)
    xcache = zeros(eltype(x0),len)
    Fcache = zeros(eltype(x0),len)
    JCache = zeros(eltype(x0),len,len)
    jconfig = ForwardDiff.JacobianConfig(f!,x0,x0)
    function j!(J,x)
        #@show J
        ForwardDiff.jacobian!(J,f!,Fcache,x,jconfig)
    end
    function fj!(F,J,x) 
        #@show J,F
        ForwardDiff.jacobian!(J,f!,F,x,jconfig)
        F,J
    end
    
    function jv!(x)
        function JacV(Fv,v)
            ForwardDiff.jacobian!(JCache,f!,Fcache,v,jconfig)
            Fv .= Jcache * v
        end
        return LinearMap(JacV,length(x))
    end
    vectorobj = NLSolvers.VectorObjective(f!,j!,fj!,jv!)
    vectorprob = NEqProblem(vectorobj)
    res = solve(vectorprob, x0,method , options)
end

@code_warntype customnlsolve(F_rosenbrock!, [0., 0.])

julia> @code_warntype customnlsolve(F_rosenbrock!, [0., 0.])
Variables
  #self#::Core.Const(customnlsolve)
  f!::Core.Const(F_rosenbrock!)
  x0::Vector{Float64}

Body::NLSolvers.ConvergenceInfo{TrustRegion{Newton{Direct, typeof(NLSolvers.DefaultNewtonLinsolve), Nothing, Nothing}, Dogleg{Nothing}, NLSolvers.BTR{Nothing}}, _A, NEqOptions{Float64, Int64, Nothing}} where _A
1 ─ %1 = Main.Newton()::Core.Const(Newton{Direct, typeof(NLSolvers.DefaultNewtonLinsolve), Nothing, Nothing}(Direct(), NLSolvers.DefaultNewtonLinsolve, nothing, nothing))
│   %2 = Main.Dogleg()::Core.Const(Dogleg{Nothing}(nothing))
│   %3 = Main.TrustRegion(%1, %2)::Core.Const(TrustRegion{Newton{Direct, typeof(NLSolvers.DefaultNewtonLinsolve), Nothing, Nothing}, Dogleg{Nothing}, NLSolvers.BTR{Nothing}}(Newton{Direct, typeof(NLSolvers.DefaultNewtonLinsolve), Nothing, Nothing}(Direct(), NLSolvers.DefaultNewtonLinsolve, nothing, nothing), Dogleg{Nothing}(nothing), NLSolvers.BTR{Nothing}(nothing)))
│   %4 = Main.NEqOptions()::Core.PartialStruct(NEqOptions{Float64, Int64, Nothing}, Any[Core.Const(0.0), Core.Const(1.0e-8), Core.Const(1.0e-12), Int64, Core.Const(nothing)])
│   %5 = (#self#)(f!, x0, %3, %4)::NLSolvers.ConvergenceInfo{TrustRegion{Newton{Direct, typeof(NLSolvers.DefaultNewtonLinsolve), Nothing, Nothing}, Dogleg{Nothing}, NLSolvers.BTR{Nothing}}, _A, NEqOptions{Float64, Int64, Nothing}} where _A
└──      return %5

but if I delete the line that defines jconfig (and modify the subsequent internal functions), then becomes type stable:

julia> @code_warntype customnlsolve(F_rosenbrock!, [0., 0.])
Variables
  #self#::Core.Const(customnlsolve)
  f!::Core.Const(F_rosenbrock!)
  x0::Vector{Float64}

Body::NLSolvers.ConvergenceInfo{TrustRegion{Newton{Direct, typeof(NLSolvers.DefaultNewtonLinsolve), Nothing, Nothing}, Dogleg{Nothing}, NLSolvers.BTR{Nothing}}, NamedTuple{(:zero, :best_residual, :ρF0, :ρ2F0, :time, :iter), Tuple{Vector{Float64}, Vector{Float64}, Float64, Float64, Float64, Int64}}, NEqOptions{Float64, Int64, Nothing}}
1 ─ %1 = Main.Newton()::Core.Const(Newton{Direct, typeof(NLSolvers.DefaultNewtonLinsolve), Nothing, Nothing}(Direct(), NLSolvers.DefaultNewtonLinsolve, nothing, nothing))
│   %2 = Main.Dogleg()::Core.Const(Dogleg{Nothing}(nothing))
│   %3 = Main.TrustRegion(%1, %2)::Core.Const(TrustRegion{Newton{Direct, typeof(NLSolvers.DefaultNewtonLinsolve), Nothing, Nothing}, Dogleg{Nothing}, NLSolvers.BTR{Nothing}}(Newton{Direct, typeof(NLSolvers.DefaultNewtonLinsolve), Nothing, Nothing}(Direct(), NLSolvers.DefaultNewtonLinsolve, nothing, nothing), Dogleg{Nothing}(nothing), NLSolvers.BTR{Nothing}(nothing)))
│   %4 = Main.NEqOptions()::Core.PartialStruct(NEqOptions{Float64, Int64, Nothing}, Any[Core.Const(0.0), Core.Const(1.0e-8), Core.Const(1.0e-12), Int64, Core.Const(nothing)])
│   %5 = (#self#)(f!, x0, %3, %4)::NLSolvers.ConvergenceInfo{TrustRegion{Newton{Direct, typeof(NLSolvers.DefaultNewtonLinsolve), Nothing, Nothing}, Dogleg{Nothing}, NLSolvers.BTR{Nothing}}, NamedTuple{(:zero, :best_residual, :ρF0, :ρ2F0, :time, :iter), Tuple{Vector{Float64}, Vector{Float64}, Float64, Float64, Float64, Int64}}, NEqOptions{Float64, Int64, Nothing}}
└──      return %5
1 Like
23

Ah, yes!, of course, I remembered that on implace functions,ForwardDiff, is type stable,good to remember it again

24

but this is coming from jconfig, no?

25

Sorry for reviving this. I have to ask, isn’t that what GMM does? That is, GMM minimizes a quadratic form of the moment equations?