sciml / nonlinearsolve.jl Goto Github PK

High-performance and differentiation-enabled nonlinear solvers (Newton methods), bracketed rootfinding (bisection, Falsi), with sparsity and Newton-Krylov support.

Home Page: https://docs.sciml.ai/NonlinearSolve/stable/

License: MIT License

Julia 100.00%

nonlinear-equations steady-state equilibrium deep-equilibrium-models bracketing newton-raphson sparse-matrices high-performance-computing differential-equations factorization

nonlinearsolve.jl's Introduction

NonlinearSolve.jl

Fast implementations of root finding algorithms in Julia that satisfy the SciML common interface.

For information on using the package, see the stable documentation. Use the in-development documentation for the version of the documentation which contains the unreleased features.

High Level Examples

using NonlinearSolve, StaticArrays

f(u, p) = u .* u .- 2
u0 = @SVector[1.0, 1.0]
prob = NonlinearProblem(f, u0)
solver = solve(prob)

## Bracketing Methods

f(u, p) = u .* u .- 2.0
u0 = (1.0, 2.0) # brackets
prob = IntervalNonlinearProblem(f, u0)
sol = solve(prob)

Citation

If you found this library to be useful in academic work, then please cite:

@article{pal2024nonlinearsolve,
  title={NonlinearSolve. jl: High-Performance and Robust Solvers for Systems of Nonlinear Equations in Julia},
  author={Pal, Avik and Holtorf, Flemming and Larsson, Axel and Loman, Torkel and Schaefer, Frank and Qu, Qingyu and Edelman, Alan and Rackauckas, Chris and others},
  journal={arXiv preprint arXiv:2403.16341},
  year={2024}
}

nonlinearsolve.jl's People

Contributors

Stargazers

Watchers

Forkers

achuchmala ranocha briceonwiley utkarsh530 mkg33 briochemc kronosthelate ericphanson frankschae arnostrouwen vpuri3 fieryswampshire oscardssmith deltadahl daviehh chriselrod yash2798 avik-pal dawbarton erikqqy fholtorf axla-io mfkiwl danielwe simonbyrne sathvikbhagavan yonatanwesen oxinabox aayushsabharwal pepijndevos rzyu45 sebastianm-c termi-official alexandrechern tansongchen ven-k spinachboul gdalle m-bossart

nonlinearsolve.jl's Issues

TagBot trigger issue

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

Incompatibility with CuArrays and ForwardDiff

MWE:

using CUDA
using NonlinearSolve
f(u,p) = u .* u .- 2
u0 = CUDA.CuArray([1.0, 1.0])
probN = NonlinearProblem{false}(f, u0)
solver = solve(probN, NewtonRaphson(), tol = 1e-9)

FixedPointAcceleration.jl

https://github.com/s-baumann/FixedPointAcceleration.jl

Deprecation warning with mutating functions

using NonlinearSolve

function f!(F, x, _)
    F[1] = (x[1] + 3) * (x[2]^3 - 7) + 18
    F[2] = sin(x[2] * exp(x[1]) - 1)
    nothing
end

x0 = [0.1; 1.2]

prob! = NonlinearProblem{true}(f!, x0)

sol! = solve(prob!, NewtonRaphson())

that last line yield this warning:

Warning: `vendor()` is deprecated, use `BLAS.get_config()` and inspect the output instead
│   caller = (::NonlinearSolve.DefaultLinSolve)(x::Vector{Float64}, A::Matrix{Float64}, b::Vector{Float64}, update_matrix::Bool; tol::Nothing, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}) at utils.jl:125
└ @ NonlinearSolve ~/.julia/packages/NonlinearSolve/9GzK0/src/utils.jl:125
u: 2-element Vector{Float64}:
 -3.645318273517821e-8
  1.0000000239478837

Make iip a compile-time constant

Trying to add SciMLNLSolve results in an error

I am following the documentation here http://nonlinearsolve.sciml.ai/dev/solvers/NonlinearSystemSolvers/#SciMLNLSolve.jl but it results in an error:

(@v1.7) pkg> add SciMLNLSolve
    Updating registry at `~/.julia/registries/General.toml`
ERROR: The following package names could not be resolved:
 * SciMLNLSolve (not found in project, manifest or registry)

wrap NLSolvers.jl?

Maybe NonlinearSolve.jl should wrap https://github.com/JuliaNLSolvers/NLSolvers.jl?

Failed to precompile

Since yesterday:

ERROR: LoadError: Failed to precompile NonlinearSolve [8913a72c-1f9b-4ce2-8d82-65094dcecaec] to /....

(@v1.8) pkg> st
Status `/usr/share/julia/environments/v1.8/Project.toml`
  [a2441757] Coverage v1.6.0
  [8913a72c] NonlinearSolve v1.2.0

chunk_size argument to NewtonRaphson not respected

It just ends up calling plain ForwardDiff.jacobian without specifying the chunk size:

https://github.com/SciML/NonlinearSolve.jl/blob/master/src/jacobian.jl#L31-L38

jacobian erroneously calls finite_difference_derivative

This call should go to FiniteDiff.finite_difference_jacobian:

https://github.com/JuliaComputing/NonlinearSolve.jl/blob/cd00b8e4e3302be2cc724de126f5ca184d0e252e/src/jacobian.jl#L35

The reason the tests still pass is that the solution to the test problem is proportional to the initial guess, and the test function works for scalar u. If you change the initial guess so that it's not proportional to [1.0, 1.0] the algorithm will fail to converge, and if you use a function that fails for scalar u (e.g., by indexing into it) an exception will be thrown when trying to compute the jacobian.

Docs not found

The "docs|stable" pill at the top of the README.md points to http://nonlinearsolve.sciml.ai/stable/, but this link does not exist, and gives a 404 error

Specializing Broyden's method for sparse Jacobians

https://www.ams.org/journals/mcom/1970-24-109/S0025-5718-1970-0258276-9/

This might be putting it in the realm of NonlinearSolve instead, to keep the SimpleNonlinearSolve.jl version simple?

tutorial Iterator interface broken

using NonlinearSolve
f(u, p) = u .* u .- 2.0
u0 = (1.0, 2.0) # brackets
probB = NonlinearProblem(f, u0)
solver = init(probB, Falsi()) # Can iterate the solver object

TrustRegion correctness

We should really implement all the tests in https://people.sc.fsu.edu/~jburkardt/m_src/test_nonlin/test_nonlin.html to validate our solvers.

Here's the result of the second test problem:

NewtonRaphson residual norm: 8.514332152287816e-16
TrustRegion residual norm: 0.0008969098117197063
NLsolve residual norm: 8.496216993896772e-16

The full code is:

using NonlinearSolve, NLsolve, LinearAlgebra

function powell_singular_function!(out, x, p = nothing)
    out[1] = x[1] + 10.0 * x[2]
    out[2] = sqrt(5.0) * ( x[3] - x[4] )
    out[3] = ( x[2] - 2.0 * x[3] )^2
    out[4] = sqrt(10.0) * ( x[1] - x[4] ) * ( x[1] - x[4] )
    nothing
end

n = 4
x = zeros(n)
x[1]   = -1.2
x[2:n] .= 1.0
nlprob = NonlinearProblem(powell_singular_function!, x)
out = similar(x)
powell_singular_function!(out, solve(nlprob, NewtonRaphson(), abstol=1e-15, reltol=1e-15).u, nothing)
println("NewtonRaphson residual norm: $(norm(out))")
powell_singular_function!(out, solve(nlprob, TrustRegion(), abstol=1e-15, reltol=1e-15).u, nothing)
println("TrustRegion residual norm: $(norm(out))")
sol = nlsolve(powell_singular_function!, x, xtol=1e-15, ftol=1e-15)
println("NLsolve residual norm: $(sol.residual_norm)")

