# Suggestions needed: speed up optimization

**URL:** https://discourse.julialang.org/t/suggestions-needed-speed-up-optimization/47449
**Category:** Performance
**Created:** [September 28, 2020, 11:25pm UTC](https://discourse.julialang.org/t/suggestions-needed-speed-up-optimization/47449 "2020-09-28T23:25:04Z")
**Posts on this page:** 12
**Page:** 2

<div class="post-metadata">

### Author: ![DShiu](https://avatars.discourse-cdn.com/v4/letter/d/8e8cbc/32.png) [@DShiu](https://discourse.julialang.org/u/DShiu)
#### Post date: [October 2, 2020, 10:08am UTC](https://discourse.julialang.org/t/suggestions-needed-speed-up-optimization/47449/21 "2020-10-02T10:08:14Z")

</div>

Agreed. After reading documentation and discussions I realized this issue, but haven’t figured out how to get it work ☹

---

<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: [October 2, 2020, 12:12pm UTC](https://discourse.julialang.org/t/suggestions-needed-speed-up-optimization/47449/22 "2020-10-02T12:12:36Z")

</div>

> [@DShiu](#):
>
> I revised the function be based on comprehension

This will allocate an array on the heap each time it is called, which is bad for performance.

---

<div class="post-metadata">

### Author: ![mikkoku](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mikkoku/32/16274_2.png) [@mikkoku](https://discourse.julialang.org/u/mikkoku)
#### Post date: [October 3, 2020, 6:13am UTC](https://discourse.julialang.org/t/suggestions-needed-speed-up-optimization/47449/23 "2020-10-03T06:13:29Z")

</div>

Could you just update negloglike inside the loop? That way you don’t need to worry about the prob vector at all.

Have you considered a non-adaptive numerical integration? In this case I see two possible problems with adaptive quadrature. 1) There is a sum and log transforms after it and thus effect of rtol in the end is not clear to me. 2) Since you are using somewhat low precision (rtol=0.001) it might cause the choice\_prob\_naive to have spurious jumps if the adaptation changes the quadrature points in a discontinuous manner. This could be bad for optimization.

---

<div class="post-metadata">

### Author: ![DShiu](https://avatars.discourse-cdn.com/v4/letter/d/8e8cbc/32.png) [@DShiu](https://discourse.julialang.org/u/DShiu)
#### Post date: [October 3, 2020, 5:00pm UTC](https://discourse.julialang.org/t/suggestions-needed-speed-up-optimization/47449/24 "2020-10-03T17:00:48Z")

</div>

What you do mean by “update negloglike inside the loop”? Compute it element by element in a loop? Doesn’t that require pre-allocating and gives me the trouble that Tamas\_Papp pointed out?

I replaced rtol with with f\_abstol. Would this help with the issue?

---

<div class="post-metadata">

### Author: ![mikkoku](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mikkoku/32/16274_2.png) [@mikkoku](https://discourse.julialang.org/u/mikkoku)
#### Post date: [October 3, 2020, 7:01pm UTC](https://discourse.julialang.org/t/suggestions-needed-speed-up-optimization/47449/25 "2020-10-03T19:01:43Z")

</div>

I mean something like this

```julia
function loglike(θ, x, y, R, c, d, p0, B)
  negloglike = 0.0
  for i = 1:size(B,1)
    prob = choice_prob_naive(θ, r[i], x[i,:], y[i,:], R[i], c[i], d[i], p0[i])
    negloglike += - (B[i]*log(prob) + (1 - B[i])*log(1 - prob))
  end
  return negloglike
end

```

I think the same issues apply to any kind of adaptive quadrature, if it is necessary to use a largeish tolerance.

I recommend checking  
[https://github.com/JuliaMath/QuadGK.jl#gaussian-quadrature-and-arbitrary-weight-functions](https://github.com/JuliaMath/QuadGK.jl#gaussian-quadrature-and-arbitrary-weight-functions)  
for non-adaptive gaussian quadrature

---

<div class="post-metadata">

### Author: ![DShiu](https://avatars.discourse-cdn.com/v4/letter/d/8e8cbc/32.png) [@DShiu](https://discourse.julialang.org/u/DShiu)
#### Post date: [October 3, 2020, 8:21pm UTC](https://discourse.julialang.org/t/suggestions-needed-speed-up-optimization/47449/26 "2020-10-03T20:21:53Z")

</div>

Thanks, this is helpful! I didn’t realize that `negloglike = 0.0` would not trigger the same issue for autodiff as the one triggered by pre-allocating using `prob = zeros(size(B,1))`.

My understanding of the tolerance problem from reading is that the adaptive integtion subdivides the space to do approximate integration until it is considered “good enough” and the number of subdivision depends on my tolerance level. If I switch to the gaussian quadrature, I impose a fixed number of subdivision and the number can be based on a tiny rtol. Then I apply this number to each integration in the loop. It this right?

---

<div class="post-metadata">

### Author: ![bashonubuntu](https://avatars.discourse-cdn.com/v4/letter/b/f19dbf/32.png) [@bashonubuntu](https://discourse.julialang.org/u/bashonubuntu)
#### Post date: [October 3, 2020, 8:50pm UTC](https://discourse.julialang.org/t/suggestions-needed-speed-up-optimization/47449/27 "2020-10-03T20:50:42Z")

</div>

One more general point is that you can pass arrays of parametric types if you wanted to i.e. you can define a function like

```julia
function f(x::Array{T}) where {T}
        return Array{T}(repeat([0], outer = [10, 1]))
       end

julia> f(fill(0.0, 2))
10×1 Array{Float64,2}:
 0.0
 0.0
 0.0
 ⋮
 0.0
 0.0
 0.0

julia> f(fill(0, 2))
10×1 Array{Int64,2}:
 0
 0
 0
 ⋮
 0
 0
 0

julia> f(Array{Real}(undef, 0, 0))
10×1 Array{Real,2}:
 0
 0
 0
 ⋮
 0
 0
 0

```

So, in your code, something like

```julia
function loglike(θ::Array{T}, r, x y, R, c, d, p0, B) where {T}
    #prob = zeros(size(B,1))
    prob = Array{T}(repeat([0], outer = [size(B)[1], 1]))
    for i = 1:size(B,1)    
        prob[i] = choice_prob_naive(θ, r[i], x[i,:], y[i,:], R[i], c[i], d[i], p0[i])
    end
    negloglike = - (B'*log.(prob) + (1.0 .- B)'*log.(1.0 .- prob))  
    return negloglike
end

```

should be fine (but also check how you have defined `choice_prob_naive()`.

---

<div class="post-metadata">

### Author: ![mikkoku](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mikkoku/32/16274_2.png) [@mikkoku](https://discourse.julialang.org/u/mikkoku)
#### Post date: [October 4, 2020, 5:39am UTC](https://discourse.julialang.org/t/suggestions-needed-speed-up-optimization/47449/28 "2020-10-04T05:39:34Z")

</div>

Computing Dual + Float64 results in Dual, which means that type of negloglike changes. In general changing types is bad, but it might not be a problem here, since the computation happens in choice\_prob\_naive. It would be better to use `zero(T)` where T is the eltype of theta.

Yes, more specifically the adaptation chooses only some intervals, where the estimated approximation error is largest and subdivides them. I think you can also pre-compute the integration points and weights and use those for each integration.

In the end you probably should experiment with adaptive and non-adaptive quadrature and autodiff and derivative free optimization, to see which combination works best.

---

<div class="post-metadata">

### Author: ![DShiu](https://avatars.discourse-cdn.com/v4/letter/d/8e8cbc/32.png) [@DShiu](https://discourse.julialang.org/u/DShiu)
#### Post date: [October 4, 2020, 6:44am UTC](https://discourse.julialang.org/t/suggestions-needed-speed-up-optimization/47449/29 "2020-10-04T06:44:54Z")

</div>

I am still confused about the type issue related to autodiff. With type annotations removed, the optimization runs fine but when I try to get the Hessian from the results, Julia complains about the type issue…

```julia
td = TwiceDifferentiable(θ -> loglike(θ, r, x, y, R, c, d, p0, B), theta; autodiff = :forward);
@time theta_naive = optimize(td, theta, BFGS(), Optim.Options(f_abstol=10^-3, iterations=100_000, show_trace=true, show_every=1))

julia> H = Optim.hessian!(td, theta_naive)
ERROR: MethodError: no method matching iterate(::Optim.MultivariateOptimizationResults{BFGS{LineSearches.InitialStatic{Float64},LineSearches.HagerZhang{Float64,Base.RefValue{Bool}},Nothing,Nothing,Flat},Float64,Array{Float64,1},Float64,Float64,Array{OptimizationState{Float64,BFGS{LineSearches.InitialStatic{Float64},LineSearches.HagerZhang{Float64,Base.RefValue{Bool}},Nothing,Nothing,Flat}},1},Bool})

```

---

<div class="post-metadata">

### Author: ![mikkoku](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mikkoku/32/16274_2.png) [@mikkoku](https://discourse.julialang.org/u/mikkoku)
#### Post date: [October 4, 2020, 8:01am UTC](https://discourse.julialang.org/t/suggestions-needed-speed-up-optimization/47449/30 "2020-10-04T08:01:55Z")

</div>

The error is implying that hessian! doesn’t expect an OptimizationResult. You are missing Optim.minimizer.

---

<div class="post-metadata">

### Author: ![DShiu](https://avatars.discourse-cdn.com/v4/letter/d/8e8cbc/32.png) [@DShiu](https://discourse.julialang.org/u/DShiu)
#### Post date: [October 4, 2020, 4:08pm UTC](https://discourse.julialang.org/t/suggestions-needed-speed-up-optimization/47449/31 "2020-10-04T16:08:13Z")

</div>

Hmmm I see. How do you get this information from reading the error message? I often find these messages hard to understand.  
Should we expect `Optim.hessian!` to run for a long time, say, 10% of the optimization?

---

<div class="post-metadata">

### Author: ![tomerarnon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tomerarnon/32/3170_2.png) [@tomerarnon](https://discourse.julialang.org/u/tomerarnon)
#### Post date: [October 4, 2020, 7:06pm UTC](https://discourse.julialang.org/t/suggestions-needed-speed-up-optimization/47449/32 "2020-10-04T19:06:19Z")

</div>

> [@DShiu](#):
>
> How do you get this information from reading the error message?

```julia
ERROR: MethodError: no method matching iterate(::Optim.MultivariateOptimizationResults{...})

```

Says that there is “no method matching” i.e. “no known way of doing:” `iterate(blah)` where blah is the type of thing it got as an input. But you didn’t call `iterate`, so that may seem strange. The stacktrace will show that `iterate` is being called from something (that is being called from something that is being called from… etc.) that is eventually called from `hessian!`. Generally with stacktraces, you look for the last line in the stack that you “own” to figure out what you did wrong.

[Previous page](https://discourse.julialang.org/t/suggestions-needed-speed-up-optimization/47449.md?page=1)
