# Status of integration, optimization, and nl solvers
**URL:** https://discourse.julialang.org/t/status-of-integration-optimization-and-nl-solvers/60025
**Category:** Numerics
**Tags:** question
**Created:** [April 26, 2021, 10:55am UTC](https://discourse.julialang.org/t/status-of-integration-optimization-and-nl-solvers/60025 "2021-04-26T10:55:07Z")
**Posts on this page:** 20
**Page:** 1
### Author:  [@amrods](https://discourse.julialang.org/u/amrods)
#### Post date: [April 26, 2021, 10:55am UTC](https://discourse.julialang.org/t/status-of-integration-optimization-and-nl-solvers/60025/1 "2021-04-26T10:55:07Z")
I’m working on a project in which, given a vector of parameters `θ`, I have to solve a nonlinear system of equations that involve some integrals. After solving the system, I need to optimize with respect to `θ` to obtain optimal parameter values that fit some data. That means that I have to solve that system of equations many times.
I have not used any packages for integration, but I have used `Optim.jl` and `NLsolve.jl` in the past. However it’s been a while since I last worked with those packages and I’m wondering if they are still the best in Julia. Also, I’m looking for suggestions of what other packages I should look at. I know that there is `GalacticOptim.jl` and `Quadrature.jl` that wrap many other optimization and integration packages, any others? And what about nl solvers? What packages should I look at that play nice with `Optim` or `GalacticOptim`, `Quadrature` and `NLsolve`?
---
### Author:  [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [April 26, 2021, 11:03am UTC](https://discourse.julialang.org/t/status-of-integration-optimization-and-nl-solvers/60025/2 "2021-04-26T11:03:54Z")
Depending on your integrals (dimension, smoothness, etc), consider approximating with a Gaussian quadrature, eg as in
> **[GitHub - JuliaApproximation/FastGaussQuadrature.jl: Julia package for Gaussian quadrature](https://github.com/JuliaApproximation/FastGaussQuadrature.jl)**
>
> Julia package for Gaussian quadrature
and from then on you effectively have a sum of function values with (potentially) parameter-dependent weights, in any case it should be easy to handle with any nonlinear solver.
If this is about structural estimation, I recently asked a similar question
> [@Optimization on a manifold](https://discourse.julialang.org/t/optimization-on-a-manifold/59850):
>
> This is both a conceptual question, and a practical one of what packages I could use to solve this. Problem statement I am trying to solve the following problem. Let a \in \mathbb{R}^M be a vector of parameters, and p \in \mathbb{R}^N prices. Market clearing F: \mathbb{R}^M \times \mathbb{R}^N \to \mathbb{R}^N requires that F(a, p) = 0 where F is differentiable in both arguments (both in theory and practice). Take a leap of faith and assume that for each a, there is a unique p(a). I am minim…
and got some good suggestions, which I am working through. Currently I find NLopt (Augmented Lagrangian for the constraint) very robust, but that’s on a toy problem.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [April 26, 2021, 11:18am UTC](https://discourse.julialang.org/t/status-of-integration-optimization-and-nl-solvers/60025/3 "2021-04-26T11:18:10Z")
> [@amrods](#):
>
> I have not used any packages for integration, but I have used `Optim.jl` and `NLsolve.jl` in the past. However it’s been a while since I last worked with those packages and I’m wondering if they are still the best in Julia. Also, I’m looking for suggestions of what other packages I should look at. I know that there is `GalacticOptim.jl` and `Quadrature.jl` that wrap many other optimization and integration packages, any others? And what about nl solvers? What packages should I look at that play nice with `Optim` or `GalacticOptim` , `Quadrature` and `NLsolve` ?
Quadrature.jl is compatible with automatic differentiation, so it should play nice with the whole ML ecosystem, etc. It just composes. See for example this example of optimization under uncertainty, mixing autodiff, Quadrature.jl, DifferentialEquations.jl, event handling, multithreading, and a bit more:
[http://tutorials.juliadiffeq.org/html/DiffEqUncertainty/02-AD\_and\_optimization.html](http://tutorials.juliadiffeq.org/html/DiffEqUncertainty/02-AD_and_optimization.html)
For a bunch of examples on ML integration, see:
[https://diffeqflux.sciml.ai/dev/](https://diffeqflux.sciml.ai/dev/)
---
### Author:  [@amrods](https://discourse.julialang.org/u/amrods)
#### Post date: [April 26, 2021, 11:21am UTC](https://discourse.julialang.org/t/status-of-integration-optimization-and-nl-solvers/60025/4 "2021-04-26T11:21:52Z")
> [@Tamas\_Papp](#):
>
> Depending on your integrals (dimension, smoothness, etc), consider approximating with a Gaussian quadrature
The integrals are fairly simple: \int\_0^T K e^{-\rho s} ds and \int\_T^\infty H(s, x) e^{-\rho s} ds. I can pick the functional form of H, but probably it will be of the form Z(s)\phi(x). I suppose these are simple enough for `FastGaussianQuadrature.jl`?
> [@Tamas\_Papp](#):
>
> If this is about structural estimation, I recently asked a similar question
It is about structural estimation, but not as general as your problem. It’s not a GE model, but rather a duration model with behavioral interpretation.
> [@Tamas\_Papp](#):
>
> Currently I find NLopt (Augmented Lagrangian for the constraint) very robust, but that’s on a toy problem.
That seems to call a library outside of Julia. I suppose that automatic differentiation does not work here. Also, last time I worked with `Optim`, it was not easily parallelizable, what about `NLopt`?
---
### Author:  [@amrods](https://discourse.julialang.org/u/amrods)
#### Post date: [April 26, 2021, 11:27am UTC](https://discourse.julialang.org/t/status-of-integration-optimization-and-nl-solvers/60025/5 "2021-04-26T11:27:08Z")
I’ll take a look at that, thank you.
---
### Author:  [@amrods](https://discourse.julialang.org/u/amrods)
#### Post date: [April 26, 2021, 11:34am UTC](https://discourse.julialang.org/t/status-of-integration-optimization-and-nl-solvers/60025/6 "2021-04-26T11:34:28Z")
> [@amrods](#):
>
> I suppose these are simple enough for `FastGaussianQuadrature.jl` ?
I see that `FastGaussianQuadrature` integrates over a closed interval, I guess `Quadrature` it is then.
---
### Author:  [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [April 26, 2021, 11:35am UTC](https://discourse.julialang.org/t/status-of-integration-optimization-and-nl-solvers/60025/7 "2021-04-26T11:35:41Z")
> [@amrods](#):
>
> I suppose these are simple enough for `FastGaussianQuadrature.jl` ?
Yes, with (shifted) Laguerre polynomials.
> [@amrods](#):
>
> That seems to call a library outside of Julia. I suppose that automatic differentiation does not work here.
You can use AD to supply the derivative, it works fine.
> [@amrods](#):
>
> last time I worked with `Optim` , it was not easily parallelizable, what about `NLopt` ?
A lot of local algorithms are sequential and not inherently parallelizable, but you can make your objective parallelizable just fine.
---
### Author:  [@amrods](https://discourse.julialang.org/u/amrods)
#### Post date: [April 26, 2021, 11:46am UTC](https://discourse.julialang.org/t/status-of-integration-optimization-and-nl-solvers/60025/8 "2021-04-26T11:46:15Z")
I’d like to play with that. I think it could be fastest way to integrate my functions. Can you provide a reference? Wikipedia doesn’t go in depth for the use of Laguerre polynomials in integrals.
EDIT: Found it in Judd, but if you have another reference I’d appreciate it.
---
### Author:  [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [April 26, 2021, 11:57am UTC](https://discourse.julialang.org/t/status-of-integration-optimization-and-nl-solvers/60025/9 "2021-04-26T11:57:54Z")
If you have the Judd book I think it is the best reference, really. Miranda and Fackler has a short subsection but few details.
---
### Author:  [@jlperla](https://discourse.julialang.org/u/jlperla)
#### Post date: [April 26, 2021, 2:31pm UTC](https://discourse.julialang.org/t/status-of-integration-optimization-and-nl-solvers/60025/10 "2021-04-26T14:31:28Z")
[GitHub - QuantEcon/Expectations.jl: Expectation operators for Distributions.jl objects](https://github.com/QuantEcon/Expectations.jl) is largely intended to wrap the (non-adaptive) Gaussian quadrature routines in FastGaussQuadrature.
Depending on your use case, it may be the most convenient (and all of the change-of-variables required to transform the spaces is implemented).
---
### Author:  [@amrods](https://discourse.julialang.org/u/amrods)
#### Post date: [April 29, 2021, 11:28pm UTC](https://discourse.julialang.org/t/status-of-integration-optimization-and-nl-solvers/60025/11 "2021-04-29T23:28:32Z")
Do you use `NLsolve` directly or through `GalacticOptim`?
---
### Author:  [@longemen3000](https://discourse.julialang.org/u/longemen3000)
#### Post date: [April 30, 2021, 2:13am UTC](https://discourse.julialang.org/t/status-of-integration-optimization-and-nl-solvers/60025/12 "2021-04-30T02:13:06Z")
maybe not ready for prime time, but i have been using [GitHub - JuliaNLSolvers/NLSolvers.jl: No bells and whistles foundation of Optim.jl](https://github.com/JuliaNLSolvers/NLSolvers.jl) as a replacement for both Optim (unconstrained) and NLSolve
---
### Author:  [@Albert\_Zevelev](https://discourse.julialang.org/u/Albert_Zevelev)
#### Post date: [April 30, 2021, 3:37am UTC](https://discourse.julialang.org/t/status-of-integration-optimization-and-nl-solvers/60025/13 "2021-04-30T03:37:26Z")
It looks like your solving an econ continuous time problem.
I’ve thought about doing this in julia quite a bit.
Can you please write out your problem?
---
### Author:  [@amrods](https://discourse.julialang.org/u/amrods)
#### Post date: [April 30, 2021, 4:08am UTC](https://discourse.julialang.org/t/status-of-integration-optimization-and-nl-solvers/60025/14 "2021-04-30T04:08:50Z")
It’s a Nash bargaining problem of time of retirement for couples (Honoré & de Paula 2018). The optimization problem is
\max \_{t\_{1}, t\_{2}}\left(\int\_{0}^{t\_{1}} K\_{1} e^{-\rho s} d s+\int\_{t\_{1}}^{\infty} H\_{1}\left(s, {x}\_{1}\right) D\left(s, t\_{2}\right) e^{-\rho s} d s-A\_{1}\right) \\ \quad \times\left(\int\_{0}^{t\_{2}} K\_{2} e^{-\rho s} d s+\int\_{t\_{2}}^{\infty} H\_{2}\left(s, {x}\_{2}\right) D\left(s, t\_{1}\right) e^{-\rho s} d s-A\_{2}\right)
where D\left(s, t\_{j}\right)=(\delta-1) \mathbb{1}\left(s \geq t\_{j}\right)+1. The objective function simplifies to
\begin{aligned} N\left(t\_{1}, t\_{2}\right)=& \overbrace{\left(K\_{1} \rho^{-1}\left(1-e^{-\rho t\_{1}}\right)+\widetilde{H}\_{1}\left(t\_{1}, {x}\_{i}\right)+(\delta-1) \widetilde{H}\_{1}\left(\max \left\{t\_{1}, t\_{2}\right\}, {x}\_{1}\right)-A\_{1}\right)}^{\equiv I} \\ & \times \underbrace{\left(K\_{2} \rho^{-1}\left(1-e^{-\rho t\_{2}}\right)+\widetilde{H}\_{2}\left(t\_{2}, {x}\_{2}\right)+(\delta-1) \widetilde{H}\_{2}\left(\max \left\{t\_{1}, t\_{2}\right\}, {x}\_{2}\right)-A\_{2}\right)}\_{\equiv I I} \end{aligned}
with \widetilde{H}\_{i}\left(t, {x}\_{i}\right)=\int\_{t}^{\infty} H\_{i}\left(s, {x}\_{i}\right) e^{-\rho s} d s.
There are 3 cases of FOC to maximize that objective function: when t^\*\_1 \< t^\*\_2, when t^\*\_1 \> t^\*\_2 and when t^\*\_1 = t^\*\_2. There are known bounds on t^\*\_1, t^\*\_2 for each of those cases.
In [another post](https://discourse.julialang.org/t/how-to-pass-constraints-in-nlsolve/60219) I ask about how to solve those FOC with `NLsolve` by supplying those theoretical bounds on the optimal values t^\*\_1, t^\*\_2. I’d appreciate any suggestion you may have.
---
### Author:  [@amrods](https://discourse.julialang.org/u/amrods)
#### Post date: [April 30, 2021, 10:50am UTC](https://discourse.julialang.org/t/status-of-integration-optimization-and-nl-solvers/60025/15 "2021-04-30T10:50:59Z")
Actually, the solution for t\_1 is independent of the solution for t\_2, so I can use the (faster?) univariate methods in `Roots.jl`. Then I can find t^\*\_2 with the same univariate approach.
---
### Author:  [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [April 30, 2021, 12:01pm UTC](https://discourse.julialang.org/t/status-of-integration-optimization-and-nl-solvers/60025/16 "2021-04-30T12:01:29Z")
\tilde{H} should be an ideal candidate for quadrature with fixed nodes (just shift by t).
Personally would just write a callable that evaluates that,
```julia
struct Hquad{T,V<:AbstractVector{T}}
nodes::V
weights::V
ρ::T
end
function (hq::Hquad)(t, x)
...
end
```
and unit test it (because I have a tendency to mess up simple transformations like adding t and dividing by \rho), and then it should be straightforward and work seamlessly with AD.
---
### Author:  [@amrods](https://discourse.julialang.org/u/amrods)
#### Post date: [April 30, 2021, 1:21pm UTC](https://discourse.julialang.org/t/status-of-integration-optimization-and-nl-solvers/60025/17 "2021-04-30T13:21:01Z")
Hmm, this is a very good suggestion. I had defined a function
```julia
function discountedintegral(f, r, T; nodes=5_000)
x, w = gausslaguerre(nodes)
exp(-r*T)/r * sum( w[i]*f(x[i]/r + T) for i in 1:nodes )
end
```
but that computes the nodes every time it is called. With your suggestion, it would be like this:
```julia
struct Hquad{T, V <: AbstractVector{T}}
nodes::V
weights::V
r::T
end
function (hq::Hquad)(f, T)
x = hq.nodes
w = hq.weights
r = hq.r
exp(-r*T)/r * sum( w[i]*f(x[i]/r + T) for i in 1:length(x) )
end
hq = Hquad(gausslaguerre(1_000)..., 0.04)
```
I’d wish there was a way to prevent hardcoding the number of nodes. Perhaps with memoization?
---
### Author:  [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [April 30, 2021, 1:25pm UTC](https://discourse.julialang.org/t/status-of-integration-optimization-and-nl-solvers/60025/18 "2021-04-30T13:25:32Z")
I would just save those nodes and weights in a constant once.
And unless your H is really nasty, I think you should be fine with much less than 1000. I would start with 20; when it works, Gaussian quadrature is amazingly accurate for most purposes.
---
### Author:  [@amrods](https://discourse.julialang.org/u/amrods)
#### Post date: [April 30, 2021, 1:36pm UTC](https://discourse.julialang.org/t/status-of-integration-optimization-and-nl-solvers/60025/19 "2021-04-30T13:36:59Z")
> [@Tamas\_Papp](#):
>
> I would just save those nodes and weights in a constant once.
I thought the point of defining a `struct` was to carry those nodes and weights within the struct itself. Perhaps you mean to store the number of nodes in a `const`? If not, I could have saved the nodes and weights and just employ a standard function that uses those consts.
---
### Author:  [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [April 30, 2021, 1:40pm UTC](https://discourse.julialang.org/t/status-of-integration-optimization-and-nl-solvers/60025/20 "2021-04-30T13:40:57Z")
I would save a version transformed by \rho in the struct already (especially if it does not change often), but the way you propose is fine too.
[Next page](https://discourse.julialang.org/t/status-of-integration-optimization-and-nl-solvers/60025.md?page=2)