JuMP Ipopt performance for MLE

Specific Domains Optimization (Mathematical)
,
1

Hello,

With much of your help, I am finally able to write a functioning code file. However, I was not able to find enough resources online for tuning the performance of Ipopt in JuMP.
Now I would like to seek your advice on how to achieve the following goals (or to find resources for achieving these goals):

  1. optimize my code for better performance
  2. incorporate codes that speed up the optimization process
  3. obtain multiple solutions if possible
  4. export optimal solution(s) to .csv file.

Below is my MWE (@odow I gave it a try!!):


using CSV, DataFrames, StatsBase
using JuMP, Ipopt, Random

# read data
const df = DataFrame(owner = ["b","n","s"], 
                    b_ln_prod = [0.5,0.25,0.1],
                    s_ln_prod=[0.1,0.25,0.5],
                    b_delta =[1,2,3],
                    s_delta=[3,2,1]);
const delta_idx = [1;2;3];

# data dimensions
const n = nrow(df);
const n_delta = length(delta_idx);
const n_parameter = n_delta + 5;

function mle()
    # set model parameters
    Random.seed!(1234)
    model = Model(Ipopt.Optimizer)
    #set_optimizer_attribute(model, "max_cpu_time", 60.0)
    set_optimizer_attribute(model, "print_level", 0)
    set_optimizer_attribute(model, "max_iter", 100)
    set_optimizer_attribute(model, "tol", 1e-4)
    set_time_limit_sec(model, 60.0)

    # add variables
    @variables(model, begin
        0<=ρ<=1, (start=0.5, base_name = "rho")
        0<=α<=1, (start=0.6, base_name = "alpha")
        0<=fB<=10, (start = 0.0, base_name = "fB")
        0<=fS<=10, (start = 0.0, base_name = "fS")
        0<=fO<=10, (start = 0.0, base_name = "fO")
        0<=δ[i=1:n_delta]<=1, (start=0.5, base_name = "delta$i")
    end)

    # constraint
    @constraint(model, con, α - ρ >= 0.001)

    # objective
    function likelihood(ρ::T, α::T, fB::T, fS::T, fO::T, δ::T...) where {T}
        sum = zero(T);
        for i = 1:n
            k = df.owner[i];
            θb = df.b_ln_prod[i];
            θs = df.s_ln_prod[i];
            b_idx = df.b_delta[i];
            s_idx = df.s_delta[i];
            δb = δ[b_idx];
            δs = δ[s_idx];
            βB = 0.5 + 0.5*δb^α;
            βS = 0.5 - 0.5*δs^α;
            βO = 0.5;
            ψB = α^(α/(1-α))*((1-α*βB)*(βB*θb)^(ρ/(1-ρ))+(1-α*(1-βB))*((1-βB)*θs)^(ρ/(1-ρ)))/((βB*θb)^(ρ/(1-ρ))+((1-βB)*θs)^(ρ/(1-ρ)))^((ρ-α)/(ρ*(1-α)));
            ψS = α^(α/(1-α))*((1-α*βS)*(βS*θb)^(ρ/(1-ρ))+(1-α*(1-βS))*((1-βS)*θs)^(ρ/(1-ρ)))/((βS*θb)^(ρ/(1-ρ))+((1-βS)*θs)^(ρ/(1-ρ)))^((ρ-α)/(ρ*(1-α)));
            ψO = α^(α/(1-α))*((1-α*βO)*(βO*θb)^(ρ/(1-ρ))+(1-α*(1-βO))*((1-βO)*θs)^(ρ/(1-ρ)))/((βO*θb)^(ρ/(1-ρ))+((1-βO)*θs)^(ρ/(1-ρ)))^((ρ-α)/(ρ*(1-α)));
            temp = ((k=="b" ? 1 : 0)*(exp(ψB-fB))/(exp(ψB-fB)+exp(ψS-fS)+exp(ψO-fO))
                    +(k=="s" ? 1 : 0)*(exp(ψS-fS))/(exp(ψB-fB)+exp(ψS-fS)+exp(ψO-fO))
                    +(k=="n" ? 1 : 0)*(exp(ψO-fO))/(exp(ψB-fB)+exp(ψS-fS)+exp(ψO-fO)))
            sum = sum + log(temp)
        end
        return -sum
    end
    register(model, :likelihood, n_parameter, likelihood; autodiff = true)
    @NLobjective(model, Min, likelihood(ρ,α,fB,fS,fO,δ...))

    # optimization
    JuMP.optimize!(model)

    # display results
    if termination_status(model) == MOI.OPTIMAL
        println("Solution is optimal")
    elseif termination_status(model) == MOI.TIME_LIMIT && has_values(model)
        println("Solution is suboptimal due to a time limit, but a primal solution is available")
    else
        error("The model was not solved correctly.")
    end
    println("  objective value = ", objective_value(model))
end
2

Your code looks okay.

You may get better performance by extracting common subexpressions into @NLexpression. It looks like there are a lot of similar terms there!

One potential issue with this is that we disable second-order derivative information for user-defined functions. Instead of writing out the summation in a function, you could build a nonlinear expression for each I, so something like this (untested. Probably has typo’s etc. But it should point you in the right direction)

expr = Any[]
for i = 1:n
    k = df.owner[i];
    θb = df.b_ln_prod[i];
    θs = df.s_ln_prod[i];
    b_idx = df.b_delta[i];
    s_idx = df.s_delta[i];
    δb = δ[b_idx];
    δs = δ[s_idx];
    βB = @NLexpression(model, 0.5 + 0.5*δb^α);
    βS = @NLexpression(model, 0.5 - 0.5*δs^α);
    βO = 0.5;
    ψB = @NLexpression(model, α^(α/(1-α))*((1-α*βB)*(βB*θb)^(ρ/(1-ρ))+(1-α*(1-βB))*((1-βB)*θs)^(ρ/(1-ρ)))/((βB*θb)^(ρ/(1-ρ))+((1-βB)*θs)^(ρ/(1-ρ)))^((ρ-α)/(ρ*(1-α))))
    ψS = @NLexpression(model, α^(α/(1-α))*((1-α*βS)*(βS*θb)^(ρ/(1-ρ))+(1-α*(1-βS))*((1-βS)*θs)^(ρ/(1-ρ)))/((βS*θb)^(ρ/(1-ρ))+((1-βS)*θs)^(ρ/(1-ρ)))^((ρ-α)/(ρ*(1-α))))
    ψO = @NLexpression(model, α^(α/(1-α))*((1-α*βO)*(βO*θb)^(ρ/(1-ρ))+(1-α*(1-βO))*((1-βO)*θs)^(ρ/(1-ρ)))/((βO*θb)^(ρ/(1-ρ))+((1-βO)*θs)^(ρ/(1-ρ)))^((ρ-α)/(ρ*(1-α))))
    if k == "b"
        push!(model, @NLexpression(model, (exp(ψB-fB))/(exp(ψB-fB)+exp(ψS-fS)+exp(ψO-fO))))
    elseif k == "s"
        # ...
    end
end
@NLobjective(model, Min, -sum(log(expr[i] for i in 1:n)))

To get better performance from Ipopt, use a different linear solver: https://github.com/jump-dev/Ipopt.jl#linear-solvers

Multiple solutions are not possible.

Use value(ρ) to obtain the solution. After that it’s just a normal Julia variable, so you can do anything you like with it. See, e.g., dictionary - Compact way to save JuMP optimization results in DataFrames - Stack Overflow

2 Likes