Failure when solving the Barrier certificates example from the Tutorials using MosekTools instead of CSDP

[hurak]

I tried to use MosekTools instead of CSDP as the solver when solving the barrier certificates example, but the code fails. I am not able to decipher if SumOfSquares.jl or MosekTools.jl or Mosek is the culprit.

using SumOfSquares
using DynamicPolynomials
using MosekTools

optimizer = optimizer_with_attributes(MosekTools.Optimizer, MOI.Silent() => true)
model = SOSModel(optimizer)
@polyvar x[1:2] 

p = 1;

f = [ x[2],
     -x[1] + (p/3)*x[1]^3 - x[2]]

g₁ = -(x[1]+1)^2 - (x[2]+1)^2 + 0.16  # 𝒳ᵤ = {x ∈ R²: g₁(x) ≥ 0}
h₁ = -(x[1]-1.5)^2 - x[2]^2 + 0.25    # 𝒳₀ = {x ∈ R²: h₁(x) ≥ 0}

X = monomials(x, 0:4)
@variable(model, B, Poly(X))

ε = 0.001
@constraint(model, B >= ε, domain = @set(g₁ >= 0))
@constraint(model, B <= 0, domain = @set(h₁ >= 0))

using LinearAlgebra # Needed for `dot`
dBdt = dot(differentiate(B, x), f)
@constraint(model, -dBdt >= 0)

optimize!(model)

