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)

  398.073 ms (0 allocations: 0 bytes)
  23.085 ms (100000 allocations: 4.58 MiB)

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{NonlinearSolve.NonlinearSolveTag, Float64}}(0.13814225311887496,1.0,0.0,0.0,0.0),), 0x000000000001bcd7)
└ @ Main faq.md:62

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:

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}:
 1.0000000030823328
 0.10000000193952718
 1.9999999964201636
 0.49999999843942466

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

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.138142253118875
 0.18692555929884616
 1.8395605384738718
 0.4300583183043974

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. 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, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, SciMLBase.StandardNonlinearProblem}, ::GeneralizedFirstOrderAlgorithm{nothing, :NewtonRaphson, NoLineSearch, Missing, NewtonDescent{Nothing, typeof(NonlinearSolve.DEFAULT_PRECS)}, Nothing, Nothing, Nothing})
  from solve(prob::NonlinearProblem, args...; sensealg, u0, p, wrap, kwargs...) @ DiffEqBase ~/.julia/packages/DiffEqBase/eLhx9/src/solve.jl:1028
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, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, SciMLBase.StandardNonlinearProblem}
  args::Tuple{GeneralizedFirstOrderAlgorithm{nothing, :NewtonRaphson, NoLineSearch, Missing, NewtonDescent{Nothing, typeof(NonlinearSolve.DEFAULT_PRECS)}, Nothing, Nothing, Nothing}}
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, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, SciMLBase.StandardNonlinearProblem}, GeneralizedFirstOrderAlgorithm{nothing, :NewtonRaphson, NoLineSearch, Missing, NewtonDescent{Nothing, typeof(NonlinearSolve.DEFAULT_PRECS)}, Nothing, Nothing, Nothing}, Nothing, Nothing, NonlinearSolve.ImmutableNLStats, NonlinearSolve.NonlinearSolveTrace{false, false, TraceMinimal, Nothing}}
1 ─ %1 = DiffEqBase.:(var"#solve#41")::Core.Const(DiffEqBase.var"#solve#41")
│   %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, SymbolicIndexingInterface.SymbolCache{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, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, SciMLBase.StandardNonlinearProblem}, GeneralizedFirstOrderAlgorithm{nothing, :NewtonRaphson, NoLineSearch, Missing, NewtonDescent{Nothing, typeof(NonlinearSolve.DEFAULT_PRECS)}, Nothing, Nothing, Nothing}, Nothing, Nothing, NonlinearSolve.ImmutableNLStats, NonlinearSolve.NonlinearSolveTrace{false, false, TraceMinimal, Nothing}}
└──      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, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, SciMLBase.StandardNonlinearProblem}, ::GeneralizedFirstOrderAlgorithm{nothing, :NewtonRaphson, NoLineSearch, Missing, NewtonDescent{Nothing, typeof(NonlinearSolve.DEFAULT_PRECS)}, Nothing, Nothing, Nothing})
  from solve(prob::NonlinearProblem, args...; sensealg, u0, p, wrap, kwargs...) @ DiffEqBase ~/.julia/packages/DiffEqBase/eLhx9/src/solve.jl:1028
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, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, SciMLBase.StandardNonlinearProblem}
  args::Tuple{GeneralizedFirstOrderAlgorithm{nothing, :NewtonRaphson, NoLineSearch, Missing, NewtonDescent{Nothing, typeof(NonlinearSolve.DEFAULT_PRECS)}, Nothing, Nothing, Nothing}}
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, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, SciMLBase.StandardNonlinearProblem}, GeneralizedFirstOrderAlgorithm{nothing, :NewtonRaphson, NoLineSearch, Missing, NewtonDescent{Nothing, typeof(NonlinearSolve.DEFAULT_PRECS)}, Nothing, Nothing, Nothing}, Nothing, Nothing, NonlinearSolve.ImmutableNLStats, NonlinearSolve.NonlinearSolveTrace{false, false, TraceMinimal, Nothing}}
1 ─ %1 = DiffEqBase.:(var"#solve#41")::Core.Const(DiffEqBase.var"#solve#41")
│   %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, SymbolicIndexingInterface.SymbolCache{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, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, SciMLBase.StandardNonlinearProblem}, GeneralizedFirstOrderAlgorithm{nothing, :NewtonRaphson, NoLineSearch, Missing, NewtonDescent{Nothing, typeof(NonlinearSolve.DEFAULT_PRECS)}, Nothing, Nothing, Nothing}, Nothing, Nothing, NonlinearSolve.ImmutableNLStats, NonlinearSolve.NonlinearSolveTrace{false, false, TraceMinimal, Nothing}}
└──      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, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, SciMLBase.StandardNonlinearProblem}, ::GeneralizedFirstOrderAlgorithm{nothing, :NewtonRaphson, NoLineSearch, Missing, NewtonDescent{Nothing, typeof(NonlinearSolve.DEFAULT_PRECS)}, Nothing, Nothing, Nothing})
  from solve(prob::NonlinearProblem, args...; sensealg, u0, p, wrap, kwargs...) @ DiffEqBase ~/.julia/packages/DiffEqBase/eLhx9/src/solve.jl:1028
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, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, SciMLBase.StandardNonlinearProblem}
  args::Tuple{GeneralizedFirstOrderAlgorithm{nothing, :NewtonRaphson, NoLineSearch, Missing, NewtonDescent{Nothing, typeof(NonlinearSolve.DEFAULT_PRECS)}, Nothing, Nothing, Nothing}}
