Getting HiGHs to return when objective value doesn't move

nilshg · September 18, 2023, 2:56pm

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:

\begin{align} \min_{w} \left| \left| Z - \sum_{j} w_jX_j \right| \right|^2 + \lambda \sum_{j}w_j \left| \left| Z - X_j \right| \right|^2 \\ s.t. \ 0 \leq w_j \leq 1 \ \forall \ j \\ \sum_{j} w_j = 1 \end{align}

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 gap or 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?

blob · September 19, 2023, 6:59am

I changed the solver to OSQP (also open source, but not for mixed integer problems - in case that matters) and I got a solution in a few seconds:

The value of the objective is different in OSQP and in HiGHS, though. When I run the code with HiGHS, I get similar results to you. However, when I check the solution summary, I see that the duals are infeasible:

So possibly a reason for not finishing is this infeasibility. I don’t know if the infeasibility is a bug in HiGHS or a property of the problem. Just for the sake of completeness, Ipopt was able to find the same solution as HiGHS:

with similar dual values.

And SCS failed complaining that the matrix is not positive definite:

I am not sure a linear solver or a solver requiring positive semi-definite matrix Q can be robust to all sorts of user provided data. Maybe the way to go is to check whether the matrices created from data satisfy the assumptions of a solver and warn the user if this is not the case?

odow · September 19, 2023, 9:12pm

There are a couple of things going on here.

The matrix has a large condition number:

λ = 0.1
model = Model()
J = size(X, 2)
@variable(model, 0 <= w[1:J] <= 1);
expr = JuMP.@expression(model,
    sum( (Z .- X*w).^2 ) + λ * sum( w[j]* sum( (Z .- X[:,j])).^2 for j in 1:J)
);
vmap = Dict(w[j] => j for j in 1:J);
Q = zeros(J, J);
for (k, v) in expr.terms
    Q[vmap[k.a], vmap[k.b]] += v
    Q[vmap[k.b], vmap[k.a]] += v
end

julia> LinearAlgebra.cond(Q)
4.079491602024465e22

The QP solver in HiGHS is immature and still seeing development/bug fixes.

I’d use Gurobi or Ipopt, and I’d write your model like this:

