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.