LCP with sparse matrix

This problem has a very efficient (linear time) solution, though it requires a bit of discussion...

Zeroth: clarifying the problem / LCP

Per clarifications in the comments, @FooBar says the original problem is elementwise min; we need to find a z (or v) such that

either the left argument is zero and the right argument is non-negative or the left argument is non-negative and the right argument is zero

As @FooBar correctly points out, a valid reparameterization leads to an LCP. However, below I show that an easier and more efficient solution to the original problem can be achieved (by exploiting structure in the original problem) without needing the LCP. Why should this be easier? Well, notice that an LCP has a quadratic term in z (Bz+q)'z, but that the original problem doesn't (only linear terms Bz+q and z). I'll exploit that fact below.

First: simplify

There is an important but key detail that simplifies this problem in a major way:

• Create four matrices A11, A12, A21, A22 using scipy.sparse.diags
• And stack them together as A = scipy.sparse.bmat([[A11, A12], [A21, A22]])

This has huge implications. Specifically, this is not a single large problem, but rather a large number of really small (2D, to be precise) problems. Notice that the block diagonal structure of this A matrix is preserved throughout all subsequent operations. And every subsequent operation is a matrix-vector multiply or an inner product. This means that really this program is separable in pairs of z (or v) variables.

To be specific, suppose each block A11,... is size n by n. Then notice critically that z_1 and z_{n+1} appear only in equations and terms with themselves, and never elsewhere. So the problem is separable into n problems, each of which is 2 dimensional. Thus, I will hereafter solve the 2D problem, and you can serialize or parallelize over n as you see fit, without sparse matrices or big opt packages.

Second: the geometry of the 2D problem

Assume we have the elementwise problem in 2D, namely:

find z such that (elementwise) min( Bz + q , z ) = 0, or declare that no such `z` exists.

Because in our setup B is now 2x2, this problem geometry corresponds to four scalar inequalities that we can manually enumerate (I’ve named them a1,a2,z1,z2):

“a1”: B11*z1 + B12*z2 + q1 >=0“a2”: B21*z1 + B22*z2 + q2 >=0“z1”: z1 >= 0“z2:” z2 >= 0

This represents a (possibly empty) polyhedra, aka the intersection of four half spaces in 2d space.

Third: solving the 2D problem

(Edit: updated this bit for clarity)

What specifically is the 2D problem then? We want to find a z where one of the following solutions (which are not all distinct, but that won’t matter) is achieved:

1. a1>=0, z1=0, a2>=0, z2=0
2. a1=0, z1>=0, a2=0, z2>=0
3. a1>=0, z1=0, a2=0, z2>=0
4. a1=0, z1>=0, a2>=0, z2=0

If one of these is achieved, we have a solution: a z where the elementwise minimum of z and Bz+q is the 0 vector. If none of the four can be achieved, the problem is infeasible and we will declare that no such z exists.

The first case occurs when the intersection point of a1,a2 is in the positive orthant. Just choose the solution z = B^-1q, and check if it is elementwise nonnegative. If it is, return B^-1q (note that, even though B is not psd, I am assuming it is invertible/full rank). So:

if B^-1q is elementwise nonnegative, return z = B^-1q.

The second case (not entirely distinct from the first) occurs when the intersection point of a1,a2 is anywhere but does contain the origin. In otherwords, whenever Bz+q >0 for z = 0. This occurs when q is elementwise positive. So:

elif q is elementwise nonnegative, return z = 0.

The third case has z1=0, which can be substituted into a2 to show that a2=0 when z2 = -q2/B22. z2 must be >=0, so -q2/B22 >=0. a1 must also be >=0, which substituting in the values for z1 and z2, gives -B22/B12*q2 + q1 >=0. So:

elif -q2/B22 >=0 and  -B22/B12*q2 + q1 >=0, return z1= 0, z2 = -q2/B22.

The fourth step is the same as the third, but swapping 1 and 2. So:

elif -q1/B11 >=0 and -B21/B11*q1 + q2 >=0, return z1 = -q1/B11, z2 =0.

The final fifth case is when the problem is infeasible. That occurs when none of the above conditions are met. So:

else return None 

Finally: put the pieces together

Solving each 2D problem is a couple of simple/efficient/trivial linear solves and compares. Each will return a pair of numbers or None. Then just do same over all n 2D problems, and concatenate the vector. If any are None, the problem is infeasible (all None). Else, you have your answer.

LCP solver for Python based on AMPLPY

as @denfromufa pointed out, there's an AMPL interface to the PATH solver. He suggested Pyomo since it's open source and able to process AMPL. However, Pyomo turned out to be slow and tedious to work with. I have ultimately written my own interface to the PATH solver in cython and hope to release that at some point, but at this moment it has no input sanitation, is quick and dirty and I don't see the time of working on it.

For now, I can share an answer that uses the python extension of AMPL. It's not as fast as a direct interface to PATH: For every LCP we want to solve, it creates a (temporary) model file, runs AMPL, and collects the output. It's somewhat quick and dirty, but I felt that I should at least report parts of the results of the several past months since asking this question.

