# Optim interface.jl type instability and efficient exception handling

**URL:** https://discourse.julialang.org/t/optim-interface-jl-type-instability-and-efficient-exception-handling/84972
**Category:** Performance
**Tags:** optim, error, exception, optimization
**Created:** [July 29, 2022, 11:32am UTC](https://discourse.julialang.org/t/optim-interface-jl-type-instability-and-efficient-exception-handling/84972 "2022-07-29T11:32:28Z")
**Posts on this page:** 1
**Page:** 1

<div class="post-metadata">

### Author: ![psterzinger](https://avatars.discourse-cdn.com/v4/letter/p/df705f/32.png) [@psterzinger](https://discourse.julialang.org/u/psterzinger)
#### Post date: [July 29, 2022, 11:32am UTC](https://discourse.julialang.org/t/optim-interface-jl-type-instability-and-efficient-exception-handling/84972/1 "2022-07-29T11:32:29Z")

</div>

Hi,

I am currently trying to optimize a simulation that involves repeated parameter estimation, which I do with Optim.jl’s `optimize` with supplied gradients and hessians.

**interjace.jl type instability**

Looking at the flamegraph of one such simulation run, it seems that there are type instabilities in the internals of optimization. In particular, I get type warnings for:

- interface.jl #optimize#89: line 82
- interface.jl value\_gradient!!: line 82
- interface.jl hessian!!: line 95

The flame graph for one such optimization run looks like this, where the intermitted red bars correspond to the type instabilities from the interface.jl calls.

 ![Bildschirmfoto 2022-07-29 um 12.08.17](https://global.discourse-cdn.com/julialang/original/3X/9/9/998a8c93c5e1ddf78759b0ac93ef6c4d73568ebc.png)

and the code for this is

```julia
using LinearAlgebra, Optim, ProfileView, Profile, Random, Distributions

n = 200; 
p = 11; 
beta_0 = randn(11); 
beta = randn(p); 
X = randn(n,p); 
X[:,1] .= 1.0; 
y = convert.(Float64, rand.(Binomial.(1, cdf.(Logistic(), X * beta_0)))); 
beta_init = fill(0.0,11); 
buff_X = similar(X); 
mus = randn(n); 
ws = similar(mus); 

function get_logistic_MLE2!(beta::Vector{Float64},X::Matrix{Float64},y::Vector{Float64},beta_init::Vector{Float64},mus::Vector{Float64},ws::Vector{Float64},buff_X::Matrix{Float64}; ϵ::Float64 = 0.0 , kwargs...) 
    opt = Optim.optimize(β -> logistic_loglikl(β,X,y,mus,ws;ϵ),
                        (g,β)->logistic_score!(g,β,X,y,mus;ϵ),
                        (H,β)->logistic_hessian!(H,β,X,mus,ws,buff_X;ϵ), 
                        beta_init;kwargs...)
    beta .= opt.minimizer
    return Optim.converged(opt) && !any(isnan.(beta)) , Optim.minimum(opt) 
end 

function logistic_loglikl(beta::Vector{Float64},X::Matrix{Float64},y::Vector{Float64},mus::Vector{Float64},buff_mus::Vector{Float64}; ϵ = 0.0)
    mul!(mus,X,beta)
    buff_mus .= y .* mus
    mus .= 1.0 .+ exp.(mus) 
    mus .= log.(mus)
    if !(ϵ == 0.0)
        buff_mus .= (1.0 + ϵ/2) .* buff_mus
        mus .= (1.0 + ϵ) .* mus 
    end 
    mus .= mus .- buff_mus 
    return sum(mus)
 end 

 function logistic_score!(g::Vector{Float64},beta::Vector{Float64},X::Matrix{Float64},y::Vector{Float64},mus::Vector{Float64}; ϵ = 0.0) 
    logistic_mean!(mus,X,beta) 
    if ϵ == 0.0
        mus .= mus .- y 
    else 
        mus .= (1.0 + ϵ) .* mus .- (1.0 + ϵ/2) .* y
    end 
    mul!(g, X', mus) 
    return g
end 

function logistic_hessian!(H::Matrix{Float64},beta::Vector{Float64},X::Matrix{Float64},mus::Vector{Float64},ws::Vector{Float64},buff_X::Matrix{Float64}; ϵ = 0.0) 
    logistic_mean!(mus,X,beta)
    ws .= mus .* (1.0 .- mus) 
    colmult!(buff_X,ws,X)
    mul!(H,X',buff_X)
    if !(ϵ == 0.0) 
        H .= (1.0 + ϵ) .* H 
    end 
    return H 
end

function logistic_score!(g::Vector{Float64},beta::Vector{Float64},X::Matrix{Float64},y::Vector{Float64},mus::Vector{Float64}; ϵ = 0.0) 
    logistic_mean!(mus,X,beta) 
    if ϵ == 0.0
        mus .= mus .- y 
    else 
        mus .= (1.0 + ϵ) .* mus .- (1.0 + ϵ/2) .* y
    end 
    mul!(g, X', mus) 
    return g
end 

function logistic_mean!(mus::Vector{Float64},X::Matrix{Float64},beta::Vector{Float64})  
    mul!(mus,X,beta,-1.0,false) 
    mus .= exp.(mus) 
    mus .= 1.0 ./ (1.0 .+ mus)
end 

function colmult!(A::Matrix{Float64},v::Vector{Float64}, B::Matrix{Float64}) 
    @inbounds for j = 1:size(A,2) 
        @inbounds for i = 1:size(A,1) 
            A[i,j] = v[i] * B[i,j]
        end
    end 
    A
end

Profile.init(10000000, 0.00001)
ProfileView.@profview get_logistic_MLE2!(beta,X,y,beta_init,mus,ws,buff_X)

```

Is there anything wrong with how I set up my functions that causes this type instability and/or is there anything I can do to get rid of it?

**Exception handling**

My simulation essentially involves repeatedly drawing some data, estimating some parameters and do some calculations, which involves solving some nonlinear systems of equations whcih I do using NLsolve.jl’s `nlsolve`. With that, there are some data configurations for which estimation will fail, and I will get some error, usually some `LAPACKException, SingularException, NLsolve.IsFiniteException,DomainError`…  
If estimation fails, all I want to do is set my parameter estimates and other things to NaN, and continue, I do not care about the type of exception.

At the moment I am wrapping all things that can fail in a try-catch block, which works, but appears to be quite costly. Is there a clever/efficient way to simply return an error without throwing it, so that I can get rid of my try-catch statements, e.g. by telling julia to expect either, say a `Vector{Float64}` or some exception as return type.

Thank you in advance for the help.