Polynomial optimization without sum of squares

Specific Domains Optimization (Mathematical)
,
1

TL-DR: I would like to solve a linear program, but the problem itself involves multivariate polynomials. So I need tools to transform polynomials to their vectors of coefficients and back. Moreoever, in the end, after solving the corresponding linear program, I would like to reinterpret the solution as a ceritificate of positivity of the polynomial I am interested in, over the region of interest (defined using linear inequalities).

While I can probably write it from scratch, yet most of the functionality I need is probably already being used in SumOfSquares or Polynomials or in some other Julia package. Some comments would be appreciated.

Suppose you have 3 variables, x, y, z, say, and linear inequalities of the form

0 <= x <= 1
0 <= y <= 1
0 <= z <= 1
1 <= x + y + z <= 2

This gives us a number of nonnegative polynomials,

x
y
z
1-x
1-y
1-z
x+y+z-1
2-x-y-z

I would like to form the vector of all polynomials of degree up to 2 in the previous polynomials. This should for example contain 1, x, y, z, 1-x, etc., but also x^2, xy, x(1-x) and so on. Let us call this vector X, say.

I am given a polynomial p

p = x*(1-y) + y*(1-z) + z*(1-x)

I am interested in the problem of writing, if possible, p as a linear combination of the polynomials in X with nonnegative coefficients and, moreover, finding such an expression with the highest possible coefficient of 1.

Note that I specifically avoided mentioning sums of squares, because in the problem I am interested in, each variable actually represents a matrix, the product is the symmetric Kronecker product and the inequality represents the Loewner order. As a result, squares of polynomials in x, y and z are not necessarily nonnegative.

That being said, I wonder whether the above problem can be solved in Julia, for example using Jump and SumOfSquares or some other packages. I looked for a while in the documentation, but could not yet find what I want. If you know how I can solve the problem above, could you please help, or maybe give some clues? Thank you.

2

Hi @joemaik, welcome to the forum :smile:

I don’t know the answer, so this is a question for @blegat.

3

Hi @odow, I wrote the code from scratch, by basically defining the cone of nonnegative polynomials and converting my problem to a linear program which I then solved, using JuMP and an LP solver (HiGHS).

2 Likes
4

Would you mind showing us the code?

1 Like
5 
using DynamicPolynomials
using JuMP
using HiGHS

@polyvar x[1:3]

P = x[1] * (1 - x[2]) + x[2] * (1 - x[3]) + x[3] * (1 - x[1])

S = [x[1], 1-x[1], x[2], 1-x[2], x[3], 1-x[3], x[1]+x[2]+x[3]-1, 2-x[1]-x[2]-x[3]]

poly_basis = monomials([x[1], x[2], x[3]], 0:2)

function coeffs(poly, poly_basis)
    return [DynamicPolynomials.coefficient(poly, monom) for monom in poly_basis]
end

function build_poly(coeffs, poly_basis)
    return sum(coeffs .* poly_basis)
end

build_poly(coeffs(P, poly_basis), poly_basis) == P

@polyvar p[1:length(S)]

constraints_short = monomials([p[i] for i in 1:length(S)], 0:2)

my_subs = tuple((p[i] for i in 1:length(S))...) => tuple((S[i] for i in 1:length(S))...)

constraints = [subs(c, my_subs) for c in constraints_short]

A = zeros((length(poly_basis), length(constraints)));

size(A)

for j in 1:length(constraints)
    A[:, j] = coeffs(constraints[j], poly_basis)
end

rhs = coeffs(P, poly_basis)

model = Model(HiGHS.Optimizer)
@variable(model, c[1:length(constraints)] .>= 0)
@objective(model, Max, c[1])
@constraint(model, A*c == rhs)

print(model)

optimize!(model)

c_opt = [value(c[i]) for i in 1:length(c)];

certificate = sum(c_opt .* constraints_short)
1 Like
6

The substitution part is slow though (not in the example above, but in some higher degree examples). I wonder if I can optimize that part of the code…

7

What you are doing (multiplying x with 1 - y to get a new nonnegative polynomial) is called the Schmüdgen certificate, see Certificate · SumOfSquares.
The quadratic case where you take pairwise products of linear functions is a special case of Schmüdgen.
In general, you can multiply by SOS polynomials but in this case, because your products are already of degree 2, the SOS polynomial should have degree 0 which mean a nonnegative constant, this is your c.
Not quite sure why the substitution is slow but you could do a code a bit more specialized to speed it up like iterating over all pairs.
Using the Schmüdgen example I mentioned above with a maxdegree = 2 might be faster, you can try.

3 Likes
8

