IterativeSolvers.jl doesn’t deal with nonlinear equations and NLsolve.jl doesn’t have Krylov methods? What about large sparse nonlinear iterative solvers?
Is your system square? Real or complex? Do you have derivatives?
Square real system. No analytical jacobian.
My package PseudoArcLengthContinuation.jl provides this. Look at the newton function. You have to pass the Jacobian yourself using ForwardDiff or an analytical expression. The docs should give you many examples of how it can be done.
A simple (trivial) example of Matrix-Free is shown here and a more difficult here. You can put IterativeSolvers if you prefer, it is also included in PseudoArcLengthContinuation.jl
2 Likes
If your Jacobian exists and your system is not too large for ForwardDiff, you could try
pkg> add NLPModels, NLPModelsIpopt
julia> using NLPModels, NLPModelsIpopt
julia> F(x) = [x[1] + x[2] - 3; x[1]^2 + x[2]^2 - 9] # your nonlinear system
julia> x0 = [1.0; 5.0] # your initial guess
julia> model = ADNLPModel(x -> 0, x0; c = F, lcon = zero(x0), ucon = zero(x0))
julia> stats = ipopt(model)
julia> stats.solution
2-element Array{Float64,1}:
3.0000000000006
-6.001258810133371e-13
I’ll make a more user-friendly interface, but the above should be equivalent to applying Newton’s method to F(x) = 0.
If your problem is too large for ForwardDiff, you could either supply an analytical Jacobian (which you apparently don’t have), or turn on IPOPT’s quasi-Newton option.
1 Like
If you don’t have an analytical jacobian but you know it is sparse, you should have a look SparseDiffTools.jl. I use it successfully for several challenging problems, it is a great tool (when AD works 
)
1 Like
For Newton-Krylov you should always use ForwardDiff! However, you shouldn’t actually build the Jacobian. In Newton-Krylov you only ever need to do Jv operations, so instead use auto_jacvec
which utilizes dual numbers in a specific way to do J*v computations without building the Jacobian. The cost of doing this is equivalent to calculating only one column of a Jacobian, so it’s much much faster than building Jacobians, and it’s O(1), i.e. it doesn’t matter how big your Jacobian is, it doesn’t scale with the size (only with the cost of f).
Note: you can do J = JacVec(f,u::AbstractArray;autodiff=true) and then put J as “the matrix” into gmres of IterativeSolvers.jl and it’ll “just work”. Then you just have to write the u updates and boom, that’s a Jacobian-Free Newton Krylov implementation.
8 Likes
Can you comment on how to get SparseDiffTools to work for NLsolve, or link an example if one exists? AFAIK, that Newton method needs the whole Jacobian, correct?
Yes, I’ve mentioned this to @pkofod a few times. NLSolvers is undergoing an overhaul right now which might add this capability.
1 Like
When you add trust-region modifications to Newton-Krylov, you have to constrain the trust-region optimization problem to the Krylov subspace. It’s just a bit of extra low-d linear algebra, but it ends up pushing the need for knowledge of the Krylov subspace into the trust-region calculation. You can’t just stick a trust-region algorithm and a Krylov linear algebra solver together and have them “just work.” So I expect a library like NLSolve will need some modification to work with Krylov methods. Unless of course NLSolve has already has that Krylov tweak incorporated.
And revelt’s PseudoArclengthContinuation.jl has already put these algorithms together.
2 Likes
If you don’t want to form the Jacobian, you can modify my previous suggestion as
pkg> add JSOSolvers
julia> nlsmodel = FeasibilityResidual(model)
julia> stats = trunk(nlsmodel)
julia> stats.solution
2-element Array{Float64,1}:
-7.03336155582655e-8
3.000000070333616
This approach solves the nonlinear system as a nonlinear least-squares problem using a variant of the Levenberg-Marquardt method.
By default, it will use LSMR to solve the trust-region subproblem, which should be less memory-hungry than GMRES, even though your system is square. This variant will only use Jacobian-vector and transposed-Jacobian-vector products computed with ForwardDiff without forming the Jacobian.
Any solver that calls the jac() method (such as IPOPT) will form the Jacobian.
This (no need for the Jacobian, only Jv) is very common in various algorithms so I am wondering if it would make sense to develop some conventions for it. Some of the Krylov libraries already work like this.
Eg J is an object that supports
AbstractMatrix(J), all methods below can fall back to thisJ * v, returning a vector of the same size and type asv,J' * v, similar- maybe
J' * J?
Yes, I’ve mentioned this to @pkofod a few times.
I guess it is mostly there. You can pass the linearsolver with the argument linsolve
Does it make sense to either change the name to be less specific or refactor some of that stuff to shared libraries? I thought the library just has continuation solver algorithms used in dynamical systems, so missed this stuff.
3 Likes
Yeah, but I don’t think it gets rid of the Jacobian calculation when you do that, so we would need a wide to get rid of it.
1 Like
Maybe LargeScaleBifurcations.jl as it is oriented towards large systems, unlike Bifurcations.jl . If you think of a better name, I am all ears !!
Also, PseudoArclengthContinuation.jl is annoying as I want to push another type of continuation which is not arclength based.
That is easier to find. But if there are a huge number of features that have nothing to do with dynamical systems, even that is too specific. But perhaps that is a better name and then you could refactor out the general code into its ownpackage in the future.
ThePackageFormerlyKnownAsPseudoArclengthContinuation.jl is always an option 
It would also be a nice name for a rock band.
2 Likes
That’s a cool package! Surprised to know there are those options. I’m mainly interested in solving steady state nonlinear DEs. Would be curious to see how it compares with the dynamic solvers in DifferentialEquations. Nlsolve.jl is not really useful for large systems.
Thank you!
Would be curious to see how it compares with the dynamic solvers in DifferentialEquations.
I have not done dynamics solvers per se, nothing like time steppers. DE is too good for that.
The closest I can think of are methods for the computation of periodic orbits as BVP. I do not use the solvers of DE for that, except for the shooting algorithm where I call DE. These BVPs are quite specific (you have to add a constraint otherwise they are not well posed) so that’s why I dont call DE. I also use a lot of tricks to make the computation faster and the solvers are optimized for large scale problems. When the periodic orbit has a large period (close to homoclinic point), the methods show their weaknesses.
Finally, the interface is less polished than for DE because I want to leave the user the possibility for not using DE. Also, I do not provide AD in PseudoArclengthContinuation.jl because the package are evolving too fast. This could be improved by a small package though but I have no use for that.
Nlsolve.jl is not really useful for large systems.
I think it is not far, you have to use the argument linsolve. If you look at this test example, you’ll get ideas. There is surely a way to pass a Matrix-Free solver in it.
You can use my newton function until they improve NLsolve.jl, especially the linear solvers part.
There is also a good newton solver in DE but it is a bit buried in the code.
In the near future, you will be able to pass your own newton solver in PseudoArclengthContinuation.jl bypassing the one used for now.
I will push a set of improvements (mainly for periodic orbits) and bug fixed for PseudoArclengthContinuation.jl next week.
1 Like