python - 具有稀疏矩阵的 LCP

Question

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

我需要为向量求解以下线性不等式系统 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 ， 一世

创建四个矩阵 A11 , A12 , A21 , A22使用 scipy.sparse.diags

并将它们堆叠在一起作为 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 从 Python 调用它存在问题。 ( my issue report here )

此外，Matlab 有一个 LCP 求解器，显然可以处理稀疏矩阵，但实现更加古怪(我真的不想为此退回 Matlab)

Answer 1

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

第零:澄清问题/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和 z_{n+1}出现 只有在方程式和术语中，而不是在其他地方。所以问题可分为 n问题，每个问题都是二维的。因此，我以后将解决二维问题，您可以在 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

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

第三:解决二维问题

(编辑:为了清楚起见更新了这一点)

那么具体的二维问题是什么？我们要找一个 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

如果实现其中之一，我们有一个解决方案: a 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 

最后:将碎片拼凑起来

解决每个 2D 问题是几个简单/高效/微不足道的线性求解和比较。每个将返回一对数字或 None .然后对所有 n 做同样的事情2D 问题，并连接向量。如果有任何是 None，则问题不可行(所有 None)。否则，你有你的答案。