稀疏矩阵LCP

Question

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

我需要解向量v的以下线性不等式组：

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

其中0是一个零向量。在

其中r是标量，u和{}是向量，A是矩阵。在

定义z = v-s，B=rI - A，q=-u + Bs，我可以将前面的问题重写为linear complementarity problem，并希望使用LCP解算器，例如来自openopt：

^{pr2}$

或者，在矩阵表示法中：

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

---

问题是我的方程组太庞大了。要创建A，我

• 使用scipy.sparse.diags创建四个矩阵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不是半正定的，所以我不能把线性互补问题（LCP）重新表述为凸二次问题（QP）
• 来自this exposition的所有QP稀疏解算器都要求问题是凸的，而我的不是
• 在Julia中，PATHSolver可以解决我的问题（只要有许可证）。但是，使用PyJulia（my issue report here）从Python调用它时出现问题
• 另外，Matlab有一个LCP解算器，显然可以处理稀疏矩阵，但它的实现更古怪（我真的不想为此放弃Matlab）

Answer 1

这个问题有一个非常有效的（线性时间）解决方案，尽管它需要一些讨论。。。在

Zeroth:澄清问题/LCP

根据评论中的澄清，@FooBar说原始问题是元素方面的min；我们需要找到一个z（或v），这样

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

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

首先：简化

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

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

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

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

第二：二维问题的几何结构

假设我们有二维的元素问题，即：

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

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

^{pr2}$

这代表一个（可能是空的）多面体，也就是二维空间中四个半空间的交集。在

第三：解决二维问题

（编辑：为了清楚起见，更新了这一点）

那么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表示z=0时。当q为元素正时发生这种情况。所以：

elif q is elementwise nonnegative, return z = 0.

第三种情况为z1=0，可以代入a2表示当z2=-q2/B22时a2=0。z2必须大于等于0，因此-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。然后对所有的n2D问题做同样的处理，并连接向量。如果有没有，这个问题是不可行的（全部没有）。否则，你有你的答案。在

Answer 2

基于AMPLPY

正如@denfromufa指出的，有一个AMPL接口到PATH解算器。他建议使用Pyomo，因为它是开源的，能够处理{}。然而，Pyomo结果是缓慢和乏味的工作。我最终在cython中编写了自己的接口到PATH解算器，并希望在某个时候发布它，但目前它没有输入卫生设施，速度快而且脏，而且我看不到时间来处理它。在

现在，我可以分享一个使用python扩展AMPL的答案。它不如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))