-
You don’t need using MathOptInterface
. It is already exported by JuMP.
-
The looks like a bug. We should be able to bridge the LogDet cone into PSD + Exponential. (Edit: here’s the issue: https://github.com/jump-dev/MathOptInterface.jl/issues/541)
Rewriting to use a LogDetConeTriangle
, here is what I get with SCS:
julia> using JuMP, SCS
julia> u = rand(3)
3-element Vector{Float64}:
0.2076471595685787
0.6577478008368793
0.4799996513530598
julia> M_SDP = Model(SCS.Optimizer)
A JuMP Model
Feasibility problem with:
Variables: 0
Model mode: AUTOMATIC
CachingOptimizer state: EMPTY_OPTIMIZER
Solver name: SCS
julia> @variable(M_SDP, -1 <= buro[1:3] <= 1)
3-element Vector{VariableRef}:
buro[1]
buro[2]
buro[3]
julia> @variable(M_SDP, P[1:3, 1:3] >= 0, PSD)
3×3 Symmetric{VariableRef, Matrix{VariableRef}}:
P[1,1] P[1,2] P[1,3]
P[1,2] P[2,2] P[2,3]
P[1,3] P[2,3] P[3,3]
julia> @constraint(M_SDP, [1; P * u + buro] in SecondOrderCone())
[1, buro[1] + 0.2076471595685787 P[1,1] + 0.6577478008368793 P[1,2] + 0.4799996513530598 P[1,3], buro[2] + 0.2076471595685787 P[1,2] + 0.6577478008368793 P[2,2] + 0.4799996513530598 P[2,3], buro[3] + 0.2076471595685787 P[1,3] + 0.6577478008368793 P[2,3] + 0.4799996513530598 P[3,3]] ∈ MathOptInterface.SecondOrderCone(4)
julia> @variable(M_SDP, t)
t
julia> lower_triangular(P) = [P[i, j] for i = 1:size(P, 1) for j = 1:i]
lower_triangular (generic function with 1 method)
julia> @constraint(M_SDP, [t; 1; lower_triangular(P)] in MOI.LogDetConeTriangle(3))
[t, 1, P[1,1], P[1,2], P[2,2], P[1,3], P[2,3], P[3,3]] ∈ MathOptInterface.LogDetConeTriangle(3)
julia> @objective(M_SDP, Max, t)
t
julia> optimize!(M_SDP)
----------------------------------------------------------------------------
SCS v2.1.3 - Splitting Conic Solver
(c) Brendan O'Donoghue, Stanford University, 2012
----------------------------------------------------------------------------
Lin-sys: sparse-indirect, nnz in A = 55, 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 = 19, constraints m = 53
Cones: linear vars: 13
soc vars: 4, soc blks: 1
sd vars: 27, sd blks: 2
exp vars: 9, dual exp vars: 0
Setup time: 6.76e-05s
----------------------------------------------------------------------------
Iter | pri res | dua res | rel gap | pri obj | dua obj | kap/tau | time (s)
----------------------------------------------------------------------------
0| 2.81e+19 2.40e+19 1.00e+00 -1.26e+20 4.84e+19 1.55e+20 7.41e-05
80| 3.56e-06 6.66e-06 6.04e-07 -4.09e+00 -4.09e+00 3.50e-17 4.42e-03
----------------------------------------------------------------------------
Status: Solved
Timing: Solve time: 4.43e-03s
Lin-sys: avg # CG iterations: 4.86, avg solve time: 1.66e-06s
Cones: avg projection time: 4.55e-05s
Acceleration: avg step time: 5.91e-06s
----------------------------------------------------------------------------
Error metrics:
dist(s, K) = 2.1945e-09, dist(y, K*) = 1.7850e-09, s'y/|s||y| = 5.9406e-12
primal res: |Ax + s - b|_2 / (1 + |b|_2) = 3.5589e-06
dual res: |A'y + c|_2 / (1 + |c|_2) = 6.6559e-06
rel gap: |c'x + b'y| / (1 + |c'x| + |b'y|) = 6.0362e-07
----------------------------------------------------------------------------
c'x = -4.0921, -b'y = -4.0921
============================================================================
julia> -objective_value(M_SDP), log(det(inv(value.(P))))
(-4.092065768181431, -4.09206772307858)