Incorrect dual variable extracted from equality constraint between Hermitian matrices

[araujoms]

When I use a constraint of the form @constraint(model, A == B), where A and B are Hermitian matrices, for some reason the dual variable associated to this constraint has the off-diagonal elements multiplied by 2. I've tested a couple of different solvers, so the bug must be inside JuMP itself. A MWE follows below. Note the text-to-last line, where I'm dividing the off-diagonal elements of the dual variable by 2 to make it match the primal result.

using JuMP
using LinearAlgebra
import Hypatia
import Clarabel
import Ket
import Random

function dual_test(d)
    Random.seed!(1)
    ρ = Ket.random_state(d^2, 1)

    model = Model()
    @variable(model, σ[1:d^2, 1:d^2] in HermitianPSDCone())
    @variable(model, λ)
    noisy_state = Hermitian(ρ + λ * I(d^2))
    @constraint(model, witness_constraint, σ == noisy_state)
    @objective(model, Min, λ)
    @constraint(model, Ket.partial_transpose(σ, 2, [d, d]) in HermitianPSDCone())

    set_optimizer(model, Hypatia.Optimizer)
    #set_optimizer(model, Clarabel.Optimizer)
    optimize!(model)

    W = dual(witness_constraint)
    W = Hermitian(Diagonal(W) + 0.5(W - Diagonal(W))) #workaround
    return objective_value(model), dot(ρ, W)
end

[odow]

Moved this to MOI because I assume it's an issue in the bridge

[odow]

Oh wait. Hypatia supports this directly, so there is no bridge.

[odow]

Using

using JuMP
using LinearAlgebra
import Hypatia
import Clarabel
import Ket
import Random

function dual_test(d)
    Random.seed!(1)
    ρ = Ket.random_state(d^2, 1)
    model = Model(Hypatia.Optimizer; add_bridges = false)
    @variable(model, σ[1:d^2, 1:d^2] in HermitianPSDCone())
    @variable(model, λ)
    noisy_state = Hermitian(ρ + λ * I(d^2))
    @constraint(model, witness_constraint, σ == noisy_state)
    @objective(model, Min, λ)
    @constraint(model, Ket.partial_transpose(σ, 2, [d, d]) in HermitianPSDCone())
    optimize!(model)
    W = dual(witness_constraint)
    display(W)
    display(MOI.get(backend(model), MOI.ConstraintDual(), index(witness_constraint)))
    FIX = 0.5 * Hermitian(W + Diagonal(W))
    return objective_value(model), dot(ρ, FIX)
end
dual_test(2)

I get

4×4 Hermitian{ComplexF64, Matrix{ComplexF64}}:
 -0.406623+0.0im       -0.155587-0.357309im   0.353364-0.16435im    0.139813+0.123812im
 -0.155587+0.357309im  -0.093377+0.0im       0.0200512+0.812999im  -0.353364+0.16435im
  0.353364+0.16435im   0.0200512-0.812999im  -0.093377+0.0im        0.155587+0.357309im
  0.139813-0.123812im  -0.353364-0.16435im    0.155587-0.357309im  -0.406623+0.0im
16-element Vector{Float64}:
 -0.40662298907232036
 -0.15558743348537818
 -0.0933770120619247
  0.35336393807275995
  0.020051174420790796
 -0.09337701931162447
  0.13981331993408974
 -0.35336394739428656
  0.15558744225289545
 -0.40662298430620564
 -0.35730868344039096
 -0.1643498428773268
  0.8129987319187173
  0.12381154008796186
  0.1643498546523043
  0.35730868667937543
(0.20836070527140493, 0.20836069727606882 + 0.0im)

so the triangle dual that get from Hypatia is the same constants as those that JuMP reports, just reshaped.

[araujoms]

I believe the issue is that what Hypatia sees is just a vector constraint; as an inner product between two vectors this dual variable is correct. JuMP, on the other hand, interprets this as an inner product between two matrices, which consists of not only the upper triangular that goes to Hypatia but also the redundant lower triangular. Well this redundant part makes the off-diagonal terms appear twice.

Not that if I use .== instead of ==, so that the redundant constraints are used, the dual matrix returned by JuMP is correct (modulo a nonzero imaginary part which is probably an independent bug).

[blegat]

Hypatia needs the MOI.Scaled{MOI.HermitianPositiveSemidefiniteConeTriangle} but its MOI wrapper says that it needs the MOI.HermitianPositiveSemidefiniteConeTriangle and hardcode the scaling in the MOI wrapper. Since it's easy to get these things wrong (as exemplified by this issue), we prefer doing everything in bridges that are shared with all the solvers rather than hardcoding this in the MOI wrapper (even if it's as simple as a rescaling). I don't know if the issue is there but updating Hypatia to use the MOI.Scaled cone might fix it.
In any case, we should build a minimal reproducible example where Hypatia return something else than other solvers and add it in the MOI conic tests because this issue shows a lack of coverage :)

[araujoms]

The issue is not restricted to Hypatia. I've also tried with Clarabel, SCS, and Mosek, and all have exactly the same problem, for both the real and complex PSD cones.

I've also tried editing Hypatia's SupportedCones to say MOI.Scaled{MOI.HermitianPositiveSemidefiniteConeTriangle} instead of MOI.HermitianPositiveSemidefiniteConeTriangle, it didn't make any difference.

[odow]