Thank you @blegat! Yes, I have actually read the Schmüdgen certificate page and the code you have mentioned. You are right in that in this specific example, there is no need to “turn off” the SOS polynomials, but I am also interested in higher degree examples, for which it is important for me not to use SOS polynomials. And I don’t understand enough the code related to the Schmüdgen certificate to modify it so that it only forms the products of the nonnegative polynomials (on the region) I have provided and no SOS polynomials (in my ring, squares of arbitrary polynomials are not necessarily nonnegative). Is it please possible to modify the code for the Schmüdgen certificate so that it does not form SOS polynomials but only products of the nonnegative polynomials I have provided?

1 Like
9

Yes, I suppose this should do the trick

function SOSC.multiplier_basis(certificate::Schmüdgen, index::SOSC.PreorderIndex, domain::SOSC.WithVariables)
    return SOSC.maxdegree_gram_basis(certificate.basis, variables(domain), 0)
end
1 Like
10

Thank you @blegat. I tried it, but unfortunately, it is giving a lower bound of 2/3, which is too big (the actual best lower bound should be 1/2, if one does not use SOS polynomials), which makes me suspect that it is still using SOS polynomials. Here is the adapted code:

using DynamicPolynomials
using Combinatorics
using SemialgebraicSets
import SCS
using Dualization

solver = dual_optimizer(SCS.Optimizer)

n = 3

@polyvar x[1:n, 1:n]

P = 0
for i in 1:n
    myprod = 1
    for j in [j for j in 1:n if j != i]
        myprod *= x[i, j]
    end
    P += myprod
end

ineqs1 = [x[i, j] for i in 1:n for j in 1:n if i != j]
ineqs2 = [[x[i, j] + x[j, k] + x[k, i] - 1, x[j, i] + x[k, j] + x[i, k] - 1] 
    for (i, j, k) in combinations(1:n, 3)]
ineqs2 = vcat(ineqs2...)
ineqs = vcat(ineqs1..., ineqs2...)
eqs = [x[i, j] + x[j, i] - 1 for i in 1:n for j in 1:n if i < j]

S = basic_semialgebraic_set(algebraic_set(eqs), ineqs)

using SumOfSquares

import MultivariateBases as MB
const SOS = SumOfSquares
const SOSC = SOS.Certificate
struct Schmüdgen{IC <: SOSC.AbstractIdealCertificate, CT <: SOS.SOSLikeCone, BT <: SOS.AbstractPolynomialBasis} <: SOSC.AbstractPreorderCertificate
    ideal_certificate::IC
    cone::CT
    basis::Type{BT}
    maxdegree::Int
end

SOSC.cone(certificate::Schmüdgen) = certificate.cone

function SOSC.preprocessed_domain(::Schmüdgen, domain::BasicSemialgebraicSet, p)
    return SOSC.with_variables(domain, p)
end

function SOSC.preorder_indices(::Schmüdgen, domain::SOSC.WithVariables)
    n = length(domain.inner.p)
    if n >= Sys.WORD_SIZE
        error("There are $(2^n - 1) products in Schmüdgen's certificate, they cannot even be indexed with `$Int`.")
    end
    return map(SOSC.PreorderIndex, 1:(2^n-1))
end

function SOSC.multiplier_basis(certificate::Schmüdgen, index::SOSC.PreorderIndex, domain::SOSC.WithVariables)
    return SOSC.maxdegree_gram_basis(certificate.basis, variables(domain), 0)
end

function SOSC.multiplier_basis_type(::Type{Schmüdgen{IC, CT, BT}}) where {IC, CT, BT}
    return BT
end

function SOSC.generator(::Schmüdgen, index::SOSC.PreorderIndex, domain::SOSC.WithVariables)
    I = [i for i in eachindex(domain.inner.p) if !iszero(index.value & (1 << (i - 1)))]
    return prod([domain.inner.p[i] for i in eachindex(domain.inner.p) if !iszero(index.value & (1 << (i - 1)))])
end

SOSC.ideal_certificate(certificate::Schmüdgen) = certificate.ideal_certificate
SOSC.ideal_certificate(::Type{<:Schmüdgen{IC}}) where {IC} = IC

SOS.matrix_cone_type(::Type{<:Schmüdgen{IC, CT}}) where {IC, CT} = SOS.matrix_cone_type(CT)

model = SOSModel(solver)
@variable(model, α)
@objective(model, Max, α)
ideal_certificate = SOSC.Newton(SOSCone(), MB.MonomialBasis, tuple())
certificate = Schmüdgen(ideal_certificate, SOSCone(), MB.MonomialBasis, maxdegree(P) + 1)
@constraint(model, c, P >= α, domain = S, certificate = certificate)
optimize!(model)
solution_summary(model)
11

@blegat, I’ve been trying to understand the code for the Schmüdgen certificate so that I can modify it, but the bitwise operations alone in multiplier_basis took me an hour to understand :). By the way, I also need to change this part, because in my products, it is ok to have repeated factors in the products. They are not needed in the original certificate, since you are including SOS polynomials, which would handle those.

12

Ah yes indeed, if you don’t include SOS you might want duplicated factors