Trying to solve a SDP problem with quadratic constraints

Specific Domains Optimization (Mathematical)
1

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
    relstr1 = rpad(pstr1[1+dist:end], k, 'I')
    relstr2 = lpad(pstr2[1:k-dist], k, 'I')
    sgn, str = anticommutator(relstr1, relstr2)
    if str == "" || all(split(str, "") .== "I")
        return -(xm[idx1] * xm[idx2])
    end
    if str[1] == 'I'
        str = rpad(str[2:end], k, '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!

2

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

1 Like
3

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
4

Brilliant!

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

2 Likes