# Solving a large, sparse, asymmetric linear complementarity problem

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.

• 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) – jlperla Oct 10 '17 at 20:26
• Correct. Note that I also added some notes about solving as a QCP. – Erwin Kalvelagen Oct 11 '17 at 1:16