具有稀疏矩阵的LCP-Java 学习之路

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

我需要为vector v 解决以下线性不等式系统：

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

其中 0 是零向量 .

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

定义 z = v-s ， B=rI - A ， q=-u + Bs ，我可以将上一个问题重写为linear complementarity problem，并希望使用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 实现适合这个问题？

我真的不认为一般"which tools to solve 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稀疏求解器都要求问题是凸的，而我的不是
在朱莉娅，PATHSolver可以解决我的问题（给定许可证） . 但是，使用PyJulia（my issue report here）从Python调用它时会出现问题
Matlab也有一个LCP求解器，显然可以处理稀疏矩阵，但实现更加古怪（我真的不想在Matlab上回退这个）

2 回答

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

Zeroth: clarifying the problem / LCP

根据评论中的澄清，@ FooBar说最初的问题是元素 min ;我们需要找到 z （或 v ）这样的

左参数为零，右参数为非负数，左参数为非负数，右参数为零

正如@FooBar正确指出的那样，有效的重新参数化导致了LCP . 但是，下面我展示了可以实现原始问题的更简单，更有效的解决方案（通过利用原始问题中的结构）而无需LCP . 为什么这会更容易？好吧，请注意，LCP在 z （Bz q）'z, but that the original problem doesn' t中只有一个二次项（只有线性项Bz q和z） . 我将在下面利用这个事实 .

First: simplify

有一个重要但关键的细节以一种主要方式简化了这个问题：

使用scipy.sparse.diags创建四个矩阵A11，A12，A21，A22并将它们堆叠在一起作为A = scipy.sparse.bmat（[[A11，A12]，[A21，A22]]）

这具有重大意义 . 具体来说，这不是问题，而是 really small （2D，准确地说）问题 . 请注意，此 A 矩阵的块对角线结构在所有后续操作中都会保留 . 并且每个后续操作都是矩阵向量乘法或内积 . 这意味着该程序实际上可以成对分离 z （或 v ）个变量 .

具体而言，假设每个块 A11,... 的大小为 n n . 然后批判性地注意 z_1 和 z_{n+1} 在他们自己的方程式和术语中出现 only ，而从未在其他地方出现过 . 所以问题可分为 n 问题，每个问题都是二维的 . 因此，我将在此后解决2D问题，您可以根据需要序列化或并行化 n ，无需稀疏矩阵或大型选择包 .

Second: the geometry of the 2D problem

假设我们在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空间中四个半空间的交集 .

Third: solving the 2D problem

（编辑：为了清晰起见，更新了这一位）

具体是什么2D问题呢？我们希望找到一个 z ，其中包含以下解决方案之一（并非所有不同，但无关紧要）：
- a1> = 0，z1 = 0，a2> = 0，z2 = 0
- a1 = 0，z1> = 0，a2 = 0，z2> = 0
- a1> = 0，z1 = 0，a2 = 0，z2> = 0
- a1 = 0，z1> = 0，a2> = 0，z2 = 0
如果其中一个实现，我们有一个解决方案：一个z在哪里z和Bz q的元素最小值是0矢量 . 如果这四个都不能实现，那么这个问题是不可行的，我们将宣布不存在这样的问题 z .

第一种情况发生在a1，a2的交点处于正orthant时 . 只需选择解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，当z2 = -q2 / B22时，可以将其代入a2以表示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
```
Finally: put the pieces together

解决每个2D问题是一些简单/有效/平凡的线性求解和比较 . 每个都会返回一对数字或 None . 然后在所有 n 2D问题上做同样的事情，并连接矢量 . 如果有的是None，则问题是不可行的（全部无） . 否则，你有答案 .
回复于 2024-05-04T13:00:57+08:00

基于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))

回复于 2024-05-04T13:00:57+08:00

具有稀疏矩阵的LCP

2 回答

基于AMPLPY的Python LCP求解器

相关问题