# Dual to PSD Constraint in JuMP

**URL:** <https://discourse.julialang.org/t/dual-to-psd-constraint-in-jump/126352>\
**Category:** Optimization (Mathematical)\
**Tags:** jump\
**Created:** [February 26, 2025, 6:55pm UTC](https://discourse.julialang.org/t/dual-to-psd-constraint-in-jump/126352 "2025-02-26T18:55:58Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![lopenguin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lopenguin/32/215482_2.png) [@lopenguin](https://discourse.julialang.org/u/lopenguin)\
**Post date:** [February 26, 2025, 6:55pm UTC](https://discourse.julialang.org/t/dual-to-psd-constraint-in-jump/126352/1 "2025-02-26T18:55:59Z")

</div>

Hi all,

This is my first forum post so feedback is welcome. In summary, finding the dual to a PSD constraint in JuMP is not particularly intuitive and certainly not well-documented. I show the approach that worked for me (which is clearly suboptimal) and solicit suggests for a more intuitive procedure. I am happy to do implementation work to make this easier.

## Finding the dual

First, I show how I find the dual for the K-means clustering example (see JuMP simple semidefinite examples; the post won’t let me include a link). Here is the optimization problem:

### Optimization Problem (K-Means Clustering)

```julia
import LinearAlgebra, SCS
using JuMP
a = [[2.0, 2.0], [2.5, 2.1], [7.0, 7.0], [2.2, 2.3], [6.8, 7.0], [7.2, 7.5]]
m = length(a)
num_clusters = 2
W = zeros(m, m)
for i in 1:m, j in i+1:m
    W[i, j] = W[j, i] = exp(-LinearAlgebra.norm(a[i] - a[j]) / 1.0)
end
model = Model(SCS.Optimizer)
set_silent(model)
@variable(model, Z[1:m, 1:m] >= 0, PSD)
@objective(model, Min, LinearAlgebra.tr(W * (LinearAlgebra.I - Z)))
@constraint(model, [i = 1:m], sum(Z[i, :]) .== 1)
@constraint(model, LinearAlgebra.tr(Z) == num_clusters)
optimize!(model)
assert_is_solved_and_feasible(model)
Z_val = value.(Z)

```

### Finding the dual to `Z`

```julia
psd_constraints = all_constraints(model, Vector{VariableRef}, MOI.PositiveSemidefiniteConeTriangle)
X = dual(psd_constraints[1])
n = Int(sqrt(2 * length(X) + 0.25) - 0.5)
y = zeros(n,n); y[LinearAlgebra.triu(trues(n,n))] = X
X = LinearAlgebra.triu(y,1)' + y

```

We can check this is the right dual variable by verifying the duality gap `tr(X'*Z_val)==0`.

To break this down: first, I query for all PSD constraints in `model`. I found the type to ask for by:

```julia
list_of_constraint_types(model)

```

Then, I take the dual of the first constraint (there is only 1 in this optimization problem). The `PositiveSemidefiniteConeTriangle` constraint is a vectorized triangular representation, so I convert it into a symmetric matrix.

## Ideas for Improvement

As a user, the most intuitive way to get the dual of a variable constraint would be to call the `dual` function on the semidefinitive variable: `dual(Z)`. I recognize this is somewhat inconsistent with the existing workflow for referencing constraints set upon variable creation:

- `BinaryRef` / `IntegerRef`
- `FixRef` / `UpperBoundRef` / `LowerBoundRef`

Perhaps we need a `PSDRef` which operates on a matrix/vector variable. It’s not clear to me whether that would be able to inherit from `GenericVariableRef`. There’s also probably a little more nuance related to using `PositiveSemidefiniteConeTriangle` vs. `PositiveSemidefiniteConeSquare` constraints from MOI.

## Conclusion

As primarily a user of JuMP, I am very interested to hear what the developer community thinks about adding a function to make finding the dual to the PSD constraint easier. I am happy to do implementation work. At the very least, this post should serve as documentation of an existing way to find duals to PSD constraints in JuMP.

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [February 26, 2025, 7:42pm UTC](https://discourse.julialang.org/t/dual-to-psd-constraint-in-jump/126352/2 "2025-02-26T19:42:46Z")

</div>

Hi @lopenguin, welcome to the forum 😄

You’re correct that we don’t have a good way to get the constraint associated with a variable set.

I’ve opened an issue [Constraint associated with a variable set · Issue #3951 · jump-dev/JuMP.jl · GitHub](https://github.com/jump-dev/JuMP.jl/issues/3951)

A work-around for now is to construct the constraint separately:

```Julia
julia> import LinearAlgebra, SCS

julia> using JuMP

julia> begin
           a = [[2.0, 2.0], [2.5, 2.1], [7.0, 7.0], [2.2, 2.3], [6.8, 7.0], [7.2, 7.5]]
           m = length(a)
           num_clusters = 2
           W = zeros(m, m)
           for i in 1:m, j in i+1:m
               W[i, j] = W[j, i] = exp(-LinearAlgebra.norm(a[i] - a[j]) / 1.0)
           end
           model = Model(SCS.Optimizer)
           set_silent(model)
           @variable(model, Z[1:m, 1:m] >= 0, Symmetric)
           @constraint(model, Z_ref, Z in PSDCone())
           @objective(model, Min, LinearAlgebra.tr(W * (LinearAlgebra.I - Z)))
           @constraint(model, [i = 1:m], sum(Z[i, :]) .== 1)
           @constraint(model, LinearAlgebra.tr(Z) == num_clusters)
           optimize!(model)
           assert_is_solved_and_feasible(model)
           Z_val = value.(Z)
           Z_dual = dual(Z_ref)
       end
6×6 LinearAlgebra.Symmetric{Float64, Matrix{Float64}}:
  0.436055 -0.185689 0.000335814 -0.250273 0.000107101 -0.000536029
 -0.185689 0.436056 0.000324706 -0.250272 6.26428e-5 -0.000482078
  0.000335814 0.000324706 0.520101 -0.000435283 -0.338518 -0.181809
 -0.250273 -0.250272 -0.000435283 0.500357 -0.000335969 0.000958925
  0.000107101 6.26428e-5 -0.338518 -0.000335969 0.482705 -0.144021
 -0.000536029 -0.000482078 -0.181809 0.000958925 -0.144021 0.325889

julia> Z_dual' * Z_val
6×6 Matrix{Float64}:
  1.13355e-5 9.63657e-6 -1.00692e-5 9.24534e-6 -1.00654e-5 -1.00675e-5
  9.81162e-6 1.15108e-5 -1.02443e-5 9.42015e-6 -1.02406e-5 -1.02421e-5
  2.43498e-5 2.43496e-5 -2.25635e-5 2.43565e-5 -2.57207e-5 -2.48098e-5
 -2.11285e-5 -2.11288e-5 2.0296e-5 -1.86478e-5 2.0287e-5 2.02912e-5
 -1.79189e-5 -1.79191e-5 1.65731e-5 -1.79273e-5 1.95066e-5 1.77142e-5
 -6.44951e-6 -6.4492e-6 6.00798e-6 -6.44691e-6 6.23308e-6 7.11415e-6

```