Porting from Python optimize.broyden2

General Usage
,
1

Still porting code from Python… What is a good alternative to the optimize.broyden2 solver from SciPy?

https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.broyden2.html

2

FastDEQ.jl has a nice Broyden implementation. @avikpal we really should move it to NonlinearSolve.jl soon :sweat_smile: .

https://fastdeq.sciml.ai/dev/api/nlsolve/

3

For now I just use NLsove:

# from scipy import optimize

# def fun(x):
#     return [x[0]  + 0.5 * (x[0] - x[1])**3 - 1.0,  0.5 * (x[1] - x[0])**3 + x[1]]
    
# def test_broyden2():
#     sol = optimize.broyden2(fun, [0, 0])
#     print(sol)
    
# test_broyden2()
# Output: [0.84116365 0.15883529], needs 10 iterations
# time for executing test_broyden2() loop in µs:  1276

using NLsolve

function fun!(F, x)
    F[1] = x[1]  + 0.5 * (x[1] - x[2])^3 - 1.0
    F[2] = 0.5 * (x[2] - x[1])^3 + x[2]
end

function test_broyden2()
    sol = nlsolve(fun!, [ 0.0; 0.0])
    sol.zero
end
# execution time: 5 µs, 255 times faster than Python

In this example it needs 5 iteration where Python needs 10…

And 255 times faster than Python…

4

And here you go, I made it 50x faster (100 ns)

using NonlinearSolve, StaticArrays, BenchmarkTools

function fun(x, p)
    SA[x[1]  + 0.5 * (x[1] - x[2])^3 - 1.0,
     0.5 * (x[2] - x[1])^3 + x[2]]
end

function test_nr()
    prob = NonlinearProblem{false}(fun,SA[ 0.0; 0.0])
    sol = solve(prob, NewtonRaphson())
    sol.u
end

@btime test_nr() # 106.760 ns (0 allocations: 0 bytes)
2 Likes
5

I already wanted to ask if there is a non-allocating way of doing this. You have been faster than me!

But what about stability? I mean, my real problems are more complex.

6

It’s a Newton-Rhapson, no trust region or line searches in that implementation. But NonlinearSolve is a uniform interface over all of the algorithms we could find, so

https://nonlinearsolve.sciml.ai/dev/solvers/NonlinearSystemSolvers/#SciMLNLSolve.jl

there’s MINPACK and NLsolve in the interface, which could be used until other implementations are added. The optimal algorithm will likely depend on the problem and we still need to do more benchmarking before really knowing the lay of the land here.