稀疏矩阵LCP

Question

我用大写字母表示矩阵，用小写字母表示向量。

对于向量v，我需要求解以下线性不等式系统：

min(rv - (u + Av), v - s) = 0

其中0是零的向量。

其中r是标量，u和s是向量，A是矩阵。

定义z = v-s、B=rI - A、q=-u + Bs，我可以将前面的问题重写为线性互补问题，并希望使用LCP求解器，例如来自openopt的

LCP(M, z): min(Bz+q, z) = 0

或者，用矩阵表示法：

z'(Bz+q) = 0
z >= 0
Bz + q >= 0

问题是我的方程组很大。要创建A，我

• 使用A11、A12、A21、A22创建四个矩阵
• 并将它们堆在一起作为A = scipy.sparse.bmat([[A11, A12], [A21, A22]])
• (这也意味着A是不对称的，因此一些有效的转换到QP问题是行不通的)

openopt.LCP显然无法处理稀疏矩阵:当我运行它时，我的计算机崩溃了。通常，A.todense()会导致内存错误。同样，compecon-python不能解决稀疏矩阵的LCP问题。

有哪些替代的LCP实现适合解决这个问题？

我真的不认为一般的“解决LCP的工具”问题需要的是样本数据，但是无论如何，我们开始了。

from numpy.random import rand
from scipy import sparse

n = 3000
r = 0.03

A = sparse.diags([-rand(n)], [0])
s = rand(n,).reshape((-1, 1))
u = rand(n,).reshape((-1, 1))

B = sparse.eye(n)*r - A
q = -u + B.dot(s)

q.shape
Out[37]: (3000, 1)
B
Out[38]: 
<3000x3000 sparse matrix of type '<class 'numpy.float64'>'
    with 3000 stored elements in Compressed Sparse Row format>

还有更多的提示：

• openopt.LCP与我的矩阵一起崩溃，我假设它在继续之前将矩阵转换为稠密。
• compecon-python完全失败了，但有一些错误，它显然需要密集的矩阵，没有稀疏性的退路。
• B不是半正定的，所以我不能把线性互补问题改为凸二次问题(QP)。
• 来自本博览会的所有QP稀疏解都要求问题是凸的，而这个问题不是
• 在朱莉娅，PATHSolver可以解决我的问题(获得许可)。但是，使用PyJulia (我这里的问题报告)从Python调用它存在问题。
• 此外，Matlab有一个LCP求解器，它显然可以处理稀疏矩阵，但是实现更古怪(我真的不想在Matlab上做这个)。

Answer 1

这个问题有一个非常有效的(线性时间)解决方案，尽管它需要一些讨论.

零:澄清问题/ LCP

根据注释中的说明，@FooBar表示，最初的问题是按元素划分的min；我们需要找到一个z (或v)，以便

左边的参数是零，右边的参数是非负的，或者左边的参数是非负的，右边的参数是零。

正如@FooBar正确指出的那样，有效的重参数化将导致LCP。然而，下面我展示了一个更简单、更有效的解决最初问题的方法(通过利用原来问题中的结构)，而不需要LCP。为什么这更容易呢？注意，LCP在z ( Bz+q )'z中有一个二次项，但是最初的问题并不是(只有线性项Bz+q和z)，下面我将利用这个事实。

优先:简化

有一个重要但关键的细节可以在很大程度上简化这个问题：

• 使用A11、A12、A21、A22创建四个矩阵
• 并将它们叠加为A= scipy.sparse.bmat([A11，A12，A21，A22])

这有着巨大的影响。具体来说，这不是一个单大问题，而是一个大数 of (确切地说，是2D)问题。注意，这个A矩阵的块对角线结构在所有后续操作中都保持不变。每一个后续的运算都是矩阵向量乘或内积。这意味着这个程序实际上是在z (或v)变量对中分离的。

具体来说，假设每个块A11,...都是按n计算的n大小。然后批判性地注意到，z_1和z_{n+1}只在方程和术语中出现，而不是在其他地方出现。因此，这个问题可以分解为n问题，每个问题都是二维的。因此，接下来我将解决2D问题，您可以在n上按您认为合适的方式序列化或并行化，而不需要稀疏矩阵或大的opt包。

第二:二维问题的几何

假设我们在2D中有一个元素问题，即：

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

因为在我们的设置中，B现在是2x2，这个问题几何对应于我们可以手动枚举的四个标量不等式(我已经将它们命名为a1、a2、z1、z2)：

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

这表示一个(可能是空的)多面体，也就是2d空间中四个半空间的交点。

第三:解决二维问题

(编辑:为清晰起见更新此位)

那么具体的2D问题是什么呢？我们希望找到一个z，其中实现了以下解决方案之一(这些解决方案并不完全不同，但这并不重要)：

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

如果其中一个实现了，我们有一个解决方案:z，其中z和Bz+q的元素极小值是0向量。如果这四个问题都无法实现，那么这个问题是不可行的，我们将宣布不存在这样的z。

第一种情况发生在a1，a2的交点为正正态分布时。只需选择解z= B^-1q，并检查它是否是元素非负的。如果是，返回B^-1q (注意，即使B不是psd，我假设它是可逆/满秩)。所以：

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

第二种情况(与第一种情况不完全不同)发生在a1、a2的交点在任何地方但包含原点时。换句话说，无论何时Bz+q >0 for z= 0。当q为正元素时，就会发生这种情况。所以：

elif q is elementwise nonnegative, return z = 0.

第三种情况是z1=0，当z2 = -q2/B22时，它可以被替换成a2来表示a2=0。z2必须是>=0，so -q2/B22 >=0.a1还必须是>=0，它在值中替代z1和z2，给出-B22/B12*q2 + q1 >=0。所以：

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

第四步与第三步相同，但交换1和2。

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

最后的第五种情况是问题不可行的时候。当不满足上述条件时，就会发生这种情况。所以：

else return None 

最后:把这些片段放在一起，

求解每个二维问题是一对简单/高效/平凡的线性求解和比较。每个人将返回一对数字或None。然后对所有的n 2D问题做同样的处理，并将向量连接起来。如果没有，这个问题是不可行的(全都没有)。否则你就有答案了。

Answer 2

基于AMPLPY的Python LCP求解器

正如@denfromufa所指出的，PATH解算器有一个PATH接口。他建议使用Pyomo，因为它是开源的，能够处理AMPL。然而，Pyomo的使用却是缓慢而乏味的。我最终在cython中为PATH解决程序编写了自己的界面，并希望在某一时刻发布它，但此时它没有输入卫生设施，速度快且不干净，我不认为有时间去做它。

现在，我可以共享一个使用AMPL的python扩展的答案。它的速度不如PATH的直接接口:对于我们想要解决的每个LCP，它创建一个(临时)模型文件，运行AMPL，并收集输出。这有点快速和肮脏，但我觉得我至少应该报告过去几个月的结果，因为我问了这个问题。

import os
# PATH license string needs to be updated
os.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, dataframe
import numpy as np
from scipy import sparse
from tempfile import mkstemp
import os

import sys
import contextlib

class DummyFile(object):
    def write(self, x): pass

@contextlib.contextmanager
def nostdout():
    save_stdout = sys.stdout
    sys.stdout = DummyFile()
    yield
    sys.stdout = save_stdout


class 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))