ERROR: Mosek.MosekError(3944, "The value (0,0) in a sparse symmetric has been specified twice at position 0 and position 1.\n")
Stacktrace:
  [1] appendsparsesymmat(task::Mosek.Task, dim::Int32, subi::Vector{Int32}, subj::Vector{Int32}, valij::Vector{Float64})
    @ Mosek ~/.julia/packages/Mosek/cMwUF/src/msk_functions.jl:3707
  [2] appendsparsesymmat(task::Mosek.Task, dim::Int64, subi::Vector{Int32}, subj::Vector{Int32}, valij::Vector{Float64})
    @ Mosek ~/.julia/packages/Mosek/cMwUF/src/msk_functions.jl:16628
  [3] (::MosekTools.var"#add_sd#24"{MosekTools.Optimizer, MosekTools.var"#25#26"{…}, Vector{…}, Vector{…}, Vector{…}})()
    @ MosekTools ~/.julia/packages/MosekTools/uICAS/src/constraint.jl:69
  [4] split_scalar_matrix(m::MosekTools.Optimizer, terms::Vector{…}, set_sd::MosekTools.var"#25#26"{…})
    @ MosekTools ~/.julia/packages/MosekTools/uICAS/src/constraint.jl:93
  [5] set_row
    @ ~/.julia/packages/MosekTools/uICAS/src/constraint.jl:103 [inlined]
  [6] add_constraint(m::MosekTools.Optimizer, axb::MathOptInterface.ScalarAffineFunction{…}, dom::MathOptInterface.EqualTo{…})
    @ MosekTools ~/.julia/packages/MosekTools/uICAS/src/constraint.jl:437
  [7] add_constraint
    @ ~/.julia/packages/MathOptInterface/gLl4d/src/Bridges/bridge_optimizer.jl:1905 [inlined]
  [8] #normalize_and_add_constraint#51
    @ ~/.julia/packages/MathOptInterface/gLl4d/src/Utilities/constraints.jl:25 [inlined]
  [9] normalize_and_add_constraint
    @ ~/.julia/packages/MathOptInterface/gLl4d/src/Utilities/constraints.jl:19 [inlined]
 [10] bridge_constraint(::Type{…}, model::MathOptInterface.Bridges.LazyBridgeOptimizer{…}, f::MathOptInterface.VectorAffineFunction{…}, set::MathOptInterface.Zeros)
    @ MathOptInterface.Bridges.Constraint ~/.julia/packages/MathOptInterface/gLl4d/src/Bridges/Constraint/bridges/scalarize.jl:46
 [11] add_bridged_constraint(b::MathOptInterface.Bridges.LazyBridgeOptimizer{…}, BridgeType::Type, f::MathOptInterface.VectorAffineFunction{…}, s::MathOptInterface.Zeros)
    @ MathOptInterface.Bridges ~/.julia/packages/MathOptInterface/gLl4d/src/Bridges/bridge_optimizer.jl:1773
 [12] add_constraint(b::MathOptInterface.Bridges.LazyBridgeOptimizer{…}, f::MathOptInterface.VectorAffineFunction{…}, s::MathOptInterface.Zeros)
    @ MathOptInterface.Bridges ~/.julia/packages/MathOptInterface/gLl4d/src/Bridges/bridge_optimizer.jl:1903
 [13] bridge_constraint(::Type{…}, model::MathOptInterface.Bridges.LazyBridgeOptimizer{…}, f::MathOptInterface.VectorAffineFunction{…}, s::PolyJuMP.ZeroPolynomialSet{…})
    @ PolyJuMP.Bridges.Constraint ~/.julia/packages/PolyJuMP/HbV7j/src/Bridges/Constraint/zero_polynomial.jl:18
 [14] add_bridged_constraint(b::MathOptInterface.Bridges.LazyBridgeOptimizer{…}, BridgeType::Type, f::MathOptInterface.VectorAffineFunction{…}, s::PolyJuMP.ZeroPolynomialSet{…})
    @ MathOptInterface.Bridges ~/.julia/packages/MathOptInterface/gLl4d/src/Bridges/bridge_optimizer.jl:1773
 [15] add_constraint(b::MathOptInterface.Bridges.LazyBridgeOptimizer{…}, f::MathOptInterface.VectorAffineFunction{…}, s::PolyJuMP.ZeroPolynomialSet{…})
    @ MathOptInterface.Bridges ~/.julia/packages/MathOptInterface/gLl4d/src/Bridges/bridge_optimizer.jl:1903
 [16] bridge_constraint(::Type{…}, model::MathOptInterface.Bridges.LazyBridgeOptimizer{…}, f::MathOptInterface.VectorAffineFunction{…}, s::SumOfSquares.SOSPolynomialSet{…})
    @ SumOfSquares.Bridges.Constraint ~/.julia/packages/SumOfSquares/xNLIR/src/Bridges/Constraint/sos_polynomial.jl:68
 [17] add_bridged_constraint(b::MathOptInterface.Bridges.LazyBridgeOptimizer{…}, BridgeType::Type, f::MathOptInterface.VectorAffineFunction{…}, s::SumOfSquares.SOSPolynomialSet{…})
    @ MathOptInterface.Bridges ~/.julia/packages/MathOptInterface/gLl4d/src/Bridges/bridge_optimizer.jl:1773
 [18] add_constraint(b::MathOptInterface.Bridges.LazyBridgeOptimizer{…}, f::MathOptInterface.VectorAffineFunction{…}, s::SumOfSquares.SOSPolynomialSet{…})
    @ MathOptInterface.Bridges ~/.julia/packages/MathOptInterface/gLl4d/src/Bridges/bridge_optimizer.jl:1903
 [19] bridge_constraint(::Type{…}, model::MathOptInterface.Bridges.LazyBridgeOptimizer{…}, f::MathOptInterface.VectorAffineFunction{…}, set::SumOfSquares.SOSPolynomialSet{…})
    @ SumOfSquares.Bridges.Constraint ~/.julia/packages/SumOfSquares/xNLIR/src/Bridges/Constraint/sos_polynomial_in_semialgebraic_set.jl:86
 [20] add_bridged_constraint(b::MathOptInterface.Bridges.LazyBridgeOptimizer{…}, BridgeType::Type, f::MathOptInterface.VectorAffineFunction{…}, s::SumOfSquares.SOSPolynomialSet{…})
    @ MathOptInterface.Bridges ~/.julia/packages/MathOptInterface/gLl4d/src/Bridges/bridge_optimizer.jl:1773
 [21] add_constraint(b::MathOptInterface.Bridges.LazyBridgeOptimizer{…}, f::MathOptInterface.VectorAffineFunction{…}, s::SumOfSquares.SOSPolynomialSet{…})
    @ MathOptInterface.Bridges ~/.julia/packages/MathOptInterface/gLl4d/src/Bridges/bridge_optimizer.jl:1903
 [22] _copy_constraints(dest::MathOptInterface.Bridges.LazyBridgeOptimizer{…}, src::MathOptInterface.Utilities.UniversalFallback{…}, index_map::MathOptInterface.Utilities.IndexMap, index_map_FS::MathOptInterface.Utilities.DoubleDicts.IndexDoubleDictInner{…}, cis_src::Vector{…})
    @ MathOptInterface.Utilities ~/.julia/packages/MathOptInterface/gLl4d/src/Utilities/copy.jl:180
 [23] _copy_constraints(dest::MathOptInterface.Bridges.LazyBridgeOptimizer{…}, src::MathOptInterface.Utilities.UniversalFallback{…}, index_map::MathOptInterface.Utilities.IndexMap, cis_src::Vector{…})
    @ MathOptInterface.Utilities ~/.julia/packages/MathOptInterface/gLl4d/src/Utilities/copy.jl:192
 [24] pass_nonvariable_constraints_fallback(dest::MathOptInterface.Bridges.LazyBridgeOptimizer{…}, src::MathOptInterface.Utilities.UniversalFallback{…}, index_map::MathOptInterface.Utilities.IndexMap, constraint_types::Vector{…})
    @ MathOptInterface.Utilities ~/.julia/packages/MathOptInterface/gLl4d/src/Utilities/copy.jl:203
 [25] pass_nonvariable_constraints(dest::MathOptInterface.Bridges.LazyBridgeOptimizer{…}, src::MathOptInterface.Utilities.UniversalFallback{…}, idxmap::MathOptInterface.Utilities.IndexMap, constraint_types::Vector{…})
    @ MathOptInterface.Bridges ~/.julia/packages/MathOptInterface/gLl4d/src/Bridges/bridge_optimizer.jl:0
 [26] _pass_constraints(dest::MathOptInterface.Bridges.LazyBridgeOptimizer{…}, src::MathOptInterface.Utilities.UniversalFallback{…}, index_map::MathOptInterface.Utilities.IndexMap, variable_constraints_not_added::Vector{…})
    @ MathOptInterface.Utilities ~/.julia/packages/MathOptInterface/gLl4d/src/Utilities/copy.jl:251
 [27] default_copy_to(dest::MathOptInterface.Bridges.LazyBridgeOptimizer{…}, src::MathOptInterface.Utilities.UniversalFallback{…})
    @ MathOptInterface.Utilities ~/.julia/packages/MathOptInterface/gLl4d/src/Utilities/copy.jl:393
 [28] copy_to
    @ ~/.julia/packages/MathOptInterface/gLl4d/src/Bridges/bridge_optimizer.jl:442 [inlined]
 [29] optimize!
    @ ~/.julia/packages/MathOptInterface/gLl4d/src/MathOptInterface.jl:121 [inlined]
 [30] optimize!(m::MathOptInterface.Utilities.CachingOptimizer{…})
    @ MathOptInterface.Utilities ~/.julia/packages/MathOptInterface/gLl4d/src/Utilities/cachingoptimizer.jl:321
 [31] optimize!(model::Model; ignore_optimize_hook::Bool, _differentiation_backend::MathOptInterface.Nonlinear.SparseReverseMode, kwargs::@Kwargs{})
    @ JuMP ~/.julia/packages/JuMP/i68GU/src/optimizer_interface.jl:595
 [32] optimize!(model::Model)
    @ JuMP ~/.julia/packages/JuMP/i68GU/src/optimizer_interface.jl:546
 [33] top-level scope
    @ ~/.julia/dev/Jhys.jl/src/barrier.jl:30
Some type information was truncated. Use `show(err)` to see complete types.

[blegat]

That's suspicious indeed, thanks for reporting it

[hurak]

It appears that the example works fine even with Mosek(Tools.jl). I do not know if any of the related packages (MosekTools v0.15.4, SumOfSquares v0.7.3, ...) has been updated in the meantime, but the optimization now runs as expected. Perhaps it is possible to close the issue now.

[blegat]

Yes, it was a but in MosekTools that has been fixed in MosekTools v0.15.4