Body::SCIMLBASE.NONLINEARSOLUTION{FLOAT64, 1, VECTOR{FLOAT64}, VECTOR{FLOAT64}, NONLINEARPROBLEM{VECTOR{FLOAT64}, FALSE, FLOAT64, NONLINEARFUNCTION{FALSE, SCIMLBASE.FULLSPECIALIZE, TYPEOF(MAIN.__ATEXAMPLE__NAMED__TYPE_UNSTABLE.F), LINEARALGEBRA.UNIFORMSCALING{BOOL}, NOTHING, NOTHING, NOTHING, NOTHING, NOTHING, NOTHING, NOTHING, NOTHING, NOTHING, NOTHING, TYPEOF(SCIMLBASE.DEFAULT_OBSERVED_NO_TIME), NOTHING, SYMBOLICINDEXINGINTERFACE.SYMBOLCACHE{NOTHING, NOTHING, NOTHING}, NOTHING}, BASE.PAIRS{SYMBOL, UNION{}, TUPLE{}, @NAMEDTUPLE{}}, SCIMLBASE.STANDARDNONLINEARPROBLEM}, GENERALIZEDFIRSTORDERALGORITHM{NOTHING, :NEWTONRAPHSON, NOLINESEARCH, MISSING, NEWTONDESCENT{NOTHING, TYPEOF(NONLINEARSOLVE.DEFAULT_PRECS)}, NOTHING, NOTHING, NOTHING}, NOTHING, NOTHING, _A, NONLINEARSOLVE.NONLINEARSOLVETRACE{FALSE, FALSE, TRACEMINIMAL, NOTHING}} WHERE _A
1 ─ %1 = DiffEqBase.:(var"#solve#41")::Core.Const(DiffEqBase.var"#solve#41")
│   %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, SymbolicIndexingInterface.SymbolCache{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.__ATEXAMPLE__NAMED__TYPE_UNSTABLE.F), LINEARALGEBRA.UNIFORMSCALING{BOOL}, NOTHING, NOTHING, NOTHING, NOTHING, NOTHING, NOTHING, NOTHING, NOTHING, NOTHING, NOTHING, TYPEOF(SCIMLBASE.DEFAULT_OBSERVED_NO_TIME), NOTHING, SYMBOLICINDEXINGINTERFACE.SYMBOLCACHE{NOTHING, NOTHING, NOTHING}, NOTHING}, BASE.PAIRS{SYMBOL, UNION{}, TUPLE{}, @NAMEDTUPLE{}}, SCIMLBASE.STANDARDNONLINEARPROBLEM}, GENERALIZEDFIRSTORDERALGORITHM{NOTHING, :NEWTONRAPHSON, NOLINESEARCH, MISSING, NEWTONDESCENT{NOTHING, TYPEOF(NONLINEARSOLVE.DEFAULT_PRECS)}, NOTHING, NOTHING, NOTHING}, NOTHING, NOTHING, _A, NONLINEARSOLVE.NONLINEARSOLVETRACE{FALSE, FALSE, TRACEMINIMAL, NOTHING}} WHERE _A
└──      return %9

Oh no! This is type unstable. This is because ForwardDiff.jl will chunk the jacobian computation and the type of this chunksize can't be statically inferred. 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, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, SciMLBase.StandardNonlinearProblem}, ::GeneralizedFirstOrderAlgorithm{nothing, :NewtonRaphson, NoLineSearch, Missing, NewtonDescent{Nothing, typeof(NonlinearSolve.DEFAULT_PRECS)}, AutoForwardDiff{2, Nothing}, AutoForwardDiff{2, Nothing}, Nothing})
  from solve(prob::NonlinearProblem, args...; sensealg, u0, p, wrap, kwargs...) @ DiffEqBase ~/.julia/packages/DiffEqBase/eLhx9/src/solve.jl:1028
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, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, SciMLBase.StandardNonlinearProblem}
  args::Tuple{GeneralizedFirstOrderAlgorithm{nothing, :NewtonRaphson, NoLineSearch, Missing, NewtonDescent{Nothing, typeof(NonlinearSolve.DEFAULT_PRECS)}, AutoForwardDiff{2, Nothing}, AutoForwardDiff{2, Nothing}, Nothing}}
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, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, SciMLBase.StandardNonlinearProblem}, GeneralizedFirstOrderAlgorithm{nothing, :NewtonRaphson, NoLineSearch, Missing, NewtonDescent{Nothing, typeof(NonlinearSolve.DEFAULT_PRECS)}, AutoForwardDiff{2, Nothing}, AutoForwardDiff{2, Nothing}, Nothing}, Nothing, Nothing, NonlinearSolve.ImmutableNLStats, NonlinearSolve.NonlinearSolveTrace{false, false, TraceMinimal, Nothing}}
1 ─ %1 = DiffEqBase.:(var"#solve#41")::Core.Const(DiffEqBase.var"#solve#41")
│   %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, SymbolicIndexingInterface.SymbolCache{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, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, SciMLBase.StandardNonlinearProblem}, GeneralizedFirstOrderAlgorithm{nothing, :NewtonRaphson, NoLineSearch, Missing, NewtonDescent{Nothing, typeof(NonlinearSolve.DEFAULT_PRECS)}, AutoForwardDiff{2, Nothing}, AutoForwardDiff{2, Nothing}, Nothing}, Nothing, Nothing, NonlinearSolve.ImmutableNLStats, NonlinearSolve.NonlinearSolveTrace{false, false, TraceMinimal, Nothing}}
└──      return %9

And boom! Type stable again. We always recommend picking the chunksize via NonlinearSolve.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.

NonlinearSolve.pickchunksize — Function

pickchunksize(x) = pickchunksize(length(x))
pickchunksize(x::Int)

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