Reformulations of a very slow problem involving absolute values

Ongoing

Carlos de Gois

• June 18, 2024 21:29

Hi everyone. I came across an interesting MIP instance that seems particularly difficult to solve. The problem is essentially the maximization of a sum of absolute values under some linear constraints. Here is a minimal working example:

from itertools import repeat
from functools import reduce
import numpy as np
import gurobipy as gp
from gurobipy import GRB
import time

A0 = np.array([[1,  1,  1,  1],
               [1,  1, -1, -1],
               [1, -1,  1, -1],
               [1, -1, -1,  1]])

B0 = np.array([[1,  0,  0,  0],
               [0,  1,  0,  0],
               [0,  0,  1,  0],
               [0,  0,  0, -1]])

# k integer controls the size of the instance
def mwe(k):
    ncoeffs = 4 ** k
    A = reduce(np.kron, repeat(A0, k))
    B = reduce(np.kron, repeat(B0, k))

    model = gp.Model()
    x = model.addMVar((ncoeffs,), lb=-1, ub=1, vtype=GRB.CONTINUOUS, name="x")

    y = A @ x
    model.addConstrs((y[i] >= 0 for i in range(ncoeffs)))
    model.addConstr(y.sum() == 1)
    tmp = A @ B @ x
    model.addConstrs((tmp[i] >= 0 for i in range(ncoeffs)))

    aux = model.addMVar(ncoeffs, vtype=GRB.CONTINUOUS, name="aux")
    model.addConstrs((aux[i] == gp.abs_(x[i]) for i in range(ncoeffs)))

    model.setObjective(aux.sum(), GRB.MAXIMIZE)
    model.optimize()
    return model

This runs almost immediately for k = 1 and k = 2, but already slows downs quite a lot for k = 3, even though this instance size does not seem to be large in absolute number of variables:

Presolved: 380 rows, 315 columns, 8694 nonzeros
Variable types: 252 continuous, 63 integer (63 binary)

Can you guide me through how to investigate why this model is so difficult to solve, and/or provide me hints on reformulations that might help?

Some things I have tried:

• To model the absolute values, I have tried using gurobipy's abs_() function as well as manually modelling it using auxiliary binary variables, which did not help. (By the way, is abs_() in general "smarter" than just doing the canonical auxiliary variables to model absolute values?)
• Tried to avoid absolute values altogether by modeling a quadratic optimization function where each term is (2 * z_i - 1) * x_i, with z_i binary and x_i in [0, 1]. Also did not seem to help.
• Used the automated parameter tuning tool. Specifically, did it for the k = 2 instances then assumed the same parameters would help the larger instances (is this usually true?). They seem to have some effect but not large enough to solve larger instances.

Thanks a lot for any help!

Related to

• General Constraints

0

• 

  Jaromił Najman

  Gurobi Staff

  • June 24, 2024 15:31

  Hi Carlos,

  To model the absolute values, I have tried using gurobipy's abs_() function as well as manually modelling it using auxiliary binary variables, which did not help. (By the way, is abs_() in general "smarter" than just doing the canonical auxiliary variables to model absolute values?)

  You did it the correct way. Using gurobipy's abs_() is usually smarter in the sense that Gurobi will internally compute (hopefully) tight big-M values for the binary formulation. If however, you can provide best possible big-M values for the binary formulation of the ABS function, then you can just write the formulation yourself.

  Tried to avoid absolute values altogether by modeling a quadratic optimization function where each term is (2 * z_i - 1) * x_i, with z_i binary and x_i in [0, 1]. Also did not seem to help.

  This is almost always a very bad idea when x is continuous, because it introduces more room for numerical inaccuracies without gaining any relaxation tightness.

  Used the automated parameter tuning tool. Specifically, did it for the k = 2 instances then assumed the same parameters would help the larger instances (is this usually true?). They seem to have some effect but not large enough to solve larger instances.

  Usually, parameters found by tuning are also useful for larger instances as long as the model structure stays the same, so what you see is correct behavior.

  I had a look at the model and the solutions found for k=1,2, and 3. It looks like you fix one variable to a given value and then in the optimal solution all other x attain the same value in absolute terms, i.e., \(|x_i| = |x_j|\) except for the one you fixed. Is this some property that may be true? 

  Altogether, I think that you should try to deduce some more properties of the model, which would help the solver, by making the formulation more precise. Our Tech-Talk on Weak vs Strong MIP formulations might be helpful here.

  Best regards, 
  Jaromił

  0