Long using time

Not sure if this is actually an issue, but it seemed like a significant time to me so thought I would ask about it.

I tried out the new @time_imports macro in v1.8.0-beta1 and found that almost half of the load time in ControlSystems came from NonlinarSolve. Also ran it on NonlinearSolve to see how it looked there, and it was faster than the time reported in ControlSystems (3.99s vs 6.58s).

I'm not sure if there are things that can be done easily in NonlinearSolve? And maybe it is really ControlSystems that needs to do something since the load time was worse there?

ControlSystems

julia> @time_imports using ControlSystems
     10.1 ms          ┌ Preferences
     24.9 ms        ┌ JLLWrappers
     30.2 ms      ┌ Rmath_jll
    147.0 ms    ┌ Rmath
      3.3 ms    ┌ Compat
      0.4 ms    ┌ NaNMath
      2.7 ms    ┌ Calculus
     83.4 ms        ┌ ChainRulesCore
     84.5 ms      ┌ ChangesOfVariables
      0.5 ms      ┌ OpenLibm_jll
      0.5 ms      ┌ InverseFunctions
      4.2 ms      ┌ DocStringExtensions
      6.8 ms      ┌ IrrationalConstants
      0.8 ms      ┌ CompilerSupportLibraries_jll
      1.6 ms        ┌ LogExpFunctions
      1.2 ms        ┌ OpenSpecFun_jll
     24.8 ms      ┌ SpecialFunctions
     21.4 ms      ┌ DualNumbers
    146.2 ms    ┌ HypergeometricFunctions
      0.4 ms    ┌ Reexport
    312.0 ms  ┌ StatsFuns
    713.7 ms    ┌ StaticArrays
    717.4 ms  ┌ DiffResults
     16.8 ms  ┌ FunctionWrappers
     12.8 ms    ┌ OrderedCollections
    109.0 ms  ┌ DataStructures
      0.4 ms  ┌ MuladdMacro
      0.5 ms    ┌ Requires
      0.1 ms      ┌ IfElse
     61.0 ms        ┌ Static
    290.3 ms      ┌ ArrayInterface
     18.3 ms      ┌ MacroTools
    345.6 ms    ┌ LabelledArrays
      0.6 ms      ┌ CommonSubexpressions
      1.1 ms      ┌ DiffRules
    114.6 ms    ┌ ForwardDiff
      0.8 ms    ┌ Adapt
    546.1 ms  ┌ PreallocationTools
     25.6 ms  ┌ KLU
    145.2 ms    ┌ Parsers
     12.7 ms    ┌ InlineStrings
      0.8 ms    ┌ ExprTools
     43.0 ms    ┌ RecipesBase
      2.6 ms    ┌ Mocking
     90.4 ms      ┌ TimeZones
    127.2 ms    ┌ Intervals
    436.6 ms    ┌ MutableArithmetics
   1003.6 ms  ┌ Polynomials
      2.3 ms    ┌ DataAPI
      0.9 ms    ┌ StatsAPI
      0.9 ms    ┌ SortingAlgorithms
     14.0 ms    ┌ Missings
     58.7 ms  ┌ StatsBase
     17.5 ms  ┌ AbstractFFTs
      0.6 ms    ┌ ZygoteRules
    845.5 ms        ┌ Hwloc_jll
    848.2 ms      ┌ Hwloc
      0.2 ms        ┌ SIMDTypes
      2.5 ms        ┌ ManualMemory
      8.0 ms      ┌ LayoutPointers
    354.0 ms      ┌ CPUSummary
      0.4 ms        ┌ BitTwiddlingConvenienceFunctions
     38.3 ms      ┌ HostCPUFeatures
   1953.7 ms    ┌ VectorizationBase
      0.2 ms      ┌ IteratorInterfaceExtensions
      0.9 ms    ┌ TableTraits
     73.8 ms      ┌ ThreadingUtilities
    101.3 ms    ┌ PolyesterWeave
      6.7 ms      ┌ SLEEFPirates
     13.0 ms    ┌ SIMDDualNumbers
      0.3 ms    ┌ UnPack
    237.9 ms    ┌ FillArrays
      0.3 ms    ┌ DataValueInterfaces
     67.9 ms    ┌ RecursiveArrayTools
      5.7 ms      ┌ CloseOpenIntervals
     25.1 ms      ┌ StrideArraysCore
     33.8 ms    ┌ Polyester
     27.0 ms      ┌ Tables
      0.3 ms      ┌ CommonSolve
      2.4 ms      ┌ ConstructionBase
      0.4 ms      ┌ TreeViews
    410.2 ms    ┌ SciMLBase
   2508.7 ms    ┌ Setfield
     73.7 ms      ┌ OffsetArrays
    658.2 ms    ┌ LoopVectorization
     21.0 ms    ┌ FiniteDiff
    318.5 ms    ┌ TriangularSolve
     32.1 ms    ┌ IterativeSolvers
    208.6 ms    ┌ RecursiveFactorization
   6583.1 ms  ┌ NonlinearSolve
      1.5 ms  ┌ LaTeXStrings
     10.8 ms    ┌ Distances
     12.8 ms    ┌ DEDataArrays
     42.8 ms    ┌ PDMats
      1.2 ms    ┌ Parameters
     17.9 ms    ┌ Krylov
     11.8 ms      ┌ NLSolversBase
     26.5 ms    ┌ LineSearches
      1.1 ms    ┌ Inflate
      5.7 ms      ┌ DensityInterface
     10.2 ms      ┌ QuadGK
    323.0 ms    ┌ Distributions
     61.0 ms    ┌ KrylovKit
    158.8 ms    ┌ LinearSolve
      0.4 ms    ┌ FastClosures
      1.0 ms    ┌ FastBroadcast
      5.3 ms    ┌ NLsolve
      2.6 ms        ┌ SimpleTraits
      9.3 ms          ┌ ArnoldiMethod
     83.1 ms        ┌ Graphs
     87.4 ms      ┌ VertexSafeGraphs
    114.6 ms    ┌ SparseDiffTools
    170.4 ms    ┌ DiffEqBase
     29.0 ms    ┌ ExponentialUtilities
   2545.5 ms  ┌ OrdinaryDiffEq
    290.4 ms  ┌ LinearMaps
    241.6 ms    ┌ FFTW_jll
    551.2 ms  ┌ FFTW
     58.9 ms    ┌ IterTools
    146.1 ms  ┌ DSP
    169.7 ms    ┌ FixedPointNumbers
    725.9 ms  ┌ ColorTypes
     21.3 ms  ┌ MatrixEquations
     28.5 ms  ┌ DelayDiffEq
     14.3 ms  ┌ DiffEqCallbacks
      7.0 ms  ┌ MatrixPencils
    578.5 ms  ┌ Colors
  14473.5 ms  ControlSystems

NonlinearSolve

