Foo*_*Bar 53 python scipy sparse-matrix
我用大写字母表示矩阵,用小写字母表示向量.
我需要解决以下矢量线性不等式系统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使用scipy.sparse.diagsA = 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)这个问题有一个非常有效(线性时间)的解决方案,虽然它需要一些讨论......
Zeroth:澄清问题/ 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}出现在他们自己的方程和术语中,而不是其他地方.所以问题可分为问题,每个问题都是二维的.因此,我将在此后解决2D问题,您可以根据需要进行序列化或并行化,无需稀疏矩阵或大型opt包.nn
第二: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
这代表一个(可能是空的)多面体,也就是2d空间中四个半空间的交集.
第三:解决2D问题
(编辑:为了清楚起见,更新了这一位)
具体是什么2D问题呢？我们希望找到z以下解决方案之一(并非所有区别,但无关紧要):
如果实现了其中之一,我们有一个解决方案:az,其中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
最后:把各个部分放在一起
解决每个2D问题是一些简单/有效/平凡的线性求解和比较.每个都会返回一对数字或None.然后在所有n2D问题上做同样的事情,并连接矢量.如果有的是None,则问题是不可行的(全部无).否则,你有答案.
正如@denfromufa 指出的，解算器有一个AMPL接口PATH。他建议，Pyomo因为它是开源的并且能够处理AMPL. 然而，Pyomo结果是工作缓慢且乏味。我最终在 cython 中编写了自己的求解器接口PATH，并希望在某个时候发布它，但目前它没有输入卫生，又快又脏，我看不到处理它的时间。
现在，我可以分享一个使用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))