Different results when registering the log-likelihood function vs pasting its contents

I’m trying to optimize this log-likelihood function:

-log(σ^2) / 2 - ϵ^2 / 2 + log( 2Φ(a * ϵ) )

Here ϵ depends on data and some more parameters.

I’m trying to optimize this in two ways:

  1. Put the sum of these likelihoods (over the data) right into @NLobjective
  2. Put the likelihood into a function, register it and put the sum over invocations of this function into @NLobjective.

I expected the results (parameter estimates) to be the same, but they’re different, even though the likelihood is exactly the same in both cases.

The code looks like this:

    import Pkg

        "CSV", "DataFrames", "Ipopt", "JuMP", "SpecialFunctions"
    ], io=devnull)

using CSV, DataFrames

using JuMP
import Ipopt

using SpecialFunctions: erf

const AV = AbstractVector{T} where T;

"Normal CDF"
Φ(x::Real) = (1 + erf(x / sqrt(2))) / 2

# 1. Load data
df = CSV.read("data_sample.csv", DataFrame)

# 2. Define log-likelihood for one observation
logLik_SF(ϵ::Real, a::Real, σ::Real) = -log(σ^2) / 2 - ϵ^2 / 2 + log( 2Φ(a * ϵ) )

# 3. Define function for estimation
function estimate_SF(Y::AV{<:Real}, L::AV{<:Real}, K::AV{<:Real}; optimizer=Ipopt.Optimizer, ll_sum::Bool=true)
    N = length(Y)
    y = log.(Y)
    l = log.(L)
    k = log.(K)
    model = Model(optimizer)
    set_optimizer_attribute(model, MOI.Silent(), true)
    @variable(model, α, start=0.5) # elasticity
    @variable(model, a, start=0) # skewed normal shape
    @variable(model, σ ≥ 0) # skewed normal scale
    ϵ = [
            (y[i] - (α * l[i] + (1 - α) * k[i])) / σ
        for i ∈ eachindex(y)
    if ll_sum
        register(model, :Φ, 1, Φ, autodiff=true)

            model, Max,
            # Pasted the content of `logLik_SF`
            # and replaced `ϵ` with `ϵ[i]`
                -log(σ^2) / 2 - ϵ[i]^2 / 2 + log( 2Φ(a * ϵ[i]) )
                for i ∈ eachindex(ϵ)
        register(model, :logLik_SF, 3, logLik_SF, autodiff=true)
            model, Max,
            # Just call `logLik_SF`
                logLik_SF(ϵ[i], a, σ)
                for i ∈ eachindex(ϵ)

println("Estimate using hard-coded expression:")
model = estimate_SF(
    df[!, :Y], df[!, :L], df[!, :C],

a, α, σ = value(model[:a]), value(model[:α]), value(model[:σ])
@show a α σ

println("\nEstimate using function call:")
model = estimate_SF(
    df[!, :Y], df[!, :L], df[!, :C],

a, α, σ = value(model[:a]), value(model[:α]), value(model[:σ])
@show a α σ

The data come from data_sample.csv: Some data - Pastebin.com.

I’m running the code like this: julia-1.7 report.jl


      Status `/private/var/folders/ys/3h0gnqns4b98zb66_vl_m35m0000gn/T/jl_vWtB7k/Project.toml`
  [336ed68f] CSV v0.9.11
  [a93c6f00] DataFrames v1.3.0
  [b6b21f68] Ipopt v0.7.0
  [4076af6c] JuMP v0.21.5
  [276daf66] SpecialFunctions v2.0.0
Estimate using hard-coded expression:
a = -0.5226232660699477
α = 0.1438052406305529
σ = 0.2727175471168534

Estimate using function call:
a = 1.5968811212735676
α = 0.18443033301359052
σ = 0.3456179651592732

As you can see, I’m getting different results depending on whether I use -log(σ^2) / 2 - ϵ[i]^2 / 2 + log( 2Φ(a * ϵ[i]) ) or logLik_SF(ϵ[i], a, σ) in the sum in the objective function. …but they’re literally the same exact function!

What could be the problem here?

Random observation (that doesn’t solve the issue):

  • installing packages like in the post installs:
    • “old” Ipopt v0.7.0
    • “old” JuMP v0.21.5
    • new SpecialFunctions v2.0.0
  • installing Ipopt and JuMP first, and then SecialFunctions installs:
    • much newer Ipopt v0.9.1
    • new JuMP v0.22.1
    • old SpecialFunctions v1.8.1

Multiple equilibria in package management strike again! The problem happens on all of these versions, though…

If you look at the Ipopt log of the two problems, they are not equivalent. In the register case, we disable second-order derivatives: Feature Request: Second Derivatives for User Defined Functions · Issue #1198 · jump-dev/JuMP.jl · GitHub

Multiple equilibria in package management strike again!

Again, this is expected. The package manager does’t know whether you want the latest JuMP or the latest SpecialFunctions. Use ] add JuMP@0.22 to force it to choose JuMP.

Yeah, I get this now - I’m just still kinda frustrated with this, but that’s beside the point.

Okay, got it, thanks!