I’m trying to solve a constrained optimization problem with JuMP and HiGHs (this is for an open source package so I want to rely on an open source solver), but for some reason the optimization never terminates although the objective value reaches a minimum very quickly.
I’m posting an MWE (thanks to the invaluable help of @jd-foster on Slack) with real world data below. The problem to solve is:
where Z is a length-k vector, and X is a (k \times J) matrix. The idea is that Z is an observation with k characteristics, and X holds J other observations with measurements on the same characteristics. The objective is to find a weighted combination of observations (columns) in X to approximate Z, with the first term in the objective penalising the overall distance between Z and the weighted average and the second term the pair-wise distance between Z and the columns of X included in the weighted average.
(For those interested this is the “Penalized Synthetic Control” estimator from Abadie and L’Hour (2021))
The runnable code is:
using Downloads, DelimitedFiles, LinearAlgebra import JuMP: JuMP, Model, @variable, @constraint, @objective import HiGHS Z = [5.17119501, 5.4708669, 2.5463325, 2.7430484, 0.0, 0.0, 0.0, 2.19362429, 2.19362429, 2.03571710] X = readdlm(Downloads.download("https://gist.githubusercontent.com/nilshg/a82b93cc08b6135fabe65400c536f6f0/raw/4b24d296e40b2efd5765db74e270a1a8b0bb5997/X.dlm")) function create_model(Z, X, λ; optimizer=HiGHS.Optimizer) model = Model(optimizer) J = size(X, 2) @variable(model, 0 <= w[1:J] <= 1) @constraint(model, sum(w[i] for i ∈ 1:J) == 1) @objective(model, Min, sum( (Z .- X*w).^2 ) + λ * sum( w[j]* sum( (Z .- X[:,j])).^2 for j in 1:J) ) return model end m = create_model(Z, X, 0.1) JuMP.set_attribute(m, "time_limit", 120.0) JuMP.optimize!(m)
which for me gives:
Iteration, Runtime, ObjVal, NullspaceDim 0, 0.356691, 24.143424, 0 1000, 86.010277, 0.587094, 7 1386, 120.059697, 0.587094, 6 Model status : Time limit reached QP ASM iterations: 1386 Objective value : 5.8709427083e-01 HiGHS run time : 120.09
Playing around with the
time_limit I found that the final objective value is actually reached after ~20 iterations or around 2 seconds, and the solver keeps running forever afterwards (I’ve let it run for about an hour before terminating).
It seems to me that there should be a way to get HiGHs to return here, but playing around with the various
tolerance related parameters doesn’t seem to do anything.
Lastly, @jd-foster solved this problem in Gurobi which suggested perturbing the diagonal term in the quadratic objective. This is not something I’d like to do in my package as default (as the diagonal is meaningful and users are likely to feed in Z and X in standardised in such a way that they expect an identity diagonal) but it might help diagnose the issue that HiGHs is facing here.
In summary, my questions are:
- How can I get HiGHs to return a value when the objective stops changing; and
- Are there any strategies to make my solution strategy more robust to whatever is happening here in the context of a package which should be able to deal with all sorts of user provided data?