function create_model2(Z, X, λ; optimizer=Ipopt.Optimizer)
    model = Model(optimizer)
    J = size(X, 2)
    @variable(model, 0 <= w[1:J] <= 1)
    @constraint(model, sum(w) == 1)
    @variable(model, residuals[1:length(Z)])
    @constraint(model, residuals .== Z .- X*w)
    coefficients = [sum(Z .- X[:,j]).^2 for j in 1:J]
    @objective(model, Min, sum(residuals.^2) + λ * coefficients' * w)
    return model
end

This moves the data into linear equality constraints (easy) and it simplifies the quadratic to a diagonal matrix.

jd-foster · September 20, 2023, 3:15am

This works very well from my testing; it has the equivalent objective and solution as the other versions, and is a few orders of magnitude quicker as bonus.

MWE for anyone interested:

import DelimitedFiles
import JuMP: JuMP, Model, @variable, @constraint, @objective
import Ipopt

Z = [5.1711950141609835, 5.470866932205469, 2.546332563474285, 2.743048443920391, 0.0, 0.0, 0.0, 2.1936242978411693, 2.1936242978411697, 2.0357171080835728]
Base.download("https://gist.githubusercontent.com/nilshg/a82b93cc08b6135fabe65400c536f6f0/raw/4b24d296e40b2efd5765db74e270a1a8b0bb5997/X.dlm", "X.dlm")
X = DelimitedFiles.readdlm("X.dlm")

function create_model2(Z, X, λ; optimizer=Ipopt.Optimizer)
    model = Model(optimizer)
    J = size(X, 2)
    @variable(model, 0 <= w[1:J] <= 1)
    @constraint(model, sum(w) == 1)
    @variable(model, residuals[1:length(Z)])
    @constraint(model, residuals .== Z .- X*w)
    coefficients = [sum(Z .- X[:,j]).^2 for j in 1:J]
    @objective(model, Min, sum(residuals.^2) + λ * coefficients' * w)
    return model
end

function run(Z, X, λ; optimizer=Ipopt.Optimizer, filter_tol = 1e-4)
    model2 = create_model2(vec(Z), X, 0.1; optimizer=optimizer)
    JuMP.optimize!(model2)

    obj_val = JuMP.objective_value(model2)
    ## Non-trivial values:
    sol_nontrivial = [(i, val) for (i, val) in enumerate(JuMP.value.(model2[:w])) if val > filter_tol]

    return obj_val, sol_nontrivial
end

@time obj_val, sol = run(Z, X, 0.1);
## 0.294944 seconds (50.57 k allocations: 9.058 MiB)

with

julia> obj_val
0.5853730907635923

julia> sol
7-element Vector{Tuple{Int64, Float64}}:
 (13, 0.18656265907917816)
 (44, 0.032810328329308656)
 (1657, 0.007132053670042325)
 (1671, 0.6827712072110597)
 (1690, 0.01430815837450065)
 (1700, 0.014476627030865307)
 (1757, 0.06191521912803008)

jajhall · December 9, 2023, 8:55pm

The HiGHS QP solver (currently) is active set, not IPM, so its dual values may be far from feasible until an optimal solution is found.

Since your Hessian sounds as if it’s only semi-definite, it looks to me as if the QP solver is floating around a space of optimal primal solutions, looking to get dual feasibility. This may be why perturbations of the Hessian diagonal lead to faster termination.

An IPM-based QP solver will be much faster. Also, as @odow says, the HiGHS QP solver isn’t of the same quality as its LP simplex/barrier or MIP solvers.

bernardofpc · January 15, 2024, 12:14pm

I want to second on this bug, with a smaller model which (it seems to me) triggers the same problem. Here is the MWE I got (from a cutting-plane algorithm, which explains the strange coefficients):

using JuMP
using HiGHS, Ipopt

m = Model(Ipopt.Optimizer)
@variable(m, x >= 0)
@variable(m, θ >= 0)
@constraint(m, 803.3590017658552x + θ >= 40880.465738232284)
@constraint(m, 397.24681679972514x + θ >= 30547.566482175083)
@objective(m, Min, 3x^2 + θ)

optimize!(m)

Solving with Ipopt returns an optimal solution. The optimal value with Ipopt is 17397.14692738473, and the solution is

julia> value.([x, θ])
2-element Vector{Float64}:
   66.20780279993312
 4246.7274726002315

Replacing with HiGHS (and a time limit of 1 second) the solver goes on an infinite loop:

Iteration, Runtime, ObjVal, NullspaceDim
0, 0.001994, 17740.001017, 0
1000, 0.008975, 17740.001017, 0
[...]
106000, 0.972836, 17740.001017, 0
107000, 0.983806, 17740.001017, 0
108000, 0.992788, 17740.001017, 0
108563, 1.000763, 22382.341885, 0
Model   status      : Time limit reached
Objective value     :  2.2382341885e+04

Notice that the solver seems “trapped” on an “ObjVal” which is actually higher than the optimal value found in Ipopt. The solution returned is also suboptimal:

julia> value.([x, θ])
2-element Vector{Float64}:
    25.44346030129327
 20440.23286911615

bernardofpc · January 15, 2024, 12:23pm

Just to mention that the coefficients do not seem to really matter here. The problem

using JuMP, HiGHS, Ipopt

m = Model(HiGHS.Optimizer)
@variable(m, x >= 0)
@variable(m, θ >= 0)
@constraint(m, 800x + θ >= 40000)
@constraint(m, 400x + θ >= 30000)
@objective(m, Min, 3x^2 + θ)

set_time_limit_sec(m, 1)
optimize!(m)

also triggers the same bug.

The last lines from HiGHS (and the solution found) are

[...]
105000, 0.977208, 16875.000000, 0
106000, 0.987248, 16875.000000, 0
107000, 0.995227, 16875.000000, 0
107452, 1.000215, 16875.000000, 0
Model   status      : Time limit reached
Objective value     :  1.6875000000e+04
HiGHS run time      :          1.01

julia> value.([x, θ])
2-element Vector{Float64}:
 75.0
  0.0

Solving with Ipopt:

julia> objective_value(m)
16666.666566669177

julia> value.([x, θ])
2-element Vector{Float64}:
   66.6666666666343
 3333.33323334879

It seems the solver is stuck at one of the corners of the polyhedral function described by θ. Both 75 here (30000/400) and 25.443460301293275 before are intersections of the lines bounding θ from below.

I’d even risk saying it is “bouncing between” two vertices of the solution on θ, instead of finding the compromise between them (and reaching the optimum in the other coordinate x).

jajhall · January 15, 2024, 12:42pm

Thanks for deriving such a small instance that exposes this strange behaviour. I’ll set this up as a HiGHS issue independent of JuMP so that our QP developer can have a look.