I’m trying to determine the most efficient form of modelling a problem for a package (QuantumNPA.jl). The natural way of writing the problem is in the dual or geometric form (see here for nomenclature); however, this gives terrible performance with Mosek. I can also write the problem in primal or standard form, which Mosek solves nicely, but then *setting up* the problem takes a long time. I’m flummoxed. I have a MATLAB package to solve the same problem with YALMIP (Moment), and there it’s easy: I write the problem in dual form, and YALMIP gives it to Mosek in primal form. Both setting up the problem and solving it are as fast as it gets.

I’d like to achieve the same result with JuMP; ideally I’d write the problem in dual form, but if that doesn’t work it’s ok, I can also write it in primal form, as long as setting up the problem is fast. I’ve written a simple example to illustrate the problem (only conceptually, the example is too simple for benchmarking to be possible):

```
using SparseArrays
using JuMP
using LinearAlgebra
using MosekTools
function build_basis()
d = 5
F = [sparse([i,j],[j,i],[1.0,1.0],d,d) for j=2:d for i=1:j-1]
pushfirst!(F,I(d))
return F
end
function npa_dual()
F = build_basis()
d = size(F[1],1)
nvars = length(F)
bigF = spzeros(d^2, nvars)
for i=1:nvars
bigF[:,i] .= vec(F[i])
end
coeffs = zeros(nvars)
coeffs[[6,7,9,10]] .= [1,1,1,-1]
model = Model(Mosek.Optimizer)
@variable(model, x[1:nvars])
@constraint(model, x[1] == 1)
moment_matrix = Symmetric(reshape(bigF*x,d,d))
@constraint(model, moment_matrix in PSDCone())
@objective(model, Max, dot(coeffs,x))
optimize!(model)
return objective_value(model)
end
function npa_primal()
F = build_basis()
d = size(F[1],1)
nvars = length(F)
coeffs = zeros(nvars)
coeffs[[6,7,9,10]] .= [1,1,1,-1]
model = Model(Mosek.Optimizer)
@variable(model, Z[1:d,1:d] in PSDCone())
for i=2:nvars
@constraint(model, dot(F[i],Z) + coeffs[i] == 0)
end
@objective(model, Min, coeffs[1] + dot(F[1],Z))
optimize!(model)
return objective_value(model)
end
```

In the primal formulation I can see that pretty much the entire time is spent in the line

```
@constraint(model, dot(F[i],Z) + coeffs[i] == 0)
```

Is there a smarter way of doing this?