# Is there a julia package that uses Newton-Krylov method for systems of nonlinear equations?

**URL:** <https://discourse.julialang.org/t/is-there-a-julia-package-that-uses-newton-krylov-method-for-systems-of-nonlinear-equations/36520>\
**Category:** General Usage\
**Created:** [March 25, 2020, 9:25pm UTC](https://discourse.julialang.org/t/is-there-a-julia-package-that-uses-newton-krylov-method-for-systems-of-nonlinear-equations/36520 "2020-03-25T21:25:20Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![JianghuiDu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jianghuidu/32/8028_2.png) [@JianghuiDu](https://discourse.julialang.org/u/JianghuiDu)\
**Post date:** [March 25, 2020, 9:25pm UTC](https://discourse.julialang.org/t/is-there-a-julia-package-that-uses-newton-krylov-method-for-systems-of-nonlinear-equations/36520/1 "2020-03-25T21:25:20Z")

</div>

IterativeSolvers.jl doesn’t deal with nonlinear equations and NLsolve.jl doesn’t have Krylov methods? What about large sparse nonlinear iterative solvers?

---

<div class="post-metadata">

**Author:** ![dpo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpo/32/3335_2.png) [@dpo](https://discourse.julialang.org/u/dpo)\
**Post date:** [March 26, 2020, 4:55am UTC](https://discourse.julialang.org/t/is-there-a-julia-package-that-uses-newton-krylov-method-for-systems-of-nonlinear-equations/36520/2 "2020-03-26T04:55:11Z")

</div>

Is your system square? Real or complex? Do you have derivatives?

---

<div class="post-metadata">

**Author:** ![JianghuiDu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jianghuidu/32/8028_2.png) [@JianghuiDu](https://discourse.julialang.org/u/JianghuiDu)\
**Post date:** [March 26, 2020, 9:55am UTC](https://discourse.julialang.org/t/is-there-a-julia-package-that-uses-newton-krylov-method-for-systems-of-nonlinear-equations/36520/3 "2020-03-26T09:55:03Z")

</div>

Square real system. No analytical jacobian.

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [March 26, 2020, 10:00am UTC](https://discourse.julialang.org/t/is-there-a-julia-package-that-uses-newton-krylov-method-for-systems-of-nonlinear-equations/36520/4 "2020-03-26T10:00:05Z")

</div>

My package [PseudoArcLengthContinuation.jl](https://www.github.com/rveltz/PseudoArcLengthContinuation.jl) provides this. Look at the `newton` function. You have to pass the Jacobian yourself using `ForwardDiff` or an analytical expression. The [docs](https://rveltz.github.io/PseudoArcLengthContinuation.jl/dev/) should give you many examples of how it can be done.

A simple (trivial) example of Matrix-Free is shown [here](https://rveltz.github.io/PseudoArcLengthContinuation.jl/dev/tutorials1/#Using-GMRES-or-another-linear-solver-1) and a more difficult [here](https://rveltz.github.io/PseudoArcLengthContinuation.jl/dev/tutorials2b/). You can put IterativeSolvers if you prefer, it is also included in [PseudoArcLengthContinuation.jl](https://www.github.com/rveltz/PseudoArcLengthContinuation.jl)

---

<div class="post-metadata">

**Author:** ![dpo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpo/32/3335_2.png) [@dpo](https://discourse.julialang.org/u/dpo)\
**Post date:** [March 26, 2020, 2:09pm UTC](https://discourse.julialang.org/t/is-there-a-julia-package-that-uses-newton-krylov-method-for-systems-of-nonlinear-equations/36520/5 "2020-03-26T14:09:25Z")

</div>

If your Jacobian exists and your system is not too large for ForwardDiff, you could try

```julia
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.

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [March 26, 2020, 2:24pm UTC](https://discourse.julialang.org/t/is-there-a-julia-package-that-uses-newton-krylov-method-for-systems-of-nonlinear-equations/36520/6 "2020-03-26T14:24:02Z")

</div>

If you don’t have an analytical jacobian but you know it is sparse, you should have a look [SparseDiffTools.jl](https://github.com/JuliaDiff/SparseDiffTools.jl). I use it successfully for several challenging problems, it is a great tool (when AD works 🤣)

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [March 26, 2020, 4:04pm UTC](https://discourse.julialang.org/t/is-there-a-julia-package-that-uses-newton-krylov-method-for-systems-of-nonlinear-equations/36520/7 "2020-03-26T16:04:48Z")

</div>

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`

> **[GitHub - JuliaDiff/SparseDiffTools.jl: Fast jacobian computation through...](https://github.com/JuliaDiff/SparseDiffTools.jl#jacobian-vector-and-hessian-vector-products)**
>
> Fast jacobian computation through sparsity exploitation and matrix coloring - GitHub - JuliaDiff/SparseDiffTools.jl: Fast jacobian computation through sparsity exploitation and matrix coloring

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.

---

<div class="post-metadata">

**Author:** ![platawiec](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/platawiec/32/31914_2.png) [@platawiec](https://discourse.julialang.org/u/platawiec)\
**Post date:** [March 26, 2020, 5:05pm UTC](https://discourse.julialang.org/t/is-there-a-julia-package-that-uses-newton-krylov-method-for-systems-of-nonlinear-equations/36520/8 "2020-03-26T17:05:19Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [March 26, 2020, 5:11pm UTC](https://discourse.julialang.org/t/is-there-a-julia-package-that-uses-newton-krylov-method-for-systems-of-nonlinear-equations/36520/9 "2020-03-26T17:11:25Z")

</div>

Yes, I’ve mentioned this to @pkofod a few times. NLSolvers is undergoing an overhaul right now which might add this capability.

---

<div class="post-metadata">

**Author:** ![John\_Gibson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/john_gibson/32/5321_2.png) [@John\_Gibson](https://discourse.julialang.org/u/John_Gibson)\
**Post date:** [March 26, 2020, 5:27pm UTC](https://discourse.julialang.org/t/is-there-a-julia-package-that-uses-newton-krylov-method-for-systems-of-nonlinear-equations/36520/10 "2020-03-26T17:27:56Z")

</div>

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](https://www.github.com/rveltz/PseudoArcLengthContinuation.jl) has already put these algorithms together.

---

<div class="post-metadata">

**Author:** ![dpo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpo/32/3335_2.png) [@dpo](https://discourse.julialang.org/u/dpo)\
**Post date:** [March 26, 2020, 7:38pm UTC](https://discourse.julialang.org/t/is-there-a-julia-package-that-uses-newton-krylov-method-for-systems-of-nonlinear-equations/36520/11 "2020-03-26T19:38:40Z")

</div>

If you don’t want to form the Jacobian, you can modify my previous suggestion as

```julia
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](https://juliasmoothoptimizers.github.io/Krylov.jl/dev/solvers/#Krylov.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()`](https://juliasmoothoptimizers.github.io/NLPModels.jl/stable/api/#NLPModels.jac) method (such as IPOPT) _will_ form the Jacobian.

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [March 27, 2020, 6:35am UTC](https://discourse.julialang.org/t/is-there-a-julia-package-that-uses-newton-krylov-method-for-systems-of-nonlinear-equations/36520/12 "2020-03-27T06:35:51Z")

</div>

> [@ChrisRackauckas](#):
>
> you shouldn’t actually build the Jacobian. In Newton-Krylov you only ever need to do `Jv` operations

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 as `v`,
- `J' * v`, similar
- maybe `J' * J`?

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [March 27, 2020, 12:27pm UTC](https://discourse.julialang.org/t/is-there-a-julia-package-that-uses-newton-krylov-method-for-systems-of-nonlinear-equations/36520/13 "2020-03-27T12:27:20Z")

</div>

> 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`

---

<div class="post-metadata">

**Author:** ![jlperla](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jlperla/32/34332_2.png) [@jlperla](https://discourse.julialang.org/u/jlperla)\
**Post date:** [March 27, 2020, 1:58pm UTC](https://discourse.julialang.org/t/is-there-a-julia-package-that-uses-newton-krylov-method-for-systems-of-nonlinear-equations/36520/14 "2020-03-27T13:58:23Z")

</div>

> [@John\_Gibson](#):
>
> And revelt’s [PseudoArclengthContinuation.jl](https://www.github.com/rveltz/PseudoArcLengthContinuation.jl) has already put these algorithms together.

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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [March 27, 2020, 2:30pm UTC](https://discourse.julialang.org/t/is-there-a-julia-package-that-uses-newton-krylov-method-for-systems-of-nonlinear-equations/36520/15 "2020-03-27T14:30:15Z")

</div>

> [@rveltz](#):
>
> I guess it is mostly there. You can pass the linearsolver with the argument `linsolve`

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.

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [March 27, 2020, 2:54pm UTC](https://discourse.julialang.org/t/is-there-a-julia-package-that-uses-newton-krylov-method-for-systems-of-nonlinear-equations/36520/16 "2020-03-27T14:54:48Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![jlperla](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jlperla/32/34332_2.png) [@jlperla](https://discourse.julialang.org/u/jlperla)\
**Post date:** [March 27, 2020, 3:11pm UTC](https://discourse.julialang.org/t/is-there-a-julia-package-that-uses-newton-krylov-method-for-systems-of-nonlinear-equations/36520/17 "2020-03-27T15:11:06Z")

</div>

> [@rveltz](#):
>
> Maybe LargeScaleBifurcations.jl as it is oriented towards large systems, unlike Bifurcations.jl . If you think of a better name, I am all ears !!

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.

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [March 28, 2020, 3:28pm UTC](https://discourse.julialang.org/t/is-there-a-julia-package-that-uses-newton-krylov-method-for-systems-of-nonlinear-equations/36520/18 "2020-03-28T15:28:38Z")

</div>

> [@rveltz](#):
>
> Also, PseudoArclengthContinuation.jl is annoying as I want to push another type of continuation which is not arclength based.

ThePackageFormerlyKnownAsPseudoArclengthContinuation.jl is always an option 😉

It would also be a nice name for a rock band.

---

<div class="post-metadata">

**Author:** ![JianghuiDu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jianghuidu/32/8028_2.png) [@JianghuiDu](https://discourse.julialang.org/u/JianghuiDu)\
**Post date:** [March 29, 2020, 3:21pm UTC](https://discourse.julialang.org/t/is-there-a-julia-package-that-uses-newton-krylov-method-for-systems-of-nonlinear-equations/36520/19 "2020-03-29T15:21:29Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [March 29, 2020, 5:43pm UTC](https://discourse.julialang.org/t/is-there-a-julia-package-that-uses-newton-krylov-method-for-systems-of-nonlinear-equations/36520/20 "2020-03-29T17:43:03Z")

</div>

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](https://github.com/JuliaNLSolvers/NLsolve.jl/blob/master/test/linsolve.jl), 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._

[Next page](https://discourse.julialang.org/t/is-there-a-julia-package-that-uses-newton-krylov-method-for-systems-of-nonlinear-equations/36520.md?page=2)