import os# PATH license string needs to be updatedos.environ['PATH_LICENSE_STRING'] = "3413119131&Courtesy&&&USR&54784&12_1_2016&1000&PATH&GEN&31_12_2017&0_0_0&5000&0_0"from amplpy import AMPL, Environment, dataframeimport numpy as npfrom scipy import sparsefrom tempfile import mkstempimport osimport sysimport contextlibclass DummyFile(object):    def write(self, x): pass@contextlib.contextmanagerdef nostdout():    save_stdout = sys.stdout    sys.stdout = DummyFile()    yield    sys.stdout = save_stdoutclass modFile:    # context manager: create temporary file, insert model data, and supply file name    # apparently, amplpy coders are inable to support StringIO    content = """        set Rn;        param B {Rn,Rn} default 0;        param q {Rn} default 0;        var x {j in Rn};        s.t. f {i in Rn}:                0 <= x[i]             complements                sum {j in Rn} B[i,j]*x[j]                 >= -q[i];    """    def __init__(self):        self.fd = None        self.temp_path = None    def __enter__(self):        fd, temp_path = mkstemp()        file = open(temp_path, 'r+')        file.write(self.content)        file.close()        self.fd = fd        self.temp_path = temp_path        return self    def __exit__(self, exc_type, exc_val, exc_tb):        os.close(self.fd)        os.remove(self.temp_path)def solveLCP(B, q, x=None, env=None, binaryDirectory=None, pathOptions={'logfile':'logpath.tmp' }, verbose=False):    # x: initial guess    if binaryDirectory is not None:        env = Environment(binaryDirectory='/home/foo/amplide.linux64/')    if verbose:        pathOptions['output'] = 'yes'    ampl = AMPL(environment=env)    # read model    with modFile() as mod:        ampl.read(mod.temp_path)    n = len(q)    # prepare dataframes    dfQ = dataframe.DataFrame('Rn', 'c')    for i in np.arange(n):        # shitty amplpy does not support np.float        dfQ.addRow(int(i)+1, np.float(q[i]))    dfB = dataframe.DataFrame(('RnRow', 'RnCol'), 'val')    if sparse.issparse(B):        if not isinstance(B, sparse.lil_matrix):            B = B.tolil()        dfB.setValues({            (i+1, j+1): B.data[i][jPos]            for i, row in enumerate(B.rows)            for jPos, j in enumerate(row)            })    else:        r = np.arange(n) + 1        Rrow, Rcol = np.meshgrid(r, r, indexing='ij')        #dfB.setColumn('RnRow', [np.float(x) for x in Rrow.reshape((-1), order='C')])        dfB.setColumn('RnRow', list(Rrow.reshape((-1), order='C').astype(float)))        dfB.setColumn('RnCol', list(Rcol.reshape((-1), order='C').astype(float)))        dfB.setColumn('val', list(B.reshape((-1), order='C').astype(float)))    # set values    ampl.getSet('Rn').setValues([int(x) for x in np.arange(n, dtype=int)+1])    if x is not None:        dfX = dataframe.DataFrame('Rn', 'x')        for i in np.arange(n):            # shitty amplpy does not support np.float            dfX.addRow(int(i)+1, np.float(x[i]))        ampl.getVariable('x').setValues(dfX)    ampl.getParameter('q').setValues(dfQ)    ampl.getParameter('B').setValues(dfB)    # solve    ampl.setOption('solver', 'pathampl')    pathOptions = ['{}={}'.format(key, val) for key, val in pathOptions.items()]    ampl.setOption('path_options', ' '.join(pathOptions))    if verbose:        ampl.solve()    else:        with nostdout():            ampl.solve()    if False:        bD = ampl.getParameter('B').getValues().toDict()        qD = ampl.getParameter('q').getValues().toDict()        xD = ampl.getVariable('x').getValues().toDict()        BB = ampl.getParameter('B').getValues().toPandas().values.reshape((n, n,), order='C')        qq = ampl.getParameter('q').getValues().toPandas().values[:, 0]        xx = ampl.getVariable('x').getValues().toPandas().values[:, 0]        ineq2 = BB.dot(xx) + qq        print((xx * ineq2).min(), (xx * ineq2).max() )    return ampl.getVariable('x').getValues().toPandas().values[:, 0]if __name__ == '__main__':    # solve problem from the Julia port at https://github.com/chkwon/PATHSolver.jl    n = 4    B = np.array([[0, 0, -1, -1], [0, 0, 1, -2], [1, -1, 2, -2], [1, 2, -2, 4]])    q = np.array([2, 2, -2, -6])    BSparse = sparse.lil_matrix(B)    env = Environment(binaryDirectory='/home/foo/amplide.linux64/')    print(solveLCP(B, q, env=env))    print(solveLCP(BSparse, q, env=env))    # to test licensing    from numpy import random    n = 1000    B = np.diag((random.randn(n)))    q = np.ones((n,))    print(solveLCP(B, q, env=env))