How to handle a non-separable bilinear objective function in the special case of decoupled constraints?

Question

I have a large number of (10000+) non-negative, real decision variables $x_i$ and $y_j$.

Let $I$ and $J$ be the index sets associated with $x$ and $y$, respectively.

Let $\bar{I}$ and $\bar{J}$ be non-empty subsets of $I$ and $J$, respectively.

The objective function I would like to minimize (subject to some constraints) is of the form \begin{aligned} \quad & \sum_{i\in{\bar{I}}}\sum_{j\in{\bar{J}}} a_{ij}x_iy_j + \sum_{i\in{I}} b_ix_i + \sum_{j\in{J}} c_jy_j \\ \end{aligned} where $a_{ij}$, $b_i$, $c_j$ are positive, real constants.

All constraints are linear and decoupled/separable in the sense that each constraint involves either only $x_i$ or only $y_j$. For example,

\begin{aligned} x_1 \leq x_2 \\ y_1 \leq y_2 \end{aligned}

are acceptable constraints, but

\begin{aligned} x_1 \leq y_1 \\ \end{aligned}

is not an acceptable constraint.

Number of constraints are also large (comparable to the number of variables).

In summary, objective function is not separable (to an $x$-part and a $y$-part), but the constraints are. There is one paper¹ from 2009 which defines such problems as "separable bilinear programs" but it does not seem to be a commonly used term.

What is a good way of solving such problems, if there is any? I have a concern that a solver may use a general method for a quadratic program and not fully make use of the separability of constraints.

Edit: Linking previously asked questions which seem related.

Linearization of the product of two real valued variables - Binary expansion approach

How to linearize the product of two continuous variables?

---

_References

_{[1] Petrik, M., & Zilberstein, S. (2009). A bilinear programming approach for multiagent planning. Journal of Artificial Intelligence Research, 35, 235-274. https://arxiv.org/pdf/1401.3461.pdf}

Answer 1

First, there are some solvers (including CPLEX) that will let you try to solve a model with linear constraints and a nonconvex quadratic objective function. In a case like that CPLEX will let you choose a "target", either a possibly local solution satisfying first order constraints or a provably optimal solution. The latter requires changing the problem to a mixed integer quadratic problem, possibly including using the McCormick envelopes in the page @wordsmithhelper linked.

I ran some tests, on problems much smaller than yours. Using the global optimality target, CPLEX found what I suspect is an optimal solution immediately and then spent considerable time trying to prove optimality. On one test problem (50 $x$ variables with 119 constraints, 75 $y$ variables with 135 constraints), CPLEX did in fact prove that the early solution was optimal. So, if your solver supports this, a possible heuristic is to tell it to solve the QP and set a tight time limit (or, if you want proven optimality, set a not-so-tight time limit and cross your fingers).

If you don't want to mess with QPs (or your solver won't let you), here is another heuristic that consistently matched the initial CPLEX solution in comparable time to what it took CPLEX to find that solution. Let $X$ and $Y$ be the feasible regions (polytopes) for $x$ and $y$. Note that for fixed (feasible) $x$, the problem becomes an LP in $y$ with objective coefficient $\sum_i a_{ij}x_i + c_j$ for $y_j$, and similarly for fixed (feasible) $y$. Let $f(x;y)$ be the linear objective function when $y$ is fixed and $g(y;x)$ the linear objective function when $x$ is fixed. The heuristic is simple.

1. Minimize $f(x;0)$ over $X$ to get $x^0.$ (For this step, it does not matter whether $0\in Y.$)
2. Minimize $g(y;x^0)$ over $Y$ to get $y^0.$ You now have an incumbent solution $(x^0,y^0)$ with cost $b^\prime x^0 + g(y^0;x^0).$ Set $t = 0.$
3. Minimize $f(x;y^t)$ over $X$ to get $x^{t+1}.$ If $f(x^{t+1};y^t) + c^\prime y^t$ is less than the current incumbent value, record $(x^{t+1}, y^t)$ as a new incumbent and keep going, else quit.
4. Minimize $g(y;x^{t+1})$ over $Y$ to get $y^{t+1}$ and compare $g(y^{t+1};x^{t+1}) + b^\prime x^{t+1}.$ If the new solution is better, record $(x^{t+1}, y^{t+1})$ as your new incumbent, increment $t$ and go back to step 3.

So you solve LPs, alternating between $x$ and $y$, until improvement stops. Note that each solution to either LP is feasible the next time you come back to that LP, since only the objective function is changing. Assuming your LP solver will "hot start" from a previous solution, this should speed up the process.