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 MatrixFree is shown here and a more difficult here. You can put IterativeSolvers if you prefer, it is also included in PseudoArcLengthContinuation.jl
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
2element Array{Float64,1}:
3.0000000000006
6.001258810133371e13
I’ll make a more userfriendly 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 quasiNewton option.
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 )
For NewtonKrylov you should always use ForwardDiff! However, you shouldn’t actually build the Jacobian. In NewtonKrylov 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 JacobianFree Newton Krylov implementation.
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.
When you add trustregion modifications to NewtonKrylov, you have to constrain the trustregion optimization problem to the Krylov subspace. It’s just a bit of extra lowd linear algebra, but it ends up pushing the need for knowledge of the Krylov subspace into the trustregion calculation. You can’t just stick a trustregion 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.
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
2element Array{Float64,1}:
7.03336155582655e8
3.000000070333616
This approach solves the nonlinear system as a nonlinear leastsquares problem using a variant of the LevenbergMarquardt method.
By default, it will use LSMR to solve the trustregion subproblem, which should be less memoryhungry than GMRES, even though your system is square. This variant will only use Jacobianvector and transposedJacobianvector 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 this 
J * 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.
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.
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.
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 MatrixFree 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.