• 

  Carlos de Gois

  • July 10, 2024 13:46

  

  

  

  Hi Jaromił,

  Thanks a lot for the help and the tech-talk recommendation! In case you or anyone else can help me further, I investigated the model a bit further and for $k \geq 2$ something similar to your suggestion seems to be true:

  1. The value of $x_1$ can be predetermined.
  2. It seems to be true that $\lvert x_i \rvert = \lvert x_j \rvert$ for $i = 2, \ldots, 4^k$

  The straightforward way of adding this to the model is to have some constraints `model.addConstrs((aux[i] == aux[i + 1] for i in range(1, ncoeffs - 1)))` (see function `mwe_equals` below).

  • Is there a smarter way to do this?

  There is also another approach: Since the coefficients are all equal, I can try to guess the optimal value of the function and then run a binary feasibility program where I only try to find an attributed of $\pm 1$ signs that satisfies all the constraints (an attempt is in function `mwe_signs`).

  • Is this a better idea than before, or are feasibility problems not recommended?

  For both cases, I can solve $k = 2, 3$ much faster, but $k=4$ still seems difficult (as in 15min were not enough to get a solution so far).

  from itertools import repeat
  from functools import reduce
  import numpy as np
  import gurobipy as gp
  from gurobipy import GRB
  
  # k > 1
  def mwe_equals(k):
      ncoeffs = 4 ** k
      first = 1 / ncoeffs
      A = reduce(np.kron, repeat(A0, k))
      B = reduce(np.kron, repeat(B0, k))
  
      model = gp.Model()
      x = model.addMVar((ncoeffs,), lb=-1, ub=1, vtype=GRB.CONTINUOUS, name="x")
      model.addConstr(x[0] == first)
  
      y = A @ x
      model.addConstrs((y[i] >= 0 for i in range(ncoeffs)))
      model.addConstr(y.sum() == 1)
      tmp = A @ B @ x
      model.addConstrs((tmp[i] >= 0 for i in range(ncoeffs)))
  
      aux = model.addMVar(ncoeffs, vtype=GRB.CONTINUOUS, name="aux")
      model.addConstrs((aux[i] == gp.abs_(x[i]) for i in range(ncoeffs)))
      model.addConstrs((aux[i] == aux[i + 1] for i in range(1, ncoeffs - 1)))
  
      model.setObjective(aux.sum(), GRB.MAXIMIZE)
      model.optimize()
      return model
  
  # k integer > 1
  # k = 2 -> try coeff = 0.0208333333
  # k = 3 -> try coeff = 0.003125
  # k = 4 -> I believe it should be feasible for coeff = 0.0005361519
  def mwe_signs(k, coeff):
      ncoeffs = 4 ** k
      first = 1 / ncoeffs
      A = reduce(np.kron, repeat(A0, k))
      B = reduce(np.kron, repeat(B0, k))
  
      model = gp.Model()
      signs = model.addMVar((ncoeffs - 1,), vtype=GRB.BINARY, name="signs")
      x = [first] + [(2 * s - 1) * coeff for s in signs]
  
      y = A @ x
      model.addConstrs((y[i] >= 0 for i in range(ncoeffs)))
      model.addConstr(y.sum() == 1)
      tmp = A @ B @ x
      model.addConstrs((tmp[i] >= 0 for i in range(ncoeffs)))
  
      model.setObjective(first + (ncoeffs - 1) * coeff, GRB.MAXIMIZE)
      model.optimize()
      return model

  Thanks in advance for any help!

  0

Please sign in to leave a comment.