Take a relatively standard linear complementarity problem for $z \in\mathbb{R}^I$. It is ultimately solving the following variational problem \begin{align} \min\{Bz + q, z\} &= 0\\ \end{align} Which corresponds to the following LCP \begin{align} z^T(B z + q) &= 0 \\ z &\geq 0\\ B z + q &\geq 0 \\ \end{align}

A few notes on this:

  • $B$ is positive-definite
  • $B$ is not symmetric. This is the key difference compared to many setups.
  • $B$ is very sparse. If the size of $B$ is $I \times I$, then there are $\approx 3 I$ non-zeros.
  • The form of $B$ is mostly the sort of banded almost-tridiagonal setup that comes out of PDE discretization, but there are a few peculiarities that make the bandwidth equal to $I$ in some formulations of the problem.
  • In practice, I want to get $I$ up to at least 2000-10000, but would love to go much larger for some other applications.

My question: in practice, what is the best way to solve this sort of problem numerically, with the goal of (1) using pre-existing libraries since I don't want to write or test my own algorithms; and (2) getting the size of the system as large as I can. (I am relatively unconstrained in terms of budget, so if there is software that can handle this with minimal work, I can pay for it.)

One approach people have suggested is to reformulate LCPs as QP problems (e.g. Relation between bounded quadratic and linear complementary problems) and then throw a large-scale QP solver such as Gurobi at it. The problems here is that $B$ is definetely not symmetric in my setup, and my understanding is that there is no general connection of QP to LCP for asymmetric matrices (every reformulation I tried just ended up with an inverse of $(I - B)$ in it, and hence useless).

Any thoughts on the best algorithm/software or reformulation of the problem?


Maybe as a MIP?

\begin{align} & z^Ty=0 \\ & y= Bz + q\\ & y, z \ge 0 \end{align}

The first condition can be written as:

\begin{align} & z_i \le \delta_i \cdot z^{UP}_i\\ & y_i \le (1-\delta_i) \cdot y^{UP}_i\\ & \delta_i \in \{0,1\} \end{align}

Sounds a bit crazy, but this actually performs quite well when applied to non-convex QPs. Of course, in this formulation, you will need bounds on $z$ and $y$. They can sometimes be obtained from the original problem. Otherwise one can solve some optimization problems first to find bounds or use SOS1 sets.

This model does not use that $B$ being pos.def. (it works with any $B$).

If $B$ is pos.def., I believe this is a conic programming problem:

\begin{align} &z^TBz \le -z^Tq\\ &z\ge 0\\ &Bz+q\ge 0 \end{align}

This convex QCP can be solved directly with solvers like Cplex and Gurobi.

  • $\begingroup$ thanks so much for the response. Very interesting approach, I am curious if anyone else has experience with MIPs for LCPs? Just to be clear here, the $z_i^{UP}$ and the $y_i^{UP}$ are just any bounds I can figure out for the variables, and then I use a MIP just looking for a feasible solution (i.e. no objective) $\endgroup$ – jlperla Oct 10 '17 at 20:26
  • $\begingroup$ Correct. Note that I also added some notes about solving as a QCP. $\endgroup$ – Erwin Kalvelagen Oct 11 '17 at 1:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.