julia> @time_imports using NonlinearSolve
     20.5 ms    ┌ MacroTools
     33.5 ms  ┌ ZygoteRules
      0.3 ms    ┌ IteratorInterfaceExtensions
      0.8 ms  ┌ TableTraits
    727.9 ms    ┌ StaticArrays
    731.4 ms  ┌ DiffResults
      1.8 ms  ┌ Compat
      0.5 ms  ┌ NaNMath
      0.4 ms  ┌ Requires
      0.2 ms  ┌ UnPack
    203.7 ms  ┌ FillArrays
      0.3 ms  ┌ DataValueInterfaces
     99.0 ms    ┌ ChainRulesCore
    100.0 ms  ┌ ChangesOfVariables
      4.1 ms    ┌ DocStringExtensions
      0.2 ms    ┌ IfElse
     30.1 ms    ┌ RecipesBase
     72.9 ms      ┌ Static
    258.0 ms    ┌ ArrayInterface
      0.6 ms    ┌ Adapt
    385.2 ms  ┌ RecursiveArrayTools
      0.5 ms  ┌ OpenLibm_jll
      0.4 ms  ┌ InverseFunctions
      0.2 ms  ┌ Reexport
      1.6 ms      ┌ DataAPI
     19.3 ms    ┌ Tables
      0.3 ms    ┌ CommonSolve
      1.3 ms    ┌ ConstructionBase
      0.3 ms    ┌ TreeViews
    358.6 ms  ┌ SciMLBase
      7.4 ms    ┌ IrrationalConstants
      0.7 ms    ┌ CompilerSupportLibraries_jll
      1.3 ms      ┌ LogExpFunctions
      9.9 ms          ┌ Preferences
     10.7 ms        ┌ JLLWrappers
    898.9 ms      ┌ OpenSpecFun_jll
    930.9 ms    ┌ SpecialFunctions
      1.1 ms    ┌ CommonSubexpressions
      2.0 ms    ┌ DiffRules
   1080.2 ms  ┌ ForwardDiff
   1038.4 ms  ┌ Setfield
     19.2 ms  ┌ FiniteDiff
     22.6 ms  ┌ IterativeSolvers
      0.6 ms  ┌ RecursiveFactorization
   3991.1 ms  NonlinearSolve

Add TrustRegion to Precompilation

Finishing #116 (comment) found a weird bug in Julia v1.6's precompilation that seems to be fixed. The best thing to do might be to just enable TrustRegion precompilation for version higher like v1.7 and call it a day. Since v1.6 doesn't have all of the extra precompilation features anyways it won't matter all that much that it's missing, but v1.9 certainly needs it.

Test of low overhead mode

using NonlinearSolve
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(NonlinearProblem{false}(NonlinearFunction{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]), NewtonRaphson()).u
    end
end
@time f(out, levels, (0.0, 2.0))
@time f2(out, levels, 1.0)

Support MVectors

see here:
https://discourse.julialang.org/t/help-with-nonlinearsolve-jl/74735

Odd numerical instability in solving a nonlinear damped pendulum BVP collocation problem

I notice an odd numerical instability in solving a simple nonlinear collocation problem. So the problem is taken from Toby Driscoll's book on numerical methods, chap 10, example 10.4.2. The equation is:

$θ'' + 0.05θ' + sinθ = 0$ with $θ(0) = 2.5, θ(8.0) = -2.0$.

The code is below. The problem is that the NonLinearSolve results are really inconsistent. I discretize the problem, but
sometimes the solver will give me the result, and other times the solver gives a Singular Matrix error. For example I might
discretize with 50 points and it gives me a Singular Matrix error, but with 49 points it gives me a solution. The numbers 49 and 50 are just make up for explanatory purposes. And the other odd thing is that sometimes at a given n, it will give a result and at other times it will either give the Singular Matrix error or just return a vector if Inf. So not sure what the issue is.

Here is the code.

using NonlinearSolve, BandedMatrices, ComponentArrays, Plots

function first_derivative_matrix(n, h)
    D = BandedMatrix{Float64}(undef, (n,n), (2,2))
    D[band(0)] .= 0
    D[band(1)] .= 0.5
    D[band(-1)] .= -0.5
    D[1, 1:3] = [-1.5, 2.0, -0.5]
    D[n, end-2:end] = [0.5, -2.0, 1.5]
    D = D/h   
end

function second_derivative_matrix(n, h)
    D = BandedMatrix{Float64}(undef, (n,n), (4,4))
    D[band(0)] .= -2.0
    D[band(1)] .= 1.0
    D[band(-1)] .= 1.0
    D[1, 1:4] = [2, -5, 4, -1]
    D[n, end-3:end] = [-1, 4, -5, 2]
    D = D/h^2   
end


function make_discrete(tspan, n)
    x = LinRange(tspan[1], tspan[2], n)
    h = (tspan[2] - tspan[1])/n
    return (x, h)
end

function objective(u, p)
    Dxx, Dx, h = p.Dxx, p.Dx, p.h
    f = Dxx*u .+ 0.05*Dx*u .+ sin.(u)
    f[1] = (u[1] - 2.5)/h^2
    f[end] = (u[end] + 2.0)/h^2
    return f
end

n = 100
tspan = (0.0, 8.0)
Θ,h = make_discrete(tspan, n)
Dx = first_derivative_matrix(n, h)
Dxx = second_derivative_matrix(n, h) 

p = ComponentArray(Dx = Dx, Dxx = Dxx, h = h) 

u0 = ones(n)
prob = NonlinearProblem{false}(objective, u0, p)
solver = solve(prob, NewtonRaphson(), tol=1e-8)

plot(Θ, solver)

Try n from a list of [10, 20, 30, 40], and at least one of them should give a solution and one of them should give Singular Matrix errors. Could this be a conditioning problem? In the objective function I divided f[1], f[end] the boundary condition by h^2 to improve the conditioning.

package loading time regression

Looks like it's due to #121 and #127

With version 1.1.1

@time_imports using NonlinearSolve

gives

142.4 ms NonlinearSolve

while with the 1.2.0 version, this shot up to ~13 seconds,

13191.1 ms NonlinearSolve

Tested w/

julia> versioninfo()
Julia Version 1.8.5
Commit 17cfb8e65ea (2023-01-08 06:45 UTC)
Platform Info:
  OS: macOS (arm64-apple-darwin21.5.0)
  CPU: 10 × Apple M1 Pro
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-13.0.1 (ORCJIT, apple-m1)
  Threads: 1 on 8 virtual cores

By reverting the changes in #127, the long using NonlinearSolve time is fixed. So maybe it would be better to revert #127 until it can be figured out where it went wrong?

Thanks!

What does `false` in `NonlinearProblem{false}` mean?

The README is a bit sparse... and the docstring doesn't help:

help?> NonlinearProblem
search: NonlinearProblem

  struct NonlinearProblem{uType, isinplace, P, F, K} <: SciMLBase.AbstractNonlinearProblem{uType, isinplace}

  ───────────────────────────────────────────────────────

  NonlinearProblem(f, u0)
  NonlinearProblem(f, u0, p; kwargs...)
  

  Define a steady state problem using an instance of
  AbstractNonlinearFunction.

  ───────────────────────────────────────────────────────

  NonlinearProblem(prob)
  

  Define a steady state problem from a standard ODE
  problem.

GPU Compatibility for 1D Arrays

It seems that (parts of) the library won't work with 1D CUDA arrays, but with higher dimensional arrays it does.

using NonlinearSolve, CUDA
CUDA.allowscalar(false) 
f(u,p) = u .* u .- 2
u0 = [1.0, 1.0]
probN = NonlinearProblem{false}(f, u0)
solver = solve(probN, NewtonRaphson(), tol = 1e-9)

works

u0 = reshape(cu(u0),:,1)
probN = NonlinearProblem{false}(f, u0)
solver = @time solve(probN, NewtonRaphson(), tol = 1e-9)

works

