# Nonlinear equations on abstract vectors?

**URL:** <https://discourse.julialang.org/t/nonlinear-equations-on-abstract-vectors/28692>\
**Category:** Optimization (Mathematical)\
**Created:** [September 12, 2019, 2:47pm UTC](https://discourse.julialang.org/t/nonlinear-equations-on-abstract-vectors/28692 "2019-09-12T14:47:15Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![tobydriscoll](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tobydriscoll/32/1843_2.png) [@tobydriscoll](https://discourse.julialang.org/u/tobydriscoll)\
**Post date:** [September 12, 2019, 2:47pm UTC](https://discourse.julialang.org/t/nonlinear-equations-on-abstract-vectors/28692/1 "2019-09-12T14:47:15Z")

</div>

I need to solve a (smooth, unconstrained) system of nonlinear equations for an unknown defined on multiple subdomains (to wit, domain decomposition for a PDE). A natural representation of the unknown might be a vector of arrays of different sizes.

My first thought was that I would need to give a true Vector to the nonlinear solver and then view it as the subdomain-aware version for my stuff. Something like a PseudoBlockVector with a resize operation within each subdomain. Doable, but I never look forward to making multiple views into the same data without creating unwanted copies (or failing to create needed copies, or deep copies…).

Then I thought a cooler Julian solution would be to subtype AbstractVector for my representation and generalize +,-,length, and so on, so that the nonlinear solver would never have to know it wasn’t working with a true Vector.

Does anyone know if the NLsolve or Optim packages, or another not on my radar, can handle that? And what the “and so on” set of necessary vector operations would be?

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [September 12, 2019, 3:30pm UTC](https://discourse.julialang.org/t/nonlinear-equations-on-abstract-vectors/28692/2 "2019-09-12T15:30:34Z")

</div>

> [@tobydriscoll](#):
>
> Then I thought a cooler Julian solution would be to subtype AbstractVector for my representation and generalize +,-,length, and so on, so that the nonlinear solver would never have to know it wasn’t working with a true Vector.

Doesn’t `BlockArray` from [GitHub - JuliaArrays/BlockArrays.jl: BlockArrays for Julia](https://github.com/JuliaArrays/BlockArrays.jl) do what you want?

---

<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:** [September 12, 2019, 3:46pm UTC](https://discourse.julialang.org/t/nonlinear-equations-on-abstract-vectors/28692/3 "2019-09-12T15:46:30Z")

</div>

You don’t need to make it an `AbstractVector`, only a subset of the methods is required for a `newton` solve. For example in my package [PseudoArcLengthContinuation.jl](https://github.com/rveltz/PseudoArcLengthContinuation.jl), I am relying on the same set of methods as [KrylovKit.jl](https://github.com/Jutho/KrylovKit.jl). It provides `newton` method and others for ``general’’ vectors.

You can also provide you own linear solver if you wish.

---

<div class="post-metadata">

**Author:** ![tobydriscoll](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tobydriscoll/32/1843_2.png) [@tobydriscoll](https://discourse.julialang.org/u/tobydriscoll)\
**Post date:** [September 12, 2019, 5:32pm UTC](https://discourse.julialang.org/t/nonlinear-equations-on-abstract-vectors/28692/4 "2019-09-12T17:32:00Z")

</div>

Again, that was my first thought, and it should be fine to start. Eventually, though, I might want different parts of the vector to live on different processes. Wouldn’t that leave me stuck if the solver wants to, say, add vectors together, unless it is willing to accept an `AbstractVector` in the first place?

---

<div class="post-metadata">

**Author:** ![tobydriscoll](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tobydriscoll/32/1843_2.png) [@tobydriscoll](https://discourse.julialang.org/u/tobydriscoll)\
**Post date:** [September 12, 2019, 6:31pm UTC](https://discourse.julialang.org/t/nonlinear-equations-on-abstract-vectors/28692/5 "2019-09-12T18:31:03Z")

</div>

@rveltz, it looks like [KrlovKit’s description of vector-like behavior](https://jutho.github.io/KrylovKit.jl/stable/) is the kind of thing I’m after. Your PALC might fit the bill, then. It looks like nice work! I must say that after reading your front-page disclaimer about performance, I was very surprised by the extent to which you seem to defy it in your examples.

---

<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:** [September 12, 2019, 7:02pm UTC](https://discourse.julialang.org/t/nonlinear-equations-on-abstract-vectors/28692/6 "2019-09-12T19:02:07Z")

</div>

It could be but I would not worry too much about this. For example: see `DArray{T,N,A} <: AbstractArray{T,N}` in [DistributedArrays.jl](https://github.com/JuliaParallel/DistributedArrays.jl). See also (WIP) [MPIArray.jl](https://github.com/barche/MPIArrays.jl/blob/master/src/MPIArrays.jl#L23) which gives an `AbstractArray`. Knowing a bit where you’re headed would help you: for example, do you want cluster of GPUs, cluster of CPU, multithreading…

For example, using `AbstractArray` allows some performance improvement by calling BLAS (search in `KrylovKit` about this).

> You might be interested in [RecursiveArrayTools.jl](https://github.com/JuliaDiffEq/RecursiveArrayTools.jl) for your partition problem

When I first wrote [PseudoArcLengthContinuation.jl](https://github.com/rveltz/PseudoArcLengthContinuation.jl), I was using `ArrayFire.jl`, `KrylovKit.jl` and `ApproxFun.jl`, that’s why the package does not require `AbstractArray`.  
I don’t use much [IterativeSolvers.jl](https://github.com/JuliaMath/IterativeSolvers.jl) because it lacks an eigensolver, it is otherwise a wonderful package.

---

<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:** [September 12, 2019, 7:39pm UTC](https://discourse.julialang.org/t/nonlinear-equations-on-abstract-vectors/28692/7 "2019-09-12T19:39:01Z")

</div>

> [@tobydriscoll](#):
>
> It looks like nice work!

Thank you!

> [@tobydriscoll](#):
>
> I was very surprised by the extent to which you seem to defy it in your examples.

Well, I do researching with this tool studying some challenging PDEs… So I guess I toned down too much the disclaimer about performance. Most of the good performance is squeezed out of the linear solver though…

That said, @ChrisRackauckas explains very well what such type of ``simple’’ package performs well in this [post](http://www.stochasticlifestyle.com/when-do-micro-optimizations-matter-in-scientific-computing/). Also, if you are interested in PALC, you should check the excellent [Bifurcations.jl](https://github.com/tkf/Bifurcations.jl)

---

<div class="post-metadata">

**Author:** ![juthohaegeman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juthohaegeman/32/8620_2.png) [@juthohaegeman](https://discourse.julialang.org/u/juthohaegeman)\
**Post date:** [September 12, 2019, 8:12pm UTC](https://discourse.julialang.org/t/nonlinear-equations-on-abstract-vectors/28692/8 "2019-09-12T20:12:51Z")

</div>

KrylovKit.jl is for linear problems and eigenvalue problems, and related iterative methods. For non-linear optimization problems, I was also working on [OptimKit.jl](https://github.com/Jutho/OptimKit.jl) , but it’s not registered and still work in progress (though the basic functionality is working). It goes a bit further than just relying you to overload `dot` and `axpy!` and these kind of functions for your custom vector-like type (which does not need to be `<:AbstractVector` or `<:AbstractArray`). You can just specify which method it has to use for those so your function for adding ‘vectors’ can even have a different name, and you don’t need to go through the trouble of defining a new custom type. If the parameters of your optimization problem are best captured as say a tuple of a scalar, a matrix, and some other thing, you can just use that, i.e. `x=(scalar, matrix, otherobj)` and then define functions for adding and scaling these, and for computing inner products.

---

<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:** [September 12, 2019, 8:52pm UTC](https://discourse.julialang.org/t/nonlinear-equations-on-abstract-vectors/28692/9 "2019-09-12T20:52:23Z")

</div>

We should tease you a bit more so we can get cool new packages more often 😉