Frequently Asked Questions

How is the performance of Julia's NonlinearSolve.jl vs MATLAB's fzero?

This is addressed in a Twitter thread with the author of the improved fzero. On the test example:

using NonlinearSolve, BenchmarkTools

const N = 100_000;
levels = 1.5 .* rand(N);
out = zeros(N);
myfun(x, lv) = x * sin(x) - lv

function f(out, levels, u0)
    for i in 1:N
        out[i] = solve(
            IntervalNonlinearProblem{false}(
                IntervalNonlinearFunction{false}(myfun), u0, levels[i]),
            Falsi()).u
    end
end

function f2(out, levels, u0)
    for i in 1:N
        out[i] = solve(
            NonlinearProblem{false}(NonlinearFunction{false}(myfun), u0, levels[i]),
            SimpleNewtonRaphson()).u
    end
end

@btime f(out, levels, (0.0, 2.0))
@btime f2(out, levels, 1.0)
  97.684 ms (0 allocations: 0 bytes)
  32.080 ms (0 allocations: 0 bytes)

MATLAB 2022a achieves 1.66s. Try this code yourself: we receive 0.009 seconds, or a 184x speedup.

For more information on performance of SciML, see the SciMLBenchmarks.

The solver tried to set a Dual Number in my Vector of Floats. How do I fix that?

This is a common problem that occurs if the code was not written to be generic based on the input types. For example, consider this example taken from this issue

using NonlinearSolve, Random

function fff_incorrect(var, p)
    v_true = [1.0, 0.1, 2.0, 0.5]
    xx = [1.0, 2.0, 3.0, 4.0]
    xx[1] = var[1] - v_true[1]
    return var - v_true
end

v_true = [1.0, 0.1, 2.0, 0.5]
v_init = v_true .+ randn!(similar(v_true)) * 0.1

prob_oop = NonlinearLeastSquaresProblem{false}(fff_incorrect, v_init)
try
    sol = solve(prob_oop, LevenbergMarquardt(); maxiters = 10000, abstol = 1e-8)
catch e
    @error e
end
┌ Error: MethodError(Float64, (Dual{ForwardDiff.Tag{NonlinearFunction{false, SciMLBase.FullSpecialize, typeof(Main.fff_incorrect), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing}, Float64}}(-0.07049586366969296,1.0,0.0,0.0,0.0),), 0x0000000000007c61)
@ Main faq.md:65

Essentially what happened was, NonlinearSolve checked that we can use ForwardDiff.jl to differentiate the function based on the input types. However, this function has xx = [1.0, 2.0, 3.0, 4.0] followed by a xx[1] = var[1] - v_true[1] where var might be a Dual number. This causes the error. To fix it:

  1. Specify the autodiff to be AutoFiniteDiff

    sol = solve(prob_oop, LevenbergMarquardt(; autodiff = AutoFiniteDiff());
        maxiters = 10000, abstol = 1e-8)
    retcode: Success
    u: 4-element Vector{Float64}:
     0.9999999984270598
     0.09999999824750039
     2.000000000793593
     0.49999999856531624

    This worked but, Finite Differencing is not the recommended approach in any scenario.

  2. Rewrite the function to use PreallocationTools.jl or write it as

    function fff_correct(var, p)
        v_true = [1.0, 0.1, 2.0, 0.5]
        xx = eltype(var)[1.0, 2.0, 3.0, 4.0]
        xx[1] = var[1] - v_true[1]
        return xx - v_true
    end
    
    prob_oop = NonlinearLeastSquaresProblem{false}(fff_correct, v_init)
    sol = solve(prob_oop, LevenbergMarquardt(); maxiters = 10000, abstol = 1e-8)
    retcode: Stalled
    u: 4-element Vector{Float64}:
     1.9999999761145675
     0.021456628282225956
     2.035565906388018
     0.4357003833710344

