Bug in dual of PositiveSemidefiniteConeSquare

[odow]

Reported on Discourse: https://discourse.julialang.org/t/rank1-psd-constraint-dual-variables/80348

I guess this is in the bridges somewhere.

julia> using JuMP, SCS, LinearAlgebra

julia> function main(use_symmetric::Bool)
           n = 3
           model = Model(SCS.Optimizer)
           @variable(model, x[1:n] >= 0)
           @variable(model, X[1:n, 1:n], Symmetric)
           A = use_symmetric ? Symmetric([1 x'; x X]) : [1 x'; x X]
           c = @constraint(model, A in PSDCone())
           @objective(model, Min, sum(X) + sum(x))
           optimize!(model)
           return dual(c)
       end
main (generic function with 1 method)

julia> main(true)
------------------------------------------------------------------
	       SCS v3.2.0 - Splitting Conic Solver
	(c) Brendan O'Donoghue, Stanford University, 2012
------------------------------------------------------------------
problem:  variables n: 9, constraints m: 13
cones: 	  l: linear vars: 3
	  s: psd vars: 10, ssize: 1
settings: eps_abs: 1.0e-04, eps_rel: 1.0e-04, eps_infeas: 1.0e-07
	  alpha: 1.50, scale: 1.00e-01, adaptive_scale: 1
	  max_iters: 100000, normalize: 1, rho_x: 1.00e-06
	  acceleration_lookback: 10, acceleration_interval: 10
lin-sys:  sparse-direct
	  nnz(A): 12, nnz(P): 0
------------------------------------------------------------------
 iter | pri res | dua res |   gap   |   obj   |  scale  | time (s)
------------------------------------------------------------------
     0| 3.60e+01  2.00e+00  2.49e+02 -1.25e+02  1.00e-01  4.09e-05 
   125| 8.20e-06  5.15e-06  1.36e-05 -6.81e-06  1.47e+00  4.65e-04 
------------------------------------------------------------------
status:  solved
timings: total: 5.75e-04s = setup: 1.05e-04s + solve: 4.70e-04s
	 lin-sys: 4.02e-05s, cones: 3.43e-04s, accel: 7.41e-06s
------------------------------------------------------------------
objective = -0.000007
------------------------------------------------------------------
4×4 Symmetric{Float64, Matrix{Float64}}:
 1.20006e-9  3.46419e-5  3.46419e-5  3.46419e-5
 3.46419e-5  1.0         1.0         1.0
 3.46419e-5  1.0         1.0         1.0
 3.46419e-5  1.0         1.0         1.0

julia> main(false)
------------------------------------------------------------------
	       SCS v3.2.0 - Splitting Conic Solver
	(c) Brendan O'Donoghue, Stanford University, 2012
------------------------------------------------------------------
problem:  variables n: 9, constraints m: 13
cones: 	  l: linear vars: 3
	  s: psd vars: 10, ssize: 1
settings: eps_abs: 1.0e-04, eps_rel: 1.0e-04, eps_infeas: 1.0e-07
	  alpha: 1.50, scale: 1.00e-01, adaptive_scale: 1
	  max_iters: 100000, normalize: 1, rho_x: 1.00e-06
	  acceleration_lookback: 10, acceleration_interval: 10
lin-sys:  sparse-direct
	  nnz(A): 12, nnz(P): 0
------------------------------------------------------------------
 iter | pri res | dua res |   gap   |   obj   |  scale  | time (s)
------------------------------------------------------------------
     0| 3.60e+01  2.00e+00  2.49e+02 -1.25e+02  1.00e-01  5.40e-05 
   125| 8.20e-06  5.15e-06  1.36e-05 -6.81e-06  1.47e+00  6.67e-04 
------------------------------------------------------------------
status:  solved
timings: total: 8.50e-04s = setup: 1.79e-04s + solve: 6.71e-04s
	 lin-sys: 5.82e-05s, cones: 4.90e-04s, accel: 1.05e-05s
------------------------------------------------------------------
objective = -0.000007
------------------------------------------------------------------
4×4 Matrix{Float64}:
 1.20006e-9  6.92838e-5  6.92838e-5  6.92838e-5
 0.0         1.0         2.0         2.0
 0.0         0.0         1.0         2.0
 0.0         0.0         0.0         1.0

If you use return MOI.get(backend(model), MOI.ConstraintDual(), index(c)), it returns:

16-element Vector{Float64}:
 1.2000592582231642e-9
 0.0
 0.0
 0.0
 6.928375571224489e-5
 0.9999999999374926
 0.0
 0.0
 6.928375571224486e-5
 1.9999999998749844
 0.9999999999374926
 0.0
 6.92837557122449e-5
 1.999999999874985
 1.9999999998749844
 0.9999999999374923

so this is not a problem in JuMP due to reshaping.

[blegat]

This is expected behavior, see e.g.,

MathOptInterface.jl/src/Bridges/Constraint/square.jl

Lines 18 to 19 in 863f74b

	# Since A+Aᵀ is symmetric and we know that S₊ is self-dual, we have shown that (P₊)*
	# is the set of matrices A such that A+Aᵀ is PSD

To get the dual with bridges you should multiply by the adjoint of the transformation. See

MathOptInterface.jl/src/Bridges/Constraint/set_map.jl

Line 174 in 863f74b

return MOIB.adjoint_map_function(typeof(bridge), value)

so the upper triangular matrix is mathematically the correct answer.

[odow]

This is expected behavior

But I returns something in JuMP that is incorrect?

[blegat]

It's not incorrect. The dual A of the PSDSquare should be such that A + A' is PSD so you should take a look at A + A', not at A directly.
Anyway, I suggested something in #1839 (comment) something that should be less suprising for users not realizing that they create a PSDSquare constraint (most user) and those not knowing what the dual of this weird cone is (even more users).