# Trying to solve a SDP problem with quadratic constraints

Hi, I’m trying to solve a SDP with quadratic constraints and I get the following error

``````Unable to transform a quadratic constraint into a second-order cone constraint because the quadratic constraint is not strongly convex.
``````

and my code is

``````using JuMP
using SCS

function anticommutator(ps1, ps2)
pstr = []
msigns = 0
m = 1
for (i, j) in zip(ps1, ps2)
if i == j
push!(pstr, 'I')
elseif i == 'I'
push!(pstr, j)
elseif j == 'I'
push!(pstr, i)
elseif i == 'X' && j == 'Y'
push!(pstr, 'Z')
msigns += 1
elseif i == 'X' && j == 'Z'
push!(pstr, 'Y')
msigns += 1
m += 1
elseif i == 'Y' && j == 'Z'
push!(pstr, 'X')
msigns += 1
elseif i == 'Y' && j == 'X'
push!(pstr, 'Z')
msigns += 1
m += 1
elseif i == 'Z' && j == 'X'
push!(pstr, 'Y')
msigns += 1
elseif i == 'Z' && j == 'Y'
push!(pstr, 'X')
msigns += 1
m += 1
end
end
if msigns % 2 == 0
m % 2, reduce(*, pstr)
else
(0, "")
end
end

function matrix_element(idx1, idx2, xm)
pstr1, s1 = idx1
pstr2, s2 = idx2
k = length(pstr1)
if s1 > s2
s1, s2 = s2, s1
pstr1, pstr2, = pstr2, pstr1
end
dist = s2 - s1
if k < dist
return -(xm[idx2] * xm[idx1])
end
sgn, str = anticommutator(relstr1, relstr2)
if str == "" || all(split(str, "") .== "I")
return -(xm[idx1] * xm[idx2])
end
if str[1] == 'I'
end
sgn = (-1) ^ sgn
sgn * xm[(str, s2)] -(xm[idx1] * xm[idx2])
end

N = 3
K = 2
J = 0.3
h = 0.1
model = Model(SCS.Optimizer)
operators = ['X', 'Y', 'Z', 'I']
sites = [i for i=1:N]
indices = collect(
Iterators.product(map(x -> reduce(*, x), collect(Iterators.product([i == 1 ? operators[1:3] : operators for i=1:K]...))), sites)
)
indices = reshape(indices, :)
n_ops = N * ((length(operators) - 1) * length(operators) ^(K-1))
if length(indices) != n_ops
ArgumentError("The indices do not match the number of operators")
end
@variable(model, X[1:n_ops])
@variable(model, M[1:n_ops, 1:n_ops], PSD)
Xm = Containers.DenseAxisArray(X, indices)
Mm = Containers.DenseAxisArray(M, indices, indices)
@constraint(model, c[i in indices, j in indices], Mm[i, j] == matrix_element(i, j, Xm))
@objective(model, Min, h * sum(Xm[("XI", i)] for i=1:N) + J * sum(Xm[("ZZ", i)] for i=1:N-1))
optimize!(model)

``````

I saw here that it is not possible to express non-convex terms as convex.

Is there a transformation I can do so that I’d be able to optimize the model?

Thanks for any help!

@blegat is our conic guru. SCS does not support quadratic equality constraints or quadratic psd matrices.

1 Like

With `N = 1`, this matrix:

``````[matrix_element(i, j, Xm) for i in indices, j in indices
``````

is

`````` -X[1]²               -X[1]*X[2]           -X[1]*X[3]           -X[1]*X[4]           …  -X[1]*X[9] - X[5]    -X[1]*X[10] - X[10]  -X[1]*X[11]          -X[1]*X[12]
-X[2]*X[1]           -X[2]²               -X[2]*X[3]           -X[2]*X[4] + X[9]       -X[2]*X[9] + X[4]    -X[2]*X[10]          -X[2]*X[11] - X[10]  -X[2]*X[12]
-X[3]*X[1]           -X[3]*X[2]           -X[3]²               -X[3]*X[4] - X[8]       -X[3]*X[9]           -X[3]*X[10]          -X[3]*X[11]          -X[3]*X[12] - X[10]
-X[4]*X[1]           -X[4]*X[2] + X[9]    -X[4]*X[3] - X[8]    -X[4]²                  -X[4]*X[9] + X[2]    -X[4]*X[10] - X[11]  -X[4]*X[11]          -X[4]*X[12]
-X[5]*X[1] - X[9]    -X[5]*X[2]           -X[5]*X[3] + X[7]    -X[5]*X[4]              -X[5]*X[9] - X[1]    -X[5]*X[10]          -X[5]*X[11] - X[11]  -X[5]*X[12]
-X[6]*X[1] + X[8]    -X[6]*X[2] - X[7]    -X[6]*X[3]           -X[6]*X[4]           …  -X[6]*X[9]           -X[6]*X[10]          -X[6]*X[11]          -X[6]*X[12] - X[11]
-X[7]*X[1]           -X[7]*X[2] - X[6]    -X[7]*X[3] + X[5]    -X[7]*X[4]              -X[7]*X[9]           -X[7]*X[10] - X[12]  -X[7]*X[11]          -X[7]*X[12]
-X[8]*X[1] + X[6]    -X[8]*X[2]           -X[8]*X[3] - X[4]    -X[8]*X[4] - X[3]       -X[8]*X[9]           -X[8]*X[10]          -X[8]*X[11] - X[12]  -X[8]*X[12]
-X[9]*X[1] - X[5]    -X[9]*X[2] + X[4]    -X[9]*X[3]           -X[9]*X[4] + X[2]       -X[9]²               -X[9]*X[10]          -X[9]*X[11]          -X[9]*X[12] - X[12]
-X[10]*X[1] - X[10]  -X[10]*X[2]          -X[10]*X[3]          -X[10]*X[4] - X[11]     -X[10]*X[9]          -X[10]²              -X[10]*X[11]         -X[10]*X[12]
-X[11]*X[1]          -X[11]*X[2] - X[10]  -X[11]*X[3]          -X[11]*X[4]          …  -X[11]*X[9]          -X[11]*X[10]         -X[11]²              -X[11]*X[12]
-X[12]*X[1]          -X[12]*X[2]          -X[12]*X[3] - X[10]  -X[12]*X[4]             -X[12]*X[9] - X[12]  -X[12]*X[10]         -X[12]*X[11]         -X[12]²
``````

This looks like `A - X * X'` where `A` is some matrix depending linearly on `X`.
According to the Schur complement, `A - X * X'` is PSD iff

``````[1 X'
X A]
``````

is PSD.
So try something like

``````@variable(model, X[1:n_ops])
Xm = Containers.DenseAxisArray(X, indices)
A = [matrix_element(i, j, Xm).aff for i in indices, j in indices]
@constraint(model, [1 X'; X A] in PSDCone())
@objective(model, Min, h * sum(Xm[("XI", i)] for i=1:N) + J * sum(Xm[("ZZ", i)] for i=1:N-1))
optimize!(model)
``````
2 Likes

Brilliant!

Thanks for the link to the Schur complement, it really helped clarify this.

2 Likes