I thought NonlinearSolve.jl was type-stable and fast. But it isn't, why?

It is hard to say why your code is not fast. Take a look at the Diagnostics API to pin-point the problem. One common issue is that there is type instability.

If you are using the defaults for the autodiff and your problem is not a scalar or using static arrays, ForwardDiff will create type unstable code and lead to dynamic dispatch internally. See this simple example:

using NonlinearSolve, InteractiveUtils

f(u, p) = @. u^2 - p

prob = NonlinearProblem{false}(f, 1.0, 2.0)

@code_warntype solve(prob, NewtonRaphson())
MethodInstance for CommonSolve.solve(::NonlinearProblem{Float64, false, Float64, NonlinearFunction{false, SciMLBase.FullSpecialize, typeof(Main.f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, SciMLBase.StandardNonlinearProblem}, ::GeneralizedFirstOrderAlgorithm{Missing, Missing, NewtonDescent{Nothing}, Nothing, Nothing, Nothing, Val{false}})
  from solve(prob::NonlinearProblem, args...; sensealg, u0, p, wrap, kwargs...) @ DiffEqBase ~/.julia/packages/DiffEqBase/JuOYZ/src/solve.jl:1104
Arguments
  #self#::Core.Const(CommonSolve.solve)
  prob::NonlinearProblem{Float64, false, Float64, NonlinearFunction{false, SciMLBase.FullSpecialize, typeof(Main.f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, SciMLBase.StandardNonlinearProblem}
  args::Tuple{GeneralizedFirstOrderAlgorithm{Missing, Missing, NewtonDescent{Nothing}, Nothing, Nothing, Nothing, Val{false}}}
Body::SciMLBase.NonlinearSolution{Float64, 0, Float64, Float64, NonlinearProblem{Float64, false, Float64, NonlinearFunction{false, SciMLBase.FullSpecialize, typeof(Main.f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, SciMLBase.StandardNonlinearProblem}, GeneralizedFirstOrderAlgorithm{Missing, Missing, NewtonDescent{Nothing}, AutoForwardDiff{nothing, Nothing}, AutoFiniteDiff{Val{:forward}, Val{:forward}, Val{:hcentral}, Nothing, Nothing, Bool}, AutoForwardDiff{nothing, Nothing}, Val{false}}, Nothing, Nothing, SciMLBase.NLStats, NonlinearSolveBase.NonlinearSolveTrace{Val{false}, Val{false}, Nothing, NonlinearSolveBase.NonlinearSolveTracing{Val{:minimal}}, NonlinearProblem{Float64, false, Float64, NonlinearFunction{false, SciMLBase.FullSpecialize, typeof(Main.f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, SciMLBase.StandardNonlinearProblem}}}
1 ─ %1 = DiffEqBase.:(var"#solve#43")::Core.Const(DiffEqBase.var"#solve#43")
 %2 = DiffEqBase.nothing::Core.Const(nothing)
 %3 = DiffEqBase.nothing::Core.Const(nothing)
 %4 = DiffEqBase.nothing::Core.Const(nothing)
 %5 = DiffEqBase.Val(true)::Core.Const(Val{true}())
 %6 = Core.NamedTuple()::Core.Const(NamedTuple())
 %7 = Base.pairs(%6)::Core.Const(Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}())
 %8 = Core.tuple(%2, %3, %4, %5, %7, #self#, prob)::Tuple{Nothing, Nothing, Nothing, Val{true}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, typeof(solve), NonlinearProblem{Float64, false, Float64, NonlinearFunction{false, SciMLBase.FullSpecialize, typeof(Main.f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, SciMLBase.StandardNonlinearProblem}}
 %9 = Core._apply_iterate(Base.iterate, %1, %8, args)::SciMLBase.NonlinearSolution{Float64, 0, Float64, Float64, NonlinearProblem{Float64, false, Float64, NonlinearFunction{false, SciMLBase.FullSpecialize, typeof(Main.f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, SciMLBase.StandardNonlinearProblem}, GeneralizedFirstOrderAlgorithm{Missing, Missing, NewtonDescent{Nothing}, AutoForwardDiff{nothing, Nothing}, AutoFiniteDiff{Val{:forward}, Val{:forward}, Val{:hcentral}, Nothing, Nothing, Bool}, AutoForwardDiff{nothing, Nothing}, Val{false}}, Nothing, Nothing, SciMLBase.NLStats, NonlinearSolveBase.NonlinearSolveTrace{Val{false}, Val{false}, Nothing, NonlinearSolveBase.NonlinearSolveTracing{Val{:minimal}}, NonlinearProblem{Float64, false, Float64, NonlinearFunction{false, SciMLBase.FullSpecialize, typeof(Main.f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, SciMLBase.StandardNonlinearProblem}}}
└──      return %9

Notice that this was type-stable, since it is a scalar problem. Now what happens for static arrays

using StaticArrays

prob = NonlinearProblem{false}(f, @SVector([1.0, 2.0]), 2.0)

@code_warntype solve(prob, NewtonRaphson())
MethodInstance for CommonSolve.solve(::NonlinearProblem{StaticArraysCore.SVector{2, Float64}, false, Float64, NonlinearFunction{false, SciMLBase.FullSpecialize, typeof(Main.f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, SciMLBase.StandardNonlinearProblem}, ::GeneralizedFirstOrderAlgorithm{Missing, Missing, NewtonDescent{Nothing}, Nothing, Nothing, Nothing, Val{false}})
  from solve(prob::NonlinearProblem, args...; sensealg, u0, p, wrap, kwargs...) @ DiffEqBase ~/.julia/packages/DiffEqBase/JuOYZ/src/solve.jl:1104
Arguments
  #self#::Core.Const(CommonSolve.solve)
  prob::NonlinearProblem{StaticArraysCore.SVector{2, Float64}, false, Float64, NonlinearFunction{false, SciMLBase.FullSpecialize, typeof(Main.f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, SciMLBase.StandardNonlinearProblem}
  args::Tuple{GeneralizedFirstOrderAlgorithm{Missing, Missing, NewtonDescent{Nothing}, Nothing, Nothing, Nothing, Val{false}}}
Body::SciMLBase.NonlinearSolution{Float64, 1, StaticArraysCore.SVector{2, Float64}, StaticArraysCore.SVector{2, Float64}, NonlinearProblem{StaticArraysCore.SVector{2, Float64}, false, Float64, NonlinearFunction{false, SciMLBase.FullSpecialize, typeof(Main.f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, SciMLBase.StandardNonlinearProblem}, GeneralizedFirstOrderAlgorithm{Missing, Missing, NewtonDescent{Nothing}, AutoForwardDiff{nothing, Nothing}, AutoFiniteDiff{Val{:forward}, Val{:forward}, Val{:hcentral}, Nothing, Nothing, Bool}, AutoForwardDiff{nothing, Nothing}, Val{false}}, Nothing, Nothing, SciMLBase.NLStats, NonlinearSolveBase.NonlinearSolveTrace{Val{false}, Val{false}, Nothing, NonlinearSolveBase.NonlinearSolveTracing{Val{:minimal}}, NonlinearProblem{StaticArraysCore.SVector{2, Float64}, false, Float64, NonlinearFunction{false, SciMLBase.FullSpecialize, typeof(Main.f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, SciMLBase.StandardNonlinearProblem}}}
1 ─ %1 = DiffEqBase.:(var"#solve#43")::Core.Const(DiffEqBase.var"#solve#43")
 %2 = DiffEqBase.nothing::Core.Const(nothing)
 %3 = DiffEqBase.nothing::Core.Const(nothing)
 %4 = DiffEqBase.nothing::Core.Const(nothing)
 %5 = DiffEqBase.Val(true)::Core.Const(Val{true}())
 %6 = Core.NamedTuple()::Core.Const(NamedTuple())
 %7 = Base.pairs(%6)::Core.Const(Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}())
 %8 = Core.tuple(%2, %3, %4, %5, %7, #self#, prob)::Tuple{Nothing, Nothing, Nothing, Val{true}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, typeof(solve), NonlinearProblem{StaticArraysCore.SVector{2, Float64}, false, Float64, NonlinearFunction{false, SciMLBase.FullSpecialize, typeof(Main.f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, SciMLBase.StandardNonlinearProblem}}
 %9 = Core._apply_iterate(Base.iterate, %1, %8, args)::SciMLBase.NonlinearSolution{Float64, 1, StaticArraysCore.SVector{2, Float64}, StaticArraysCore.SVector{2, Float64}, NonlinearProblem{StaticArraysCore.SVector{2, Float64}, false, Float64, NonlinearFunction{false, SciMLBase.FullSpecialize, typeof(Main.f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, SciMLBase.StandardNonlinearProblem}, GeneralizedFirstOrderAlgorithm{Missing, Missing, NewtonDescent{Nothing}, AutoForwardDiff{nothing, Nothing}, AutoFiniteDiff{Val{:forward}, Val{:forward}, Val{:hcentral}, Nothing, Nothing, Bool}, AutoForwardDiff{nothing, Nothing}, Val{false}}, Nothing, Nothing, SciMLBase.NLStats, NonlinearSolveBase.NonlinearSolveTrace{Val{false}, Val{false}, Nothing, NonlinearSolveBase.NonlinearSolveTracing{Val{:minimal}}, NonlinearProblem{StaticArraysCore.SVector{2, Float64}, false, Float64, NonlinearFunction{false, SciMLBase.FullSpecialize, typeof(Main.f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, SciMLBase.StandardNonlinearProblem}}}
└──      return %9

Again Type-Stable! Now let's try using a regular array:

prob = NonlinearProblem(f, [1.0, 2.0], 2.0)

@code_warntype solve(prob, NewtonRaphson())
MethodInstance for CommonSolve.solve(::NonlinearProblem{Vector{Float64}, false, Float64, NonlinearFunction{false, SciMLBase.FullSpecialize, typeof(Main.f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, SciMLBase.StandardNonlinearProblem}, ::GeneralizedFirstOrderAlgorithm{Missing, Missing, NewtonDescent{Nothing}, Nothing, Nothing, Nothing, Val{false}})
  from solve(prob::NonlinearProblem, args...; sensealg, u0, p, wrap, kwargs...) @ DiffEqBase ~/.julia/packages/DiffEqBase/JuOYZ/src/solve.jl:1104
Arguments
  #self#::Core.Const(CommonSolve.solve)
  prob::NonlinearProblem{Vector{Float64}, false, Float64, NonlinearFunction{false, SciMLBase.FullSpecialize, typeof(Main.f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, SciMLBase.StandardNonlinearProblem}
  args::Tuple{GeneralizedFirstOrderAlgorithm{Missing, Missing, NewtonDescent{Nothing}, Nothing, Nothing, Nothing, Val{false}}}
Body::SciMLBase.NonlinearSolution{Float64, 1, Vector{Float64}, Vector{Float64}, NonlinearProblem{Vector{Float64}, false, Float64, NonlinearFunction{false, SciMLBase.FullSpecialize, typeof(Main.f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, SciMLBase.StandardNonlinearProblem}, GeneralizedFirstOrderAlgorithm{Missing, Missing, NewtonDescent{Nothing}, AutoForwardDiff{nothing, Nothing}, AutoFiniteDiff{Val{:forward}, Val{:forward}, Val{:hcentral}, Nothing, Nothing, Bool}, AutoForwardDiff{nothing, Nothing}, Val{false}}, Nothing, Nothing, SciMLBase.NLStats, NonlinearSolveBase.NonlinearSolveTrace{Val{false}, Val{false}, Nothing, NonlinearSolveBase.NonlinearSolveTracing{Val{:minimal}}, NonlinearProblem{Vector{Float64}, false, Float64, NonlinearFunction{false, SciMLBase.FullSpecialize, typeof(Main.f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, SciMLBase.StandardNonlinearProblem}}}
1 ─ %1 = DiffEqBase.:(var"#solve#43")::Core.Const(DiffEqBase.var"#solve#43")
 %2 = DiffEqBase.nothing::Core.Const(nothing)
 %3 = DiffEqBase.nothing::Core.Const(nothing)
 %4 = DiffEqBase.nothing::Core.Const(nothing)
 %5 = DiffEqBase.Val(true)::Core.Const(Val{true}())
 %6 = Core.NamedTuple()::Core.Const(NamedTuple())
 %7 = Base.pairs(%6)::Core.Const(Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}())
 %8 = Core.tuple(%2, %3, %4, %5, %7, #self#, prob)::Tuple{Nothing, Nothing, Nothing, Val{true}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, typeof(solve), NonlinearProblem{Vector{Float64}, false, Float64, NonlinearFunction{false, SciMLBase.FullSpecialize, typeof(Main.f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, SciMLBase.StandardNonlinearProblem}}
 %9 = Core._apply_iterate(Base.iterate, %1, %8, args)::SciMLBase.NonlinearSolution{Float64, 1, Vector{Float64}, Vector{Float64}, NonlinearProblem{Vector{Float64}, false, Float64, NonlinearFunction{false, SciMLBase.FullSpecialize, typeof(Main.f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, SciMLBase.StandardNonlinearProblem}, GeneralizedFirstOrderAlgorithm{Missing, Missing, NewtonDescent{Nothing}, AutoForwardDiff{nothing, Nothing}, AutoFiniteDiff{Val{:forward}, Val{:forward}, Val{:hcentral}, Nothing, Nothing, Bool}, AutoForwardDiff{nothing, Nothing}, Val{false}}, Nothing, Nothing, SciMLBase.NLStats, NonlinearSolveBase.NonlinearSolveTrace{Val{false}, Val{false}, Nothing, NonlinearSolveBase.NonlinearSolveTracing{Val{:minimal}}, NonlinearProblem{Vector{Float64}, false, Float64, NonlinearFunction{false, SciMLBase.FullSpecialize, typeof(Main.f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, SciMLBase.StandardNonlinearProblem}}}
└──      return %9

Ah it is still type stable. But internally since the chunksize is not statically inferred, it will be dynamic and lead to dynamic dispatch. To fix this, we directly specify the chunksize:

@code_warntype solve(
    prob,
    NewtonRaphson(;
        autodiff = AutoForwardDiff(; chunksize = NonlinearSolve.pickchunksize(prob.u0))
    )
)
MethodInstance for CommonSolve.solve(::NonlinearProblem{Vector{Float64}, false, Float64, NonlinearFunction{false, SciMLBase.FullSpecialize, typeof(Main.f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, SciMLBase.StandardNonlinearProblem}, ::GeneralizedFirstOrderAlgorithm{Missing, Missing, NewtonDescent{Nothing}, AutoForwardDiff{2, Nothing}, Nothing, Nothing, Val{false}})
  from solve(prob::NonlinearProblem, args...; sensealg, u0, p, wrap, kwargs...) @ DiffEqBase ~/.julia/packages/DiffEqBase/JuOYZ/src/solve.jl:1104
Arguments
  #self#::Core.Const(CommonSolve.solve)
  prob::NonlinearProblem{Vector{Float64}, false, Float64, NonlinearFunction{false, SciMLBase.FullSpecialize, typeof(Main.f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, SciMLBase.StandardNonlinearProblem}
  args::Tuple{GeneralizedFirstOrderAlgorithm{Missing, Missing, NewtonDescent{Nothing}, AutoForwardDiff{2, Nothing}, Nothing, Nothing, Val{false}}}
Body::SciMLBase.NonlinearSolution{Float64, 1, Vector{Float64}, Vector{Float64}, NonlinearProblem{Vector{Float64}, false, Float64, NonlinearFunction{false, SciMLBase.FullSpecialize, typeof(Main.f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, SciMLBase.StandardNonlinearProblem}, GeneralizedFirstOrderAlgorithm{Missing, Missing, NewtonDescent{Nothing}, AutoForwardDiff{2, Nothing}, AutoFiniteDiff{Val{:forward}, Val{:forward}, Val{:hcentral}, Nothing, Nothing, Bool}, AutoForwardDiff{2, Nothing}, Val{false}}, Nothing, Nothing, SciMLBase.NLStats, NonlinearSolveBase.NonlinearSolveTrace{Val{false}, Val{false}, Nothing, NonlinearSolveBase.NonlinearSolveTracing{Val{:minimal}}, NonlinearProblem{Vector{Float64}, false, Float64, NonlinearFunction{false, SciMLBase.FullSpecialize, typeof(Main.f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, SciMLBase.StandardNonlinearProblem}}}
1 ─ %1 = DiffEqBase.:(var"#solve#43")::Core.Const(DiffEqBase.var"#solve#43")
 %2 = DiffEqBase.nothing::Core.Const(nothing)
 %3 = DiffEqBase.nothing::Core.Const(nothing)
 %4 = DiffEqBase.nothing::Core.Const(nothing)
 %5 = DiffEqBase.Val(true)::Core.Const(Val{true}())
 %6 = Core.NamedTuple()::Core.Const(NamedTuple())
 %7 = Base.pairs(%6)::Core.Const(Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}())
 %8 = Core.tuple(%2, %3, %4, %5, %7, #self#, prob)::Tuple{Nothing, Nothing, Nothing, Val{true}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, typeof(solve), NonlinearProblem{Vector{Float64}, false, Float64, NonlinearFunction{false, SciMLBase.FullSpecialize, typeof(Main.f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, SciMLBase.StandardNonlinearProblem}}
 %9 = Core._apply_iterate(Base.iterate, %1, %8, args)::SciMLBase.NonlinearSolution{Float64, 1, Vector{Float64}, Vector{Float64}, NonlinearProblem{Vector{Float64}, false, Float64, NonlinearFunction{false, SciMLBase.FullSpecialize, typeof(Main.f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, SciMLBase.StandardNonlinearProblem}, GeneralizedFirstOrderAlgorithm{Missing, Missing, NewtonDescent{Nothing}, AutoForwardDiff{2, Nothing}, AutoFiniteDiff{Val{:forward}, Val{:forward}, Val{:hcentral}, Nothing, Nothing, Bool}, AutoForwardDiff{2, Nothing}, Val{false}}, Nothing, Nothing, SciMLBase.NLStats, NonlinearSolveBase.NonlinearSolveTrace{Val{false}, Val{false}, Nothing, NonlinearSolveBase.NonlinearSolveTracing{Val{:minimal}}, NonlinearProblem{Vector{Float64}, false, Float64, NonlinearFunction{false, SciMLBase.FullSpecialize, typeof(Main.f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, Nothing, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, SciMLBase.StandardNonlinearProblem}}}
└──      return %9

And boom! Type stable again. We always recommend picking the chunksize via NonlinearSolveBase.pickchunksize, however, if you manually specify the chunksize, it must be ≤ length of input. However, a very large chunksize can lead to excessive compilation times and slowdown.

NonlinearSolveBase.pickchunksizeFunction
pickchunksize(x) = pickchunksize(length(x))
pickchunksize(x::Int)

Determine the chunk size for ForwardDiff and PolyesterForwardDiff based on the input length.

source