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，我执行以下操作:

• 使用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）问题。
• 所有来自 这个展示 的 QP 稀疏求解器都要求问题是凸的，而我的问题则不是。
• 在 Julia 中，PATHSolver 可以解决我的问题（如果有许可证）。但是，使用 PyJulia从 Python 调用它存在问题（我的问题报告在这里）。
• 此外，Matlab 有一个 LCP 求解器，显然可以处理稀疏矩阵，但是它们的实现更加奇怪（我真的不想因此回退到 Matlab）。

---

- FooBar

---

@FooBar 你需要提供样本输入数据，你甚至没有指定矩阵的大小。无论如何，尝试使用cvxpy而不是OpenOpt，或者更好地尝试使用scipy：https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.linprog.html - denfromufa

@FooBar，你尝试过Scipy中这两个稀疏矩阵最小二乘求解器吗？https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.lsqr.html#scipy.sparse.linalg.lsqr https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.lsmr.html#scipy.sparse.linalg.lsmr - denfromufa

显然，你可以使用PATH求解器与Pyomo一起使用，因为后者默认支持AMPL格式和Pyomo输出。https://stackoverflow.com/q/36809878/2230844 - denfromufa

@Leon提出了一个需要澄清的关键问题；实际问题陈述是含糊不清的。当您编写原始问题 min(rv - (u + Av), v - s) = 0时，您具体指的是什么：(a) min_v (rv - u - Av).dot(v-s)（如最小化二次型）; (b) 查找v，使得(min_v (rv - u - Av).dot(v-s)) == 0; 或者(c/d)上述两者之一但是按元素运算（不是内积.dot，而是按元素最小值，因此0是一个由零向量组成的向量？所有这些都有不同的答案; 有些相当好（有些不太好）。 - muskrat

@muskrat 后者，0 是一个向量。对于每个点，要么左参数为零且右参数为非负数，要么左参数为非负数且右参数为零。 - FooBar

Answer 1

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

第零步：澄清问题 / LCP

根据评论中的澄清，@FooBar表示原始问题是逐元素的min；我们需要找到一个z（或v），使得

左边的参数为零且右边的参数为非负数或左边的参数为非负数且右边的参数为零

正如@FooBar正确指出的那样，有效的重新参数化会导致LCP。 然而，下面我将展示如何在不需要LCP的情况下通过利用原始问题中的结构来实现更简单和更高效的解决方案。 为什么这样更容易？ 好吧，请注意，LCP在z中具有二次项（Bz + q）'z，但原始问题没有（仅具有线性项Bz + q和z）。 我将在下面利用这个事实。

第一步：简化

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

• 使用scipy.sparse.diags创建四个矩阵A11、A12、A21、A22
• 并将它们堆叠在一起，如A = scipy.sparse.bmat([[A11，A12]，[A21，A22]])

这有巨大的影响。 具体来说，这不是一个单一的大型问题，而是一系列非常小的（确切地说是2D）问题。 请注意，此A矩阵的块对角线结构在所有后续操作中都得以保留。 每个后续操作都是矩阵向量乘积或内积。 这意味着实际上这个程序在z（或v）变量对的情况下是可分离的。

具体来说，假设每个块A11，...的大小为n乘n。 然后，请注意关键的一点，即z_1和z_{n+1}仅出现在自身的方程和术语中，从未出现在其他地方。 因此，该问题可以分解为n个问题，每个问题都是2维的。 因此，我将在此之后解决2D问题，您可以根据需要对n进行串行化或并行化，而无需使用稀疏矩阵或大型opt软件包。

第二步：2D问题的几何

假设我们在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

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

第三步：解决二维问题

那么具体的二维问题是什么呢？我们要找到一个 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 不是半正定的，我也假设它是可逆/满秩的）。所以：

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

第二种情况（与第一种情况并不完全不同）发生在a1，a2的交点位于任何位置但确实包含原点的情况下。换言之，每当z = 0时Bz+q>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 

最后：将所有部分整合起来

解决每个2D问题都是几个简单/高效/平凡的线性求解和比较。每次求解都会返回一对数字或None。然后在所有n个2D问题上执行相同操作，并连接向量。如果有任何一个是None，则该问题无法解决（全部为None）。否则，您就得到了答案。

- muskrat

---

我并没有声称最初的问题 min(...) 是一个 LCP。事实上，只有重构才是LCP。关于PSDB，来自维基百科的解释是：“如果M是正定的，则任何用于求解（严格）凸QPs的算法都可以解决LCP。”。我认为这个声明隐含地意味着LCP确实**不依赖于B是否为PSD。 - FooBar

我完全同意；我并不是想暗示你声称了什么。我的意思是，原始问题可以非常高效地解决，而无需担心lcp。至于第二点，是的，你是对的，lcp不一定是psd M，但我的意思是许多解决方案都依赖于M是psd。 - muskrat

@FooBar 我编辑了我的答案中的LCP部分，以纠正我之前的陈述错误；感谢你的指出，希望这能澄清问题。 - muskrat

我正在尝试慢慢地进行，以免错过任何东西，所以请原谅所有的“琐事”。接下来，我不确定我理解如何将4个A11，...矩阵堆叠使A成为块对角矩阵？ - FooBar

有没有任何理由按照这个特定的顺序测试不同的情况？否则我建议先测试最后一个情况，因为它需要求逆矩阵。当然，这只是一个2x2的矩阵，但如果可能的话，我认为没有理由不这样做。尤其是第二种情况测试起来会更快。 - RcCookie

此外，如果矩阵不可逆，其他情况是否仍然能得到有效的解？ - RcCookie

Answer 2

基于AMPLPY的Python LCP求解器

正如@denfromufa所指出的那样，PATH求解器有一个AMPL接口。他建议使用开源的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))

- FooBar