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)
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
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:
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/IntegerRefFixRef/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.
