SDP-relaxation in polynomial optimization

[kalmarek]

This is just a note for the documentation, namely the example here:
https://ericphanson.github.io/SDPAFamily.jl/dev/examples/#SDP-relaxation-in-polynomial-optimization-problem-1

The fact that SCS couldn't solve the problem could be the consequence of how Convex formulation. I just tried JuMPed version and it works just fine:

julia> using Jump, SCS;
julia> function relaxed_pop(r::Int)
                  m = Model()
                  v = @variable m v[1:2*r]
                  M1 = v[1:1+r]'
                  for i in 2:r
                      M1 = vcat(M1, v[i:i+r]')
                  end
                  t = [1 v[1:r]']
                  M1 = vcat(t, M1)
                  @constraint m M1 in PSDCone()
                  @constraint m M1[2:end, 1:end-1] in PSDCone()
                  @constraint m M1[2:end, 2:end] - M1[1:end-1, 1:end-1] in PSDCone()
                  return m
              end;

julia> m = relaxed_pop(5)
A JuMP Model
Feasibility problem with:
Variables: 10
`Array{GenericAffExpr{Float64,VariableRef},1}`-in-`MathOptInterface.PositiveSemidefiniteConeSquare`: 3 constraints
Model mode: AUTOMATIC
CachingOptimizer state: NO_OPTIMIZER
Solver name: No optimizer attached.
Names registered in the model: v

julia> set_optimizer(m, SCS.Optimizer)

julia> optimize!(m)
----------------------------------------------------------------------------
        SCS v2.1.1 - Splitting Conic Solver
        (c) Brendan O'Donoghue, Stanford University, 2012
----------------------------------------------------------------------------
Lin-sys: sparse-indirect, nnz in A = 64, CG tol ~ 1/iter^(2.00)
eps = 1.00e-05, alpha = 1.50, max_iters = 5000, normalize = 1, scale = 1.00
acceleration_lookback = 10, rho_x = 1.00e-03
Variables n = 10, constraints m = 51
Cones:  sd vars: 51, sd blks: 3
Setup time: 5.15e-03s
----------------------------------------------------------------------------
 Iter | pri res | dua res | rel gap | pri obj | dua obj | kap/tau | time (s)
----------------------------------------------------------------------------
     0| 1.13e+19  2.45e+13  1.00e+00  0.00e+00  4.24e+13  3.57e+13  6.73e-03 
   100| 5.09e-03  2.32e-08  1.33e-09  0.00e+00  1.33e-09  6.40e-23  1.48e-02 
   200| 5.55e-03  3.74e-08  7.50e-09  0.00e+00  7.50e-09  4.30e-24  2.03e-02 
   300| 1.98e-01  1.44e-06  9.96e-08  0.00e+00  9.96e-08  2.54e-23  2.58e-02 
   400| 3.09e-03  1.11e-08  1.00e-10  0.00e+00 -1.00e-10  8.10e-23  3.37e-02 
   500| 4.08e-03  2.93e-08  1.59e-09  0.00e+00  1.59e-09  8.51e-23  3.99e-02 
   600| 3.57e-02  3.34e-07  4.91e-08  0.00e+00  4.91e-08  2.35e-23  4.51e-02 
   700| 7.44e-04  6.29e-09  5.78e-10  0.00e+00 -5.78e-10  3.66e-23  5.05e-02 
   800| 6.82e-04  4.00e-09  1.35e-09  0.00e+00 -1.35e-09  4.67e-23  5.54e-02 
   900| 3.54e-03  1.56e-08  2.34e-09  0.00e+00  2.34e-09  2.33e-23  6.06e-02 
  1000| 3.04e-02  2.40e-07  9.32e-09  0.00e+00  9.32e-09  3.74e-23  6.45e-02 
  1100| 3.65e-03  3.46e-08  4.10e-09  0.00e+00  4.10e-09  3.49e-23  6.92e-02 
  1200| 9.83e-02  5.10e-07  5.31e-08  0.00e+00 -5.31e-08  2.79e-23  7.93e-02 
  1300| 7.09e-01  3.72e-06  1.96e-07  0.00e+00 -1.96e-07  3.90e-23  8.86e-02 
  1400| 3.62e-03  2.02e-08  6.16e-09  0.00e+00 -6.16e-09  4.36e-23  9.64e-02 
  1500| 3.12e-04  2.41e-09  2.19e-10  0.00e+00 -2.19e-10  1.01e-22  1.05e-01 
  1600| 3.05e-04  9.67e-10  1.01e-10  0.00e+00 -1.01e-10  2.87e-23  1.25e-01 
  1700| 1.75e-02  9.18e-08  8.22e-09  0.00e+00 -8.22e-09  5.27e-23  1.36e-01 
  1800| 1.01e-04  7.51e-10  2.97e-11  0.00e+00 -2.97e-11  3.52e-23  1.64e-01 
  1900| 1.37e-04  1.11e-09  2.87e-10  0.00e+00  2.87e-10  3.99e-23  1.74e-01 
  2000| 5.62e-03  2.71e-08  2.48e-09  0.00e+00 -2.48e-09  4.73e-23  1.92e-01 
  2100| 9.71e-03  5.09e-08  5.37e-09  0.00e+00 -5.37e-09  2.73e-23  2.08e-01 
  2200| 1.02e-04  5.56e-10  1.13e-10  0.00e+00  1.13e-10  6.61e-23  2.46e-01 
  2300| 1.19e-04  9.54e-10  1.09e-10  0.00e+00  1.09e-10  1.79e-23  2.54e-01 
  2400| 2.21e-04  1.06e-09  7.14e-11  0.00e+00 -7.14e-11  1.82e-23  2.81e-01 
  2500| 7.06e-05  6.86e-10  6.78e-11  0.00e+00 -6.78e-11  1.40e-23  2.89e-01 
  2600| 2.09e-02  1.12e-07  1.16e-08  0.00e+00 -1.16e-08  3.99e-23  3.04e-01 
  2700| 9.40e-01  1.88e-06  8.24e-08  0.00e+00 -8.24e-08  4.80e-24  3.12e-01 
  2800| 9.16e-05  6.87e-10  8.57e-11  0.00e+00 -8.57e-11  1.46e-23  3.25e-01 
  2900| 7.13e-05  3.54e-10  6.99e-11  0.00e+00  6.99e-11  2.50e-23  3.32e-01 
  3000| 1.84e-01  0.00e+00  0.00e+00  0.00e+00 -0.00e+00  4.17e-24  3.38e-01 
  3040| 4.84e-06  2.08e-11  6.16e-12  0.00e+00  6.16e-12  9.66e-23  3.40e-01 
----------------------------------------------------------------------------
Status: Solved
Timing: Solve time: 3.40e-01s
        Lin-sys: avg # CG iterations: 8.48, avg solve time: 9.66e-06s
        Cones: avg projection time: 7.89e-05s
        Acceleration: avg step time: 1.85e-05s
----------------------------------------------------------------------------
Error metrics:
dist(s, K) = 1.9366e-09, dist(y, K*) = 2.9910e-09, s'y/|s||y| = -4.1375e-11
primal res: |Ax + s - b|_2 / (1 + |b|_2) = 4.8386e-06
dual res:   |A'y + c|_2 / (1 + |c|_2) = 2.0795e-11
rel gap:    |c'x + b'y| / (1 + |c'x| + |b'y|) = 6.1619e-12
----------------------------------------------------------------------------
c'x = 0.0000, -b'y = 0.0000
============================================================================

julia> objective_value(m)
0.0

[ericphanson]

Interesting, thanks for the note; we should update the docs. The original paper (linked in the docs) used SDPA and SeDuMi, not SCS, although I'm not sure exactly how they formulated the problem.

With SDPA (via SDPA.jl), I find that r=5 gives the correct answer but r=6 gives warnings and the incorrect answer. However, SDPA.jl does not implement a presolve and MOI can sometimes emit linearly dependent constraints violating one of the assumptions of SDPA. However, even using the presolve via SDPAFamily I get the same warnings and with larger r, I get e.g.

julia> r=10
10

julia> m = relaxed_pop(r)
A JuMP Model
Feasibility problem with:
Variables: 20
`Array{GenericAffExpr{Float64,VariableRef},1}`-in-`MathOptInterface.PositiveSemidefiniteConeSquare`: 3 constraints
Model mode: AUTOMATIC
CachingOptimizer state: NO_OPTIMIZER
Solver name: No optimizer attached.
Names registered in the model: v

julia> set_optimizer(m, () -> SDPAFamily.Optimizer{Float64}(presolve=true, variant=:sdpa))

julia> optimize!(m)
Warning: Step length is too small.  :: line 198 in sdpa_dataset.cpp
Warning: cannot move: step length is too short :: line 176 in sdpa_solve.cpp

julia> objective_value(m)
0.0005971314478685485

(Note the above code calls plain SDPA, not SDPA-GMP).

It's not a silent failure though. I wonder if maybe the paper is just a bit out of date and the solvers can now more-or-less handle these problems? Or if maybe the paper authors also used some kind of non-optimal formulation like Convex emits here.

By the way, I suspect Convex's formulation is inefficient but not truly incorrect here (not that you were implying it was!), although it is definitely not good that it leads to the wrong answer in this case. I think the solution is to use MOI more directly in Convex as JuMP does instead of essentially assembling the MPB standard form and passing that to MOI. It'll take some work though.

[kalmarek]

I remember seeing different formulations of the same problem that performed differently with SCS, so that was my guess ;)

just tried with r up to 25, SCS gets numerically unstable (known thing with acceleration_lookback!=0), but does the job and delivers sound primal = 0.0 and dual <= 1e-11 in ~60000 iterations in under 2 minutes here.

[kalmarek]

btw I tried the same using MOI and SDPAFamily and hit

julia> let m = relaxed_pop(25);
           set_optimizer(m, SDPAFamily.Optimizer);
           optimize!(m)
           #@show termination_status(m)
           #objective_value(m)
       end
ERROR: TypeError: in typeassert, expected MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64},MathOptInterface.EqualTo{Float64}}, got MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{BigFloat},MathOptInterface.EqualTo{BigFloat}}
Stacktrace:
 [1] getindex at /home/kalmar/.julia/packages/MathOptInterface/RmalA/src/Utilities/copy.jl:75 [inlined]
 [2] #112 at /home/kalmar/.julia/packages/MathOptInterface/RmalA/src/Utilities/copy.jl:123 [inlined]
 [3] iterate at ./generator.jl:47 [inlined]
 [4] _collect(::Array{MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64},MathOptInterface.EqualTo{Float64}},1}, ::Base.Generator{Array{MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64},MathOptInterface.EqualTo{Float64}},1},MathOptInterface.Utilities.var"#112#113"{MathOptInterface.Utilities.IndexMap}}, ::Base.EltypeUnknown, ::Base.HasShape{1}) at ./array.jl:678
 [5] collect_similar at ./array.jl:607 [inlined]
 [6] map at ./abstractarray.jl:2072 [inlined]
 [7] pass_attributes(::SDPAFamily.Optimizer{BigFloat}, ::MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}}, ::Bool, ::MathOptInterface.Utilities.IndexMap, ::Array{MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64},MathOptInterface.EqualTo{Float64}},1}, ::Function) at /home/kalmar/.julia/packages/MathOptInterface/RmalA/src/Utilities/copy.jl:123
 [8] pass_constraints(::SDPAFamily.Optimizer{BigFloat}, ::MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}}, ::Bool, ::MathOptInterface.Utilities.IndexMap, ::Array{DataType,1}, ::Array{Array{#s118,1} where #s118<:(MathOptInterface.ConstraintIndex{MathOptInterface.SingleVariable,S} where S),1}, ::Array{DataType,1}, ::Array{Array{#s13,1} where #s13<:(MathOptInterface.ConstraintIndex{MathOptInterface.VectorOfVariables,S} where S),1}, ::typeof(MathOptInterface.Utilities.allocate_constraints), ::typeof(MathOptInterface.Utilities.allocate)) at /home/kalmar/.julia/packages/MathOptInterface/RmalA/src/Utilities/copy.jl:266
 [9] allocate_load(::SDPAFamily.Optimizer{BigFloat}, ::MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}}, ::Bool) at /home/kalmar/.julia/packages/MathOptInterface/RmalA/src/Utilities/copy.jl:684
 [10] #automatic_copy_to#109 at /home/kalmar/.julia/packages/MathOptInterface/RmalA/src/Utilities/copy.jl:17 [inlined]
 [11] #copy_to#22 at /home/kalmar/.julia/packages/SDPAFamily/rG0xe/src/MOI_wrapper.jl:211 [inlined]
 [12] attach_optimizer(::MathOptInterface.Utilities.CachingOptimizer{SDPAFamily.Optimizer{BigFloat},MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}}}) at /home/kalmar/.julia/packages/MathOptInterface/RmalA/src/Utilities/cachingoptimizer.jl:149
 [13] optimize!(::MathOptInterface.Utilities.CachingOptimizer{SDPAFamily.Optimizer{BigFloat},MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}}}) at /home/kalmar/.julia/packages/MathOptInterface/RmalA/src/Utilities/cachingoptimizer.jl:185
 [14] optimize!(::MathOptInterface.Bridges.LazyBridgeOptimizer{MathOptInterface.Utilities.CachingOptimizer{SDPAFamily.Optimizer{BigFloat},MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}}}}) at /home/kalmar/.julia/packages/MathOptInterface/RmalA/src/Bridges/bridge_optimizer.jl:239
 [15] optimize!(::MathOptInterface.Utilities.CachingOptimizer{MathOptInterface.AbstractOptimizer,MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}}}) at /home/kalmar/.julia/packages/MathOptInterface/RmalA/src/Utilities/cachingoptimizer.jl:189
 [16] optimize!(::Model, ::Nothing; bridge_constraints::Bool, ignore_optimize_hook::Bool, kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at /home/kalmar/.julia/packages/JuMP/CZ8vV/src/optimizer_interface.jl:131
 [17] optimize! at /home/kalmar/.julia/packages/JuMP/CZ8vV/src/optimizer_interface.jl:107 [inlined] (repeats 2 times)
 [18] top-level scope at REPL[10]:3

[ericphanson]

Ah we should figure out how to throw a better error message for that. The fix is to use SDPAFamily.Optimizer{Float64} instead; it defaults to BigFloat since it’s a high precision solver but JuMP can only treat Float64’s (so far).

[kalmarek]

indeed that's a tricky question

[ericphanson]

The simplest thing actually might be to just not choose a default at all and throw a clear error message when one passes just SDPA.Optimizer telling the user to choose a numeric type.

Maybe a nicer solution is if MOI can support numeric promotion (some kind of MOI.ScalarAffineFunction{Float64} to MOI.ScalarAffineFunction{BigFloat} bridge), but actually we will end up needing to bridge both directions to load the solution back after solving and that's not lossless.

It's a little unfortunate that SDPAFamily.Optimizer needs to choose a numeric type and the modelling language's model (either Convex or JuMP, the latter of which is hard-coded to Float64) also needs to choose a numeric type, and these choices must agree but one cannot determine the other. Even in Convex, one has to do problem = minimize(objective; numeric_type = BigFloat) for example, to use the right model. Ideally then, if Convex already knows it is using a BigFloat model, it could pass that info along to the optimizer, but there doesn't really seem to be a generic way to do that.