Optimization of a multivariate polynomial - affine in each variable

Question

I'm trying to maximize a polynomial which has ${n\choose 3}$ variables: For starters I'm planning to do this for $n=7$.

$\displaystyle\sum_{1\leq i<j<k<l\le n}x_{ijk}x_{ikl}+x_{ijl}x_{jkl}-x_{ijk}x_{ijl}-x_{ijl}x_{ikl}-x_{ikl}x_{jkl}-x_{jkl}x_{ijk}+x_{ijk}x_{ijl}x_{ikl}x_{jkl}$

over $\{-1,1\}^{n\choose3}$.

My arguments so far:

Note that this is a discrete set. But since the polynomial is linear in each variable (checked by taking second derivative w.r.t. each variable), we can maximize over the entire cube $[-1,1]^{n\choose3}$ and deduce by working our way from an interior point to a corner by replacing each co-ordinate with a 1 or a -1 depending on which one gives a bigger value. Thus optimizing over a cube will automatically give a maximum on a corner.

This is where I'm stuck. What else should I do before running an optimization problem (if that)? Has anyone used Gloptipoly to optimize polynomials?

Answer 1

If more than one variable is fractional, your rounding procedure no longer works since the first rounding affects the second one. Although indeed the value increases in each stage, the greedy approach does not necessarily find the optimal value. Otherwise you might as well start the rounding procedure from a random point.

You have a difficult function (from an optimization perspective), yet only have $2^{35} \approx 34 \cdot 10^9$ feasible points. My suggestion is to compute the objective value for each feasible point.

Since Python is a scripted language, $2^{35}$ could still be a challenge. If you do not want to use a compiled language, look into Numba and @jit. You may have to write out the full polynomial. The following code should solve your problem in a few hours:

import itertools
import datetime
from numba import jit

n=7;
subscripts = [(i,j,k) for i in range(n) for j in range(n) for k in range(n) if i<j and j<k]
summation_index = [(i,j,k,l) for i in range(n) for j in range(n) for k in range(n) for l in range(n) if i<j and j<k and k<l]

def phase1():
    print(print_polynomial(summation_index))

def phase2():
    max_value = -float('inf')
    max_x = None
    count = 0
    for x in itertools.product([-1, 1], repeat=35):
        count += 1
        value = evaluate_polynomial(x)
        if value>max_value:
            max_value = value
            max_x = x
        if count%1000000 == 0:
            progress= 100 * count / 2**35
            print("{} {:.1f}%".format(datetime.datetime.now(), progress))
    print(max_x,max_value)

def get_x_str(i,j,k):
    return 'x[' + str(subscripts.index((i,j,k))) + ']'

def print_polynomial(summation_index):
    poly = ""
    for idx in summation_index:
        (i,j,k,l) = idx
        xijk = get_x_str(i,j,k)
        xikl = get_x_str(i,k,l)
        xijl = get_x_str(i,j,l)
        xjkl = get_x_str(j,k,l)
        poly += "+ {}*({}-{}-{}+{}*{}*{}) + {}*({}-{}) - {}*{}".format(xijk, xikl, xijl, xjkl, xijl, xikl, xjkl, xijl, xjkl, xikl, xikl, xjkl)
    return poly[2:]

@jit
def evaluate_polynomial(x):
    return 0

# run phase 1, copy the output into evaluate_polynomial, then run phase 2
phase1();