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.
and the code for this is
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.