v1.4.0 is giving incorrect solutions

[kamilkhanlab]

After upgrading PATHSolver.jl from v1.3.0 to v1.4.0, it now seems to produce incorrect solutions via JuMP v1.6.0. We have seen consistently incorrect results across a few setups: a 2016 MacBook Pro (Intel), a more recent MacBook Pro (M1), and a Windows laptop. If we revert PATHSolver.jl back to v1.3.0, then we get correct solutions again.

Here are summaries of two examples illustrating this behavior. The corresponding REPL output for each example is pasted subsequently.

1. When attempting to replicate the LCP example in PATHSolver.jl's README, PATHSolver.jl v1.4.0 instead reports the incorrect solution: x = [6.0, 0.0, 0.0, 0.0].

2. When attempting to solve the simple MLCP: -1.0 + 2.0*u ⟂ u with u unbounded, PATHSolver.jl v1.4.0 reports the solution as u = 0.0, which is incorrect. (Since u is otherwise unbounded, this MLCP should instead be equivalent to the equation -1.0 + 2.0*u == 0.0, by my reading of JuMP's documentation on complementarity constraints.)

Here is the corresponding REPL output for Example 1 above (attempting to replicate the LCP in the README):

julia> using JuMP, PATHSolver

julia> M = [
            0 0 -1 -1
            0 0 1 -2
            1 -1 2 -2
            1 2 -2 4
        ]
4×4 Matrix{Int64}:
 0   0  -1  -1
 0   0   1  -2
 1  -1   2  -2
 1   2  -2   4

julia> q = [2, 2, -2, -6]
4-element Vector{Int64}:
  2
  2
 -2
 -6

julia> model = Model(PATHSolver.Optimizer)
A JuMP Model
Feasibility problem with:
Variables: 0
Model mode: AUTOMATIC
CachingOptimizer state: EMPTY_OPTIMIZER
Solver name: Path 5.0.03

julia> set_optimizer_attribute(model, "output", "no")

julia> @variable(model, x[1:4] >= 0)
4-element Vector{VariableRef}:
 x[1]
 x[2]
 x[3]
 x[4]

julia> @constraint(model, M * x .+ q ⟂ x)
[-x[3] - x[4] + 2, x[3] - 2 x[4] + 2, x[1] - x[2] + 2 x[3] - 2 x[4] - 2, x[1] + 2 x[2] - 2 x[3] + 4 x[4] - 6, x[1], x[2], x[3], x[4]] ∈ MathOptInterface.Complements(8)

julia> optimize!(model)
Reading options file /var/folders/sw/8x_srljj1wn1njw0jsp4vyqm0000gn/T/jl_T2HyeD
Read of options file complete.

Path 5.0.03 (Fri Jun 26 09:58:07 2020)
Written by Todd Munson, Steven Dirkse, Youngdae Kim, and Michael Ferris

Zero:     1 Single:     0 Double:     3
Postsolved residual: 1.6754e+00

julia> value.(x)
4-element Vector{Float64}:
 6.0
 0.0
 0.0
 0.0

julia> termination_status(model)
LOCALLY_SOLVED::TerminationStatusCode = 4

And here is the REPL output for Example 2 above (a simple MCP):

julia> import JuMP, PATHSolver

julia> optimizer = PATHSolver.Optimizer
MathOptInterface.Utilities.GenericOptimizer{T, MathOptInterface.Utilities.ObjectiveContainer{T}, MathOptInterface.Utilities.VariablesContainer{T}, MathOptInterface.Utilities.VectorOfConstraints{MathOptInterface.VectorAffineFunction{T}, MathOptInterface.Complements}} where T

julia> solverAttributes = ("output" => "no",)
("output" => "no",)

julia> model = JuMP.Model(JuMP.optimizer_with_attributes(optimizer, solverAttributes...))
A JuMP Model
Feasibility problem with:
Variables: 0
Model mode: AUTOMATIC
CachingOptimizer state: EMPTY_OPTIMIZER
Solver name: Path 5.0.03

julia> JuMP.@variable(model, u)
u

julia> JuMP.@constraint(model, complements(-1.0 + 2.0*u, u))
[2 u - 1, u] ∈ MathOptInterface.Complements(2)

julia> JuMP.optimize!(model)
Reading options file /var/folders/sw/8x_srljj1wn1njw0jsp4vyqm0000gn/T/jl_5O5srS
Read of options file complete.

Path 5.0.03 (Fri Jun 26 09:58:07 2020)
Written by Todd Munson, Steven Dirkse, Youngdae Kim, and Michael Ferris

MCPR: One:    0 Two:    0 Thr:    0 Qua:    0 WBd:    0 Fix:    1 IFx:    0
Postsolved residual: 1.0000e+00
julia> uStar = JuMP.value.(u)
-0.0

julia> terminationStatus = JuMP.termination_status(model)
LOCALLY_SOLVED::TerminationStatusCode = 4

[odow]

Let me take a look

[odow]

It's to do with the starting point.

The quickest work-around is to explicitly set a valid starting point:

julia> model = Model(PATHSolver.Optimizer)
A JuMP Model
Feasibility problem with:
Variables: 0
Model mode: AUTOMATIC
CachingOptimizer state: EMPTY_OPTIMIZER
Solver name: Path 5.0.03

julia> set_optimizer_attribute(model, "output", "no")

julia> @variable(model, x[1:4] >= 0, start = 0)
4-element Vector{VariableRef}:
 x[1]
 x[2]
 x[3]
 x[4]

julia> @constraint(model, M * x .+ q ⟂ x)
[-x[3] - x[4] + 2, x[3] - 2 x[4] + 2, x[1] - x[2] + 2 x[3] - 2 x[4] - 2, x[1] + 2 x[2] - 2 x[3] + 4 x[4] - 6, x[1], x[2], x[3], x[4]] ∈ MathOptInterface.Complements(8)

julia> optimize!(model)
(M, q, lower, upper, initial) = (
  ⋅     ⋅   -1.0  -1.0
  ⋅     ⋅    1.0  -2.0
 1.0  -1.0   2.0  -2.0
 1.0   2.0  -2.0   4.0, [2.0, 2.0, -2.0, -6.0], [0.0, 0.0, 0.0, 0.0], [Inf, Inf, Inf, Inf], [0.0, 0.0, 0.0, 0.0])
Reading options file /var/folders/bg/dzq_hhvx1dxgy6gb5510pxj80000gn/T/jl_tgiBJB
Read of options file complete.

Path 5.0.03 (Fri Jun 26 09:58:07 2020)
Written by Todd Munson, Steven Dirkse, Youngdae Kim, and Michael Ferris

julia> value.(x)
4-element Vector{Float64}:
 2.8
 0.0
 0.8
 1.2

I'll make a PR with a better fix.

[odow]

Fix incoming: https://github.com/chkwon/PATHSolver.jl/pull/81.

Should be released in an hour or so. Once a new release has been tagged, https://github.com/chkwon/PATHSolver.jl/releases, try updating your packages. You'll be after HiGHS.jl v1.4.1.