u0 = cu([1.0, 1.0])
probN = NonlinearProblem{false}(f, u0)
solver = solve(probN, NewtonRaphson(), tol = 1e-9)

doesn't work, with error:

MethodError: no method matching copy_cublasfloat(::CuArray{Float32, 1, CUDA.Mem.DeviceBuffer})
Closest candidates are:
  copy_cublasfloat(::CuArray{T, 2}) where T at ~/.julia/packages/CUDA/Uurn4/lib/cusolver/linalg.jl:6

Stacktrace:
 [1] \(_A::CuArray{Float32, 2, CUDA.Mem.DeviceBuffer}, _B::CuArray{Float32, 1, CUDA.Mem.DeviceBuffer})
   @ CUDA.CUSOLVER ~/.julia/packages/CUDA/Uurn4/lib/cusolver/linalg.jl:24
 [2] macro expansion
   @ ~/.julia/packages/Setfield/AS2xF/src/sugar.jl:197 [inlined]
 [3] perform_step(solver::NonlinearSolve.NewtonImmutableSolver{NonlinearFunction{false, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing}, NewtonRaphson{12, true, DataType, NonlinearSolve.DefaultLinSolve}, CuArray{Float32, 1, CUDA.Mem.DeviceBuffer}, CuArray{Float32, 1, CUDA.Mem.DeviceBuffer}, SciMLBase.NullParameters, typeof(NonlinearSolve.DEFAULT_NORM), Float64, Nothing, NonlinearProblem{CuArray{Float32, 1, CUDA.Mem.DeviceBuffer}, false, SciMLBase.NullParameters, NonlinearFunction{false, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}}, alg::NewtonRaphson{12, true, DataType, NonlinearSolve.DefaultLinSolve}, #unused#::Val{false})
   @ NonlinearSolve ~/.julia/packages/NonlinearSolve/ZxbEJ/src/raphson.jl:58
 [4] solve!(solver::NonlinearSolve.NewtonImmutableSolver{NonlinearFunction{false, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing}, NewtonRaphson{12, true, DataType, NonlinearSolve.DefaultLinSolve}, CuArray{Float32, 1, CUDA.Mem.DeviceBuffer}, CuArray{Float32, 1, CUDA.Mem.DeviceBuffer}, SciMLBase.NullParameters, typeof(NonlinearSolve.DEFAULT_NORM), Float64, Nothing, NonlinearProblem{CuArray{Float32, 1, CUDA.Mem.DeviceBuffer}, false, SciMLBase.NullParameters, NonlinearFunction{false, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}})
   @ NonlinearSolve ~/.julia/packages/NonlinearSolve/ZxbEJ/src/solve.jl:63
 [5] #solve#6
   @ ~/.julia/packages/NonlinearSolve/ZxbEJ/src/solve.jl:5 [inlined]
 [6] top-level scope
   @ In[7]:3
 [7] eval
   @ ./boot.jl:373 [inlined]
 [8] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
   @ Base ./loading.jl:1196

Error in TrustRegion when inplace is true

Found an error in the new TrustRegion solver when iip = true.
This went trough testing because the tests where @test all(sol.u .* sol.u .- 2 .< 1e-9) and not @test all(abs.(sol.u .* sol.u .- 2) .< 1e-9).
Trying to fix it as I'm writing this, and hopefully have a fix within the hour.

SpeedMapping.jl wrapper

Incompatibility with units

Here's what I tried to do:

using ModelingToolkit,Unitful,NonlinearSolve
vars = @variables τ_lawson 
pars = @parameters τ_ratio τ_slow
eqs = [
    τ_ratio ~ τ_slow/τ_lawson,
]

ns = NonlinearSystem(eqs, vars, pars)
ps =[τ_ratio => 1.0,
     τ_slow=>1.0u"s"]
guess = [τ_lawson=>0.1u"s"]
prob = NonlinearProblem(ns,guess,ps)
sol = solve(prob,NewtonRaphson())

and then this happened:

julia> sol = solve(prob,NewtonRaphson())
ERROR: DimensionError: s and 10.0 are not dimensionally compatible.
Stacktrace:
  [1] convert(#unused#::Type{Quantity{Float64, 𝐓, Unitful.FreeUnits{(s,), 𝐓, nothing}}}, x::Float64)
    @ Unitful ~/.julia/packages/Unitful/yI5QN/src/conversion.jl:112
  [2] setindex!(A::Vector{Quantity{Float64, 𝐓, Unitful.FreeUnits{(s,), 𝐓, nothing}}}, x::Float64, i1::Int64)
    @ Base /Applications/Julia-1.6.app/Contents/Resources/julia/share/julia/base/array.jl:839
  [3] macro expansion
    @ ~/.julia/packages/SymbolicUtils/SXYR9/src/code.jl:325 [inlined]
  [4] macro expansion
    @ ~/.julia/packages/Symbolics/H8dtg/src/build_function.jl:359 [inlined]
  [5] macro expansion
    @ ~/.julia/packages/SymbolicUtils/SXYR9/src/code.jl:283 [inlined]
  [6] macro expansion
    @ ~/.julia/packages/RuntimeGeneratedFunctions/3SZ1T/src/RuntimeGeneratedFunctions.jl:124 [inlined]
  [7] macro expansion
    @ ./none:0 [inlined]
  [8] generated_callfunc
    @ ./none:0 [inlined]
  [9] (::RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x73c43a35, 0x06c47ba7, 0xb1d644a0, 0xa2947d04, 0xc57c3159)})(::Vector{Quantity{Float64, 𝐓, Unitful.FreeUnits{(s,), 𝐓, nothing}}}, ::Vector{Quantity{Float64, 𝐓, Unitful.FreeUnits{(s,), 𝐓, nothing}}}, ::Vector{Quantity{Float64, D, U} where {D, U}})
    @ RuntimeGeneratedFunctions ~/.julia/packages/RuntimeGeneratedFunctions/3SZ1T/src/RuntimeGeneratedFunctions.jl:112
 [10] f
    @ ~/.julia/packages/ModelingToolkit/VigaA/src/systems/nonlinear/nonlinearsystem.jl:156 [inlined]
 [11] NonlinearFunction
    @ ~/.julia/packages/SciMLBase/1aTqd/src/scimlfunctions.jl:342 [inlined]
 [12] init(::NonlinearProblem{Vector{Quantity{Float64, 𝐓, Unitful.FreeUnits{(s,), 𝐓, nothing}}}, true, Vector{Quantity{Float64, D, U} where {D, U}}, NonlinearFunction{true, ModelingToolkit.var"#f#302"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xd260deed, 0x7dc90363, 0x08fe8abd, 0xd742811b, 0x0b2c6f74)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x73c43a35, 0x06c47ba7, 0xb1d644a0, 0xa2947d04, 0xc57c3159)}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, ModelingToolkit.var"#generated_observed#305"{NonlinearSystem, Dict{Any, Any}}, Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, ::NewtonRaphson{12, true, DataType, NonlinearSolve.DefaultLinSolve}; alias_u0::Bool, maxiters::Int64, tol::Float64, internalnorm::Function, kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ NonlinearSolve ~/.julia/packages/NonlinearSolve/WjsHK/src/solve.jl:52
 [13] init
    @ ~/.julia/packages/NonlinearSolve/WjsHK/src/solve.jl:43 [inlined]
 [14] solve(::NonlinearProblem{Vector{Quantity{Float64, 𝐓, Unitful.FreeUnits{(s,), 𝐓, nothing}}}, true, Vector{Quantity{Float64, D, U} where {D, U}}, NonlinearFunction{true, ModelingToolkit.var"#f#302"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xd260deed, 0x7dc90363, 0x08fe8abd, 0xd742811b, 0x0b2c6f74)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x73c43a35, 0x06c47ba7, 0xb1d644a0, 0xa2947d04, 0xc57c3159)}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, ModelingToolkit.var"#generated_observed#305"{NonlinearSystem, Dict{Any, Any}}, Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, ::NewtonRaphson{12, true, DataType, NonlinearSolve.DefaultLinSolve}; kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ NonlinearSolve ~/.julia/packages/NonlinearSolve/WjsHK/src/solve.jl:4
 [15] solve(::NonlinearProblem{Vector{Quantity{Float64, 𝐓, Unitful.FreeUnits{(s,), 𝐓, nothing}}}, true, Vector{Quantity{Float64, D, U} where {D, U}}, NonlinearFunction{true, ModelingToolkit.var"#f#302"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xd260deed, 0x7dc90363, 0x08fe8abd, 0xd742811b, 0x0b2c6f74)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x73c43a35, 0x06c47ba7, 0xb1d644a0, 0xa2947d04, 0xc57c3159)}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, ModelingToolkit.var"#generated_observed#305"{NonlinearSystem, Dict{Any, Any}}, Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, ::NewtonRaphson{12, true, DataType, NonlinearSolve.DefaultLinSolve})
    @ NonlinearSolve ~/.julia/packages/NonlinearSolve/WjsHK/src/solve.jl:4
 [16] top-level scope
    @ REPL[87]:1

It looks like the culprit is in #init#10(alias_u0, maxiters, tol, internalnorm, kwargs, , prob, alg, args) at /Users/lmorton/.julia/packages/NonlinearSolve/WjsHK/src/solve.jl:35

   50    if iip
> 51      fu = zero(u)
   52      f(fu, u, p)

Here fu is:

Quantity{Float64, 𝐓, Unitful.FreeUnits{(s,), 𝐓, nothing}}[0.0 s]

which propagates all the way down to the setindex! call that's trying to store the dimensionless result there.

Adjoint overload

It's just the implicit function theorem

https://github.com/mitmath/18335/blob/spring20/notes/adjoint/adjoint.pdf

While there is a generic one for SteadyStateProblem in DiffEqSensitivity (that could be extended to NonlinearProblem). SciML/SciMLSensitivity.jl#389 will work on anything with __solve

Implement Levenburg-Marquedt

https://docs.scipy.org/doc/scipy/reference/optimize.root-lm.html

https://en.wikipedia.org/wiki/Levenberg%E2%80%93Marquardt_algorithm

Use DiffEqBase DEFAULT_LINSOLVE

This isn't quite possible because it would be a circular dependency, so we'd to do something a little trickier. But we should plan to keep the two swappable, and in the future find a way to share the same default.

Consistent segfault when precompiling on Julia 1.8, Mac M1

When precompiling DifferentialEquations, I get the following error, which seems to come from NonlinearSolve.jl. I'm not sure if it's a NonlinearSolve issue or a Julia issue, but I thought I'd post it here for reference.

Error message here

[ Info: Precompiling DifferentialEquations [0c46a032-eb83-5123-abaf-570d42b7fbaa]

signal (11): Segmentation fault: 11
in expression starting at /Users/anshul/.julia/packages/NonlinearSolve/7hPLq/src/NonlinearSolve.jl:15
jl_decode_value at /Applications/Julia-1.8.app/Contents/Resources/julia/lib/julia/libjulia-internal.1.8.dylib (unknown line)
jl_decode_value at /Applications/Julia-1.8.app/Contents/Resources/julia/lib/julia/libjulia-internal.1.8.dylib (unknown line)
jl_decode_value_array at /Applications/Julia-1.8.app/Contents/Resources/julia/lib/julia/libjulia-internal.1.8.dylib (unknown line)
ijl_uncompress_ir at /Applications/Julia-1.8.app/Contents/Resources/julia/lib/julia/libjulia-internal.1.8.dylib (unknown line)
_Z16jl_emit_codeinstP19_jl_code_instance_tP15_jl_code_info_tR20_jl_codegen_params_t at /Applications/Julia-1.8.app/Contents/Resources/julia/lib/julia/libjulia-codegen.1.8.dylib (unknown line)
_Z20jl_compile_workqueueRNSt3__13mapIP19_jl_code_instance_tNS_5tupleIJNS_10unique_ptrIN4llvm6ModuleENS_14default_deleteIS6_EEEE20_jl_llvm_functions_tEEENS_4lessIS2_EENS_9allocatorINS_4pairIKS2_SB_EEEEEER20_jl_codegen_params_t17CompilationPolicy at /Applications/Julia-1.8.app/Contents/Resources/julia/lib/julia/libjulia-codegen.1.8.dylib (unknown line)
_ZL20_jl_compile_codeinstP19_jl_code_instance_tP15_jl_code_info_tm at /Applications/Julia-1.8.app/Contents/Resources/julia/lib/julia/libjulia-codegen.1.8.dylib (unknown line)
jl_generate_fptr_impl at /Applications/Julia-1.8.app/Contents/Resources/julia/lib/julia/libjulia-codegen.1.8.dylib (unknown line)
jl_compile_method_internal at /Applications/Julia-1.8.app/Contents/Resources/julia/lib/julia/libjulia-internal.1.8.dylib (unknown line)
ijl_apply_generic at /Applications/Julia-1.8.app/Contents/Resources/julia/lib/julia/libjulia-internal.1.8.dylib (unknown line)
jl_f__call_in_world_total at /Applications/Julia-1.8.app/Contents/Resources/julia/lib/julia/libjulia-internal.1.8.dylib (unknown line)
ijl_apply_generic at /Applications/Julia-1.8.app/Contents/Resources/julia/lib/julia/libjulia-internal.1.8.dylib (unknown line)
do_apply at /Applications/Julia-1.8.app/Contents/Resources/julia/lib/julia/libjulia-internal.1.8.dylib (unknown line)
concrete_eval_call at ./compiler/abstractinterpretation.jl:737
abstract_call_method_with_const_args at ./compiler/abstractinterpretation.jl:780
abstract_call_gf_by_type at ./compiler/abstractinterpretation.jl:165
abstract_call_known at ./compiler/abstractinterpretation.jl:1666
abstract_call at ./compiler/abstractinterpretation.jl:1724
abstract_call at ./compiler/abstractinterpretation.jl:1703
abstract_eval_statement at ./compiler/abstractinterpretation.jl:1845
typeinf_local at ./compiler/abstractinterpretation.jl:2310
typeinf_nocycle at ./compiler/abstractinterpretation.jl:2406
_typeinf at ./compiler/typeinfer.jl:230
typeinf at ./compiler/typeinfer.jl:213
typeinf_edge at ./compiler/typeinfer.jl:876
abstract_call_method at ./compiler/abstractinterpretation.jl:632
abstract_call_gf_by_type at ./compiler/abstractinterpretation.jl:156
abstract_call_known at ./compiler/abstractinterpretation.jl:1666
abstract_call at ./compiler/abstractinterpretation.jl:1724
abstract_call at ./compiler/abstractinterpretation.jl:1703
abstract_eval_statement at ./compiler/abstractinterpretation.jl:1845
typeinf_local at ./compiler/abstractinterpretation.jl:2310
typeinf_nocycle at ./compiler/abstractinterpretation.jl:2406
_typeinf at ./compiler/typeinfer.jl:230
typeinf at ./compiler/typeinfer.jl:213
typeinf_edge at ./compiler/typeinfer.jl:876
abstract_call_method at ./compiler/abstractinterpretation.jl:632
abstract_call_gf_by_type at ./compiler/abstractinterpretation.jl:156
abstract_call_known at ./compiler/abstractinterpretation.jl:1666
abstract_call at ./compiler/abstractinterpretation.jl:1724
abstract_call at ./compiler/abstractinterpretation.jl:1703
abstract_eval_statement at ./compiler/abstractinterpretation.jl:1845
typeinf_local at ./compiler/abstractinterpretation.jl:2284
typeinf_nocycle at ./compiler/abstractinterpretation.jl:2406
_typeinf at ./compiler/typeinfer.jl:230
typeinf at ./compiler/typeinfer.jl:213
typeinf_edge at ./compiler/typeinfer.jl:876
abstract_call_method at ./compiler/abstractinterpretation.jl:632
abstract_call_gf_by_type at ./compiler/abstractinterpretation.jl:156
abstract_call_known at ./compiler/abstractinterpretation.jl:1666
abstract_call at ./compiler/abstractinterpretation.jl:1724
abstract_call at ./compiler/abstractinterpretation.jl:1703
abstract_eval_statement at ./compiler/abstractinterpretation.jl:1845
typeinf_local at ./compiler/abstractinterpretation.jl:2284
typeinf_nocycle at ./compiler/abstractinterpretation.jl:2406
_typeinf at ./compiler/typeinfer.jl:230
typeinf at ./compiler/typeinfer.jl:213
typeinf_ext at ./compiler/typeinfer.jl:957
typeinf_ext_toplevel at ./compiler/typeinfer.jl:990
typeinf_ext_toplevel at ./compiler/typeinfer.jl:986
jfptr_typeinf_ext_toplevel_15079 at /Applications/Julia-1.8.app/Contents/Resources/julia/lib/julia/sys.dylib (unknown line)
ijl_apply_generic at /Applications/Julia-1.8.app/Contents/Resources/julia/lib/julia/libjulia-internal.1.8.dylib (unknown line)
jl_type_infer at /Applications/Julia-1.8.app/Contents/Resources/julia/lib/julia/libjulia-internal.1.8.dylib (unknown line)
jl_generate_fptr_impl at /Applications/Julia-1.8.app/Contents/Resources/julia/lib/julia/libjulia-codegen.1.8.dylib (unknown line)
jl_compile_method_internal at /Applications/Julia-1.8.app/Contents/Resources/julia/lib/julia/libjulia-internal.1.8.dylib (unknown line)
ijl_apply_generic at /Applications/Julia-1.8.app/Contents/Resources/julia/lib/julia/libjulia-internal.1.8.dylib (unknown line)
@reexport at /Users/anshul/.julia/packages/Reexport/OxbHO/src/Reexport.jl:4
jl_invoke_julia_macro at /Applications/Julia-1.8.app/Contents/Resources/julia/lib/julia/libjulia-internal.1.8.dylib (unknown line)
jl_expand_macros at /Applications/Julia-1.8.app/Contents/Resources/julia/lib/julia/libjulia-internal.1.8.dylib (unknown line)
ijl_expand_stmt_with_loc at /Applications/Julia-1.8.app/Contents/Resources/julia/lib/julia/libjulia-internal.1.8.dylib (unknown line)
jl_toplevel_eval_flex at /Applications/Julia-1.8.app/Contents/Resources/julia/lib/julia/libjulia-internal.1.8.dylib (unknown line)
jl_toplevel_eval_flex at /Applications/Julia-1.8.app/Contents/Resources/julia/lib/julia/libjulia-internal.1.8.dylib (unknown line)
ijl_toplevel_eval_in at /Applications/Julia-1.8.app/Contents/Resources/julia/lib/julia/libjulia-internal.1.8.dylib (unknown line)
eval at ./boot.jl:368 [inlined]
include_string at ./loading.jl:1277
ijl_apply_generic at /Applications/Julia-1.8.app/Contents/Resources/julia/lib/julia/libjulia-internal.1.8.dylib (unknown line)
_include at ./loading.jl:1334
include at ./Base.jl:422 [inlined]
include_package_for_output at ./loading.jl:1400
jfptr_include_package_for_output_46348 at /Applications/Julia-1.8.app/Contents/Resources/julia/lib/julia/sys.dylib (unknown line)
ijl_apply_generic at /Applications/Julia-1.8.app/Contents/Resources/julia/lib/julia/libjulia-internal.1.8.dylib (unknown line)
do_call at /Applications/Julia-1.8.app/Contents/Resources/julia/lib/julia/libjulia-internal.1.8.dylib (unknown line)
eval_body at /Applications/Julia-1.8.app/Contents/Resources/julia/lib/julia/libjulia-internal.1.8.dylib (unknown line)
jl_interpret_toplevel_thunk at /Applications/Julia-1.8.app/Contents/Resources/julia/lib/julia/libjulia-internal.1.8.dylib (unknown line)
jl_toplevel_eval_flex at /Applications/Julia-1.8.app/Contents/Resources/julia/lib/julia/libjulia-internal.1.8.dylib (unknown line)
jl_toplevel_eval_flex at /Applications/Julia-1.8.app/Contents/Resources/julia/lib/julia/libjulia-internal.1.8.dylib (unknown line)
ijl_toplevel_eval_in at /Applications/Julia-1.8.app/Contents/Resources/julia/lib/julia/libjulia-internal.1.8.dylib (unknown line)
jlplt_ijl_toplevel_eval_in_13971 at /Applications/Julia-1.8.app/Contents/Resources/julia/lib/julia/sys.dylib (unknown line)
eval at ./boot.jl:368 [inlined]
include_string at ./loading.jl:1277
include_string at ./loading.jl:1287
ijl_apply_generic at /Applications/Julia-1.8.app/Contents/Resources/julia/lib/julia/libjulia-internal.1.8.dylib (unknown line)
exec_options at ./client.jl:301
_start at ./client.jl:522
jfptr__start_58807 at /Applications/Julia-1.8.app/Contents/Resources/julia/lib/julia/sys.dylib (unknown line)
ijl_apply_generic at /Applications/Julia-1.8.app/Contents/Resources/julia/lib/julia/libjulia-internal.1.8.dylib (unknown line)
true_main at /Applications/Julia-1.8.app/Contents/Resources/julia/lib/julia/libjulia-internal.1.8.dylib (unknown line)
jl_repl_entrypoint at /Applications/Julia-1.8.app/Contents/Resources/julia/lib/julia/libjulia-internal.1.8.dylib (unknown line)
Allocations: 5929456 (Pool: 5927121; Big: 2335); GC: 3
ERROR: LoadError: Failed to precompile NonlinearSolve [8913a72c-1f9b-4ce2-8d82-65094dcecaec] to /Users/anshul/.julia/compiled/v1.8/NonlinearSolve/jl_tbYC6N.
Stacktrace:
  [1] error(s::String)
    @ Base ./error.jl:35
  [2] compilecache(pkg::Base.PkgId, path::String, internal_stderr::IO, internal_stdout::IO, ignore_loaded_modules::Bool)
    @ Base ./loading.jl:1551
  [3] compilecache
    @ ./loading.jl:1495 [inlined]
  [4] _require(pkg::Base.PkgId)
    @ Base ./loading.jl:1199
  [5] _require_prelocked(uuidkey::Base.PkgId)
    @ Base ./loading.jl:1084
  [6] macro expansion
    @ ./loading.jl:1064 [inlined]
  [7] macro expansion
    @ ./lock.jl:223 [inlined]
  [8] require(into::Module, mod::Symbol)
    @ Base ./loading.jl:1028
  [9] include
    @ ./Base.jl:422 [inlined]
 [10] include_package_for_output(pkg::Base.PkgId, input::String, depot_path::Vector{String}, dl_load_path::Vector{String}, load_path::Vector{String}, concrete_deps::Vector{Pair{Base.PkgId, UInt64}}, source::String)
    @ Base ./loading.jl:1400
 [11] top-level scope
    @ stdin:1
in expression starting at /Users/anshul/.julia/packages/DiffEqBase/U3LtB/src/DiffEqBase.jl:1

julia> versioninfo()
Julia Version 1.8.0-beta3
Commit 3e092a2521 (2022-03-29 15:42 UTC)
Platform Info:
  OS: macOS (arm64-apple-darwin21.3.0)
  CPU: 10 × Apple M1 Max
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-13.0.1 (ORCJIT, apple-m1)
  Threads: 8 on 8 virtual cores
Environment:
  JULIA_NUM_THREADS = 8

(@v1.8) pkg> st -m NonlinearSolve
Status `~/.julia/environments/v1.8/Manifest.toml`
  [8913a72c] NonlinearSolve v0.3.16

Jacobian and Hessian Free Halley's Method

https://www.autodiff.org/Docs/euroad/Trond_Steihaug.pdf describes a bit how it would be done via jvps and hvps, which would mix in AD tricks.

`DEFAULT_PRECS` is not exported by NonlinearSolve

In #86 you refer to DEFAULT_PRECS but that does not seem to be exported https://github.com/SciML/NonlinearSolve.jl/search?q=DEFAULT_PRECS . Found this in CI of another package.

Maxiters seems to be offset by 1

Check out the following example:

julia> using NonlinearSolve, LinearAlgebra

julia> f(u,p) = u .* u .- p;

julia> u0 = [1.0, 1.0];

julia> p = 2.0;

julia> probN = NonlinearProblem{false}(f, u0, p);

julia> bla0 = solve(probN, NewtonRaphson(), maxiters = 0).u;

julia> bla1 = solve(probN, NewtonRaphson(), maxiters = 1).u;

julia> bla2 = solve(probN, NewtonRaphson(), maxiters = 2).u;

julia> norm(f(u0, p))
1.4142135623730951

julia> norm(f(bla0, p))
1.4142135623730951

julia> norm(f(bla1, p))
1.4142135623730951

julia> norm(f(bla2, p))
0.3535533905932738

It seems as if setting maxter to 1 acutally results in zero iterations occuring. And I am left to wonder if maxter=2 then means 2 iterations, or actually 1 iteration...

This is the same for my actual application, which makes me more certain that the issue is not specific to this problem. The other example is:

julia> using NonlinearSolve, LinearAlgebra

julia> v = 120;  #kg

julia> k = 2.5;  #m

julia> w = 4.0;  #km/m

julia> α = 2e-7; #kg⁻¹

julia> function F(x⃗, parameters)
           L₀, L, p, x, θ, φ, a, H = x⃗
           (;d, n) = parameters
           return [
           a*(cosh(x/a)-1) - p,
           2a * sinh(x/a) - L ,
           2x + 2k*cos(θ) - d,
           p+k*sin(θ) - n,
           sinh(x/a) - tan(φ),
           (1+v/(w*L₀)) * tan(φ) - tan(θ),
           L₀ * (1 + α*H) - L,
           w*L₀ / 2sin(φ) - H,
           ]
       end;

julia> params = (d=30, n=5);

julia> u0 =    [29, 27, 1, params.d/2, deg2rad(45),  deg2rad(22.5), 40, 500];

julia> probN = NonlinearProblem{false}(F, u0, params);

julia> bla0 = solve(probN, NewtonRaphson(), maxiters = 0).u;

julia> bla1 = solve(probN, NewtonRaphson(), maxiters = 1).u;

julia> bla2 = solve(probN, NewtonRaphson(), maxiters = 2).u;

julia> norm(F(u0, params))
348.4941911206226

julia> norm(F(bla0, params))
348.4941911206226

julia> norm(F(bla1, params))
348.4941911206226

julia> norm(F(bla2, params))
3.865006154302766

advanced docs custom linear solver

similar to LinearSolve

Error when using a neural network inside a system of nonlinear equations.

I've set up the following very simple system of nonlinear equations to solve.

using NonlinearSolve, Lux, LinearAlgebra, StableRNGs

rng = StableRNG(1111)

function func(u, p)
    [u[1].*u[1] .- p[1],
    u[2].*u[2] .- p[2]]
end

guess = [1., 1.]
params = [2.,3.]

prob = NonlinearProblem{false}(func, guess, params)
solution = solve(prob, NewtonRaphson())

This works fine. I am now trying to learn one of the terms from this system of equations using a neural network (similar to the UDE approach https://docs.sciml.ai/Overview/dev/showcase/missing_physics/).

I define a network.

ann = Lux.Chain(
    Lux.Dense(2,20,tanh), Lux.Dense(20,20,tanh), Lux.Dense(20,1)
)
p, st = Lux.setup(rng, ann)

And define the grey-box model, where I replace the u[2]*u[2] term with the output of the network.

function GreyBox(u, p)
    z = ann(u, p, st)[1]
    [u[1].*u[1] .- params[1],
    z[1] .- params[2]]
end

NN_prob = NonlinearProblem{false}(GreyBox, guess, p)
NN_solution = solve(NN_prob NewtonRaphson())

In this case, solve isn't able to find the solution and the following error is returned.

ERROR: SingularException(2)

Was wondering if there would be a way to work around this issue.

Wrap SIAMFANLEquations

https://github.com/ctkelley/SIAMFANLEquations.jl

DF-Sane method

https://www.rdocumentation.org/packages/BB/versions/2019.10-1/topics/dfsane
https://docs.scipy.org/doc/scipy/reference/optimize.root-dfsane.html

J Barzilai, and JM Borwein (1988), Two-point step size gradient methods, IMA J Numerical Analysis, 8, 141-148.

L Grippo, F Lampariello, and S Lucidi (1986), A nonmonotone line search technique for Newton's method, SIAM J on Numerical Analysis, 23, 707-716.

W LaCruz, JM Martinez, and M Raydan (2006), Spectral residual mathod without gradient information for solving large-scale nonlinear systems of equations, Mathematics of Computation, 75, 1429-1448.

R Varadhan and C Roland (2008), Simple and globally-convergent methods for accelerating the convergence of any EM algorithm, Scandinavian J Statistics.

R Varadhan and PD Gilbert (2009), BB: An R Package for Solving a Large System of Nonlinear Equations and for Optimizing a High-Dimensional Nonlinear Objective Function, J. Statistical Software, 32:4, http://www.jstatsoft.org/v32/i04/

warning on 1.7 about vendor() is depreciated

On 1.7 rc1:

┌ Warning: `vendor()` is deprecated, use `BLAS.get_config()` and inspect the output instead
│   caller = (::NonlinearSolve.DefaultLinSolve)(x::Vector{Float64}, A::Matrix{Float64}, b::Vector{Float64}, update_matrix::Bool; tol::Nothing, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}) at utils.jl:125
└ @ NonlinearSolve ~/.julia/packages/NonlinearSolve/37aUT/src/utils.jl:125

Dead links trust region, northwestern.edu

https://docs.sciml.ai/NonlinearSolve/dev/api/nonlinearsolve/
Trust-region methods

Add tests on Regula Falsi

Implement Exciting-Mixing, Scalar and Diagonal Broyden approximations

https://docs.scipy.org/doc/scipy/reference/optimize.root-excitingmixing.html
https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.diagbroyden.html#scipy.optimize.diagbroyden
https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.linearmixing.html#scipy.optimize.linearmixing

Broken link in ReadMe.md

The link to the stable docs in the ReadMe, https://nonlinearsolve.sciml.ai/stable/, is broken (404).

The dev link works, as does https://nonlinearsolve.sciml.ai/.

Newton trust region

Our NewtonRaphson implementation isn't robust against functions that look like this

julia> nl(t, p) = 0.010000000000000002 + 10.000000000000002 / (1 + (0.21640425613334457 + 216.40425613334457 / (1 + (0.21640425613334457 + 216.40425613334457 / (1 + 0.0006250000000000001(t^2.0)))^2.0))^2.0) - 0.0011552453009332421t
nl (generic function with 2 methods)

julia> solve(NonlinearProblem(nl, 10.0), NewtonRaphson()).retcode
:MAXITERS_EXCEED

use LinearSolve.jl for linsolve

https://github.com/SciML/LinearSolve.jl

Can not iterate solver before solving

When attempting to iterate a solver that has not been solved, the user is told that there is no method matching iterate for the solver. Same goes for getindex:

julia> using NonlinearSolve

julia> f(u, p) = u .* u .- 2.0
f (generic function with 1 method)

julia> u0 = (1.0, 2.0) # brackets
(1.0, 2.0)

julia> probB = NonlinearProblem(f, u0)
NonlinearProblem with uType Tuple{Float64, Float64}. In-place: false
u0: (1.0, 2.0)

julia> solver = init(probB, Falsi()) # Can iterate the solver object
NonlinearSolve.BracketingImmutableSolver{NonlinearFunction{false, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing}, Falsi, Float64, Float64, SciMLBase.NullParameters, Nothing, NonlinearProblem{Tuple{Float64, Float64}, false, SciMLBase.NullParameters, NonlinearFunction{false, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}}(1, NonlinearFunction{false, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing}(f, LinearAlgebra.UniformScaling{Bool}(true), nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, SciMLBase.DEFAULT_OBSERVED_NO_TIME, nothing), Falsi(), 1.0, 2.0, -1.0, 2.0, SciMLBase.NullParameters(), false, 1000, NonlinearSolve.DEFAULT, nothing, false, NonlinearProblem{Tuple{Float64, Float64}, false, SciMLBase.NullParameters, NonlinearFunction{false, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}(NonlinearFunction{false, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing}(f, LinearAlgebra.UniformScaling{Bool}(true), nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, SciMLBase.DEFAULT_OBSERVED_NO_TIME, nothing), (1.0, 2.0), SciMLBase.NullParameters(), Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}()))

julia> solver[1]
ERROR: MethodError: no method matching getindex(::NonlinearSolve.BracketingImmutableSolver{NonlinearFunction{false, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing}, Falsi, Float64, Float64, SciMLBase.NullParameters, Nothing, NonlinearProblem{Tuple{Float64, Float64}, false, SciMLBase.NullParameters, NonlinearFunction{false, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}}, ::Int64)
Stacktrace:
 [1] top-level scope
   @ REPL[6]:1

julia> for s in solver
           @show s
       end
ERROR: MethodError: no method matching iterate(::NonlinearSolve.BracketingImmutableSolver{NonlinearFunction{false, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing}, Falsi, Float64, Float64, SciMLBase.NullParameters, Nothing, NonlinearProblem{Tuple{Float64, Float64}, false, SciMLBase.NullParameters, NonlinearFunction{false, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}})
Closest candidates are:
  iterate(::Union{LinRange, StepRangeLen}) at C:\Users\Dennis Bal\.julia\juliaup\julia-1.7.1+0~x64\share\julia\base\range.jl:826
  iterate(::Union{LinRange, StepRangeLen}, ::Integer) at C:\Users\Dennis Bal\.julia\juliaup\julia-1.7.1+0~x64\share\julia\base\range.jl:826
  iterate(::T) where T<:Union{Base.KeySet{<:Any, <:Dict}, Base.ValueIterator{<:Dict}} at C:\Users\Dennis Bal\.julia\juliaup\julia-1.7.1+0~x64\share\julia\base\dict.jl:695
  ...
Stacktrace:
 [1] top-level scope
   @ .\REPL[7]:1

julia> solver = solve!(solver)
u: 1.414213562373095

julia> solver[1]
1.414213562373095

julia> for s in solver
           @show s
       end
s = 1.414213562373095

However, should the unsolved solver not be an empty collection?

This is especially confusing due to the documentation, which says that you can iterate the solver in the line before it is solved. From http://nonlinearsolve.sciml.ai/dev/tutorials/iterator_interface/ :

f(u, p) = u .* u .- 2.0
u0 = (1.0, 2.0) # brackets
probB = NonlinearProblem(f, u0)
solver = init(probB, Falsi()) # Can iterate the solver object
solver = solve!(solver)

Storing trace

I want to store the trace during solving. I would love to use the native NonlinearSolve methods for that, but it appears to me that no such thing is implemented. So I want to use the NLSolveJL interface. But from reading the documentation, I am still not able to see - how do i use it?

So far I have tried the following

#The functional version:
julia> solver = solve(probN, NewtonRaphson(), tol = 1e-9)
u: 8-element Vector{Float64}:
  27.52327682597094
  27.523962383296873
   3.2064613688401113
  13.258386042046885
   0.8000852068277563
   0.4578205135984268
  27.929656903088112
 124.54137097709595

#Failed attempts:
julia> solver = solve(probN, NLSolveJL(method=:newton,), tol = 1e-9)
ERROR: UndefVarError: NLSolveJL not defined
Stacktrace:
 [1] top-level scope
   @ REPL[9]:1

julia> solver = NLSolveJL(probN, tol = 1e-9, method=:newton)        
ERROR: UndefVarError: NLSolveJL not defined
Stacktrace:
 [1] top-level scope
   @ REPL[11]:1

With NLSolveJL not being defined, I was wondering if it is simply not exported. But this seems not to be the case:

julia> NonlinearSolve.NLSolveJL
ERROR: UndefVarError: NLSolveJL not defined
Stacktrace:
 [1] getproperty(x::Module, f::Symbol)
   @ Base .\Base.jl:35
 [2] top-level scope
   @ REPL[12]:1

I think that a section with working examples in the docs of how to access non-native solvers would go a long ways.

Non-allocating stepping for TrustRegion

#116 merged without getting it fully non-allocating, mostly was about finishing correctness. In particular, this line https://github.com/SciML/NonlinearSolve.jl/pull/116/files#diff-d4ded13aeff73c93c2977a4227ff0ba4a434dc462ff8cc433761133aff882405R299 needs a bunch of mul! stuff so that's a good place to start.

[Docs] Stopping conditions

It would be nice to have a section in the documentation for different stopping conditions. The only way I was able to dig the information out was with methods(solve).

Things that should be documented are (possibly incomplete list):

Register

@JuliaRegistrator register()

Double linear solve setups for Levenberg–Marquardt algorithm

There are two separate linear solves, so it would be good to allow linsolve and precs for the two different ones, defaulting to using the same choice for the two. I don't know how helpful this would end up really being in practice, though for iLU you may want to cache different sparsities or something?

Specialize Jacobians when doing Krylov methods to avoid building the matrices

Just doing the AD tricks for Newton-Krylov. Should probably wait for SciMLOperators.jl integration.