# Trying to solve a SDP problem with quadratic constraints

**URL:** <https://discourse.julialang.org/t/trying-to-solve-a-sdp-problem-with-quadratic-constraints/104461>\
**Category:** Optimization (Mathematical)\
**Created:** [October 1, 2023, 10:15am UTC](https://discourse.julialang.org/t/trying-to-solve-a-sdp-problem-with-quadratic-constraints/104461 "2023-10-01T10:15:11Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![golanor](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/golanor/32/108140_2.png) [@golanor](https://discourse.julialang.org/u/golanor)\
**Post date:** [October 1, 2023, 10:15am UTC](https://discourse.julialang.org/t/trying-to-solve-a-sdp-problem-with-quadratic-constraints/104461/1 "2023-10-01T10:15:11Z")

</div>

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

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

```

and my code is

```julia
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](https://discourse.julialang.org/t/unable-to-transform-a-quadratic-constraint-into-a-second-order-cone-constraint-because-the-quadratic-constraint-is-not-strongly-convex/79767) 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!

---

<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:** [October 2, 2023, 2:23pm UTC](https://discourse.julialang.org/t/trying-to-solve-a-sdp-problem-with-quadratic-constraints/104461/2 "2023-10-02T14:23:24Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![blegat](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/blegat/32/217090_2.png) [@blegat](https://discourse.julialang.org/u/blegat)\
**Post date:** [October 3, 2023, 9:31am UTC](https://discourse.julialang.org/t/trying-to-solve-a-sdp-problem-with-quadratic-constraints/104461/3 "2023-10-03T09:31:55Z")

</div>

With `N = 1`, this matrix:

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

```

is

```julia
 -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](https://en.wikipedia.org/wiki/Schur_complement), `A - X * X'` is PSD iff

```julia
[1 X'
 X A]

```

is PSD.  
So try something like

```julia
@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)

```

---

<div class="post-metadata">

**Author:** ![golanor](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/golanor/32/108140_2.png) [@golanor](https://discourse.julialang.org/u/golanor)\
**Post date:** [October 3, 2023, 5:06pm UTC](https://discourse.julialang.org/t/trying-to-solve-a-sdp-problem-with-quadratic-constraints/104461/4 "2023-10-03T17:06:43Z")

</div>

Brilliant!

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