I believe the issue is that what Hypatia sees is just a vector constraint; as an inner product between two vectors this dual variable is correct. JuMP, on the other hand, interprets this as an inner product between two matrices, which consists of not only the upper triangular that goes to Hypatia but also the redundant lower triangular. Well this redundant part makes the off-diagonal terms appear twice.

Yes, I think this is correct.

[odow]

Here's a better example of what is happening:

julia> using JuMP

julia> using Clarabel

julia> using Hypatia

julia> using SCS

julia> using LinearAlgebra

julia> function hermitian_equality(optimizer)
           model = Model(optimizer)
           set_silent(model)
           @variable(model, x[1:2, 1:2])
           @objective(model, Min, sum(x))
           y = [x[1,1] (x[1,2]+x[2,1]im); (x[1,2]-x[2,1]im) x[2,2]]
           z = [1 2+3im; 2-3im 4]
           @constraint(model, c, Hermitian(y - z) == 0)
           optimize!(model)
           @assert is_solved_and_feasible(model)
           return dual(c)
       end
hermitian_equality (generic function with 1 method)

julia> function hermitian_broadcast(optimizer)
           model = Model(optimizer)
           set_silent(model)
           @variable(model, x[1:2, 1:2])
           @objective(model, Min, sum(x))
           y = [x[1,1] (x[1,2]+x[2,1]im); (x[1,2]-x[2,1]im) x[2,2]]
           z = [1 2+3im; 2-3im 4]
           @constraint(model, c, y .== z)
           optimize!(model)
           @assert is_solved_and_feasible(model)
           return dual.(c)
       end
hermitian_broadcast (generic function with 1 method)

julia> hermitian_equality(Clarabel.Optimizer)
2×2 Hermitian{ComplexF64, Matrix{ComplexF64}}:
 1.0+0.0im  1.0+1.0im
 1.0-1.0im  1.0+0.0im

julia> hermitian_broadcast(Clarabel.Optimizer)
2×2 Matrix{ComplexF64}:
 1.0+0.0im  0.5+0.5im
 0.5-0.5im  1.0+0.0im

julia> hermitian_equality(Hypatia.Optimizer)
2×2 Hermitian{ComplexF64, Matrix{ComplexF64}}:
 1.0+0.0im  1.0+1.0im
 1.0-1.0im  1.0+0.0im

julia> hermitian_broadcast(Hypatia.Optimizer)
2×2 Matrix{ComplexF64}:
 1.0+0.0im  0.0+0.0im
 1.0-1.0im  1.0+0.0im

julia> hermitian_equality(SCS.Optimizer)
2×2 Hermitian{ComplexF64, Matrix{ComplexF64}}:
 1.0+0.0im  1.0+1.0im
 1.0-1.0im  1.0+0.0im

julia> hermitian_broadcast(SCS.Optimizer)
2×2 Matrix{ComplexF64}:
 1.0+0.0im  0.5+0.5im
 0.5-0.5im  1.0+0.0im

[odow]

This also happens for Symmetric:

julia> function hermitian_equality(optimizer)
           model = Model(optimizer)
           set_silent(model)
           @variable(model, x[1:2, 1:2])
           @objective(model, Min, sum(x))
           y = [x[1,1] (x[1,2]+x[2,1]im); (x[1,2]+x[2,1]im) x[2,2]]
           z = [1 2+3im; 2+3im 4]
           @constraint(model, c, Symmetric(y - z) == 0)
           optimize!(model)
           @assert is_solved_and_feasible(model)
           return dual(c)
       end
hermitian_equality (generic function with 1 method)

julia> function hermitian_broadcast(optimizer)
           model = Model(optimizer)
           set_silent(model)
           @variable(model, x[1:2, 1:2])
           @objective(model, Min, sum(x))
           y = [x[1,1] (x[1,2]+x[2,1]im); (x[1,2]+x[2,1]im) x[2,2]]
           z = [1 2+3im; 2+3im 4]
           @constraint(model, c, y .== z)
           optimize!(model)
           @assert is_solved_and_feasible(model)
           return dual.(c)
       end
hermitian_broadcast (generic function with 1 method)

julia> hermitian_equality(Clarabel.Optimizer)
2×2 Symmetric{ComplexF64, Matrix{ComplexF64}}:
 1.0+0.0im  1.0+1.0im
 1.0+1.0im  1.0+0.0im

julia> hermitian_broadcast(Clarabel.Optimizer)
2×2 Matrix{ComplexF64}:
 1.0+0.0im  0.5+0.5im
 0.5+0.5im  1.0+0.0im

julia> hermitian_equality(Hypatia.Optimizer)
2×2 Symmetric{ComplexF64, Matrix{ComplexF64}}:
 1.0+0.0im  1.0+1.0im
 1.0+1.0im  1.0+0.0im

julia> hermitian_broadcast(Hypatia.Optimizer)
2×2 Matrix{ComplexF64}:
 1.0+0.0im  0.0+0.0im
 1.0+1.0im  1.0+0.0im

julia> hermitian_equality(SCS.Optimizer)
2×2 Symmetric{ComplexF64, Matrix{ComplexF64}}:
 1.0+0.0im  1.0+1.0im
 1.0+1.0im  1.0+0.0im

julia> hermitian_broadcast(SCS.Optimizer)
2×2 Matrix{ComplexF64}:
 1.0+0.0im  0.5+0.5im
 0.5+0.5im  1.0+0.0im

So it is just the case that when we add these symmetric equality constraints we're doubling up on the duals.