如何用Python解决整数优化问题？

问题描述 投票：0回答：1

我有一个线性方程组 Ax=b，A - 矩阵 (m x n)，x - (n x 1)，b - (m x 1)。矩阵 A 几乎全为 0，只包含 ~ n*m^(1/2) 个非零元素。 m 和 n 可以很大，例如 n=28680，m=22500 或 n=28680，m=62500。很明显，系统的高概率精确解根本不存在。因此，我需要找到实现 sum( (Ax-b)^2 ) 函数的最小值的向量 x0。同时，非常希望向量 x0 仅由 0 和 1 组成。因此，我需要解决问题 sum( (Ax-b)^2 ) -> min, x 来自 {0,1}

你是如何尝试解决这个问题的？ 首先，我解决了一个简化的问题。 sum( (Ax-b)^2 ) -> min, x 来自实数 使用tensoflow方法很快就解决了这个问题。我汇总了最终的解决方案，但这导致质量很差。

然后我试图解决完整的问题。向我推荐了 CVXPY 库。通过研究 CVXPY 官方网站，我发现它适合我对混合整数二次类使用求解器。我使用 SCIP 求解器编写了一个程序。由于等待最优解的时间太长，我将求解时间限制在 200 秒。该程序有效，但仅适用于较小的 m 和 n，例如 n=4950，m=1600。如果进一步增加 m 和 n，则 SCIP 会给出一个简单的零向量作为解决方案，尽管实际上任何包含一个 1 的向量都更接近最佳值。这是程序代码：

import cvxpy as cp
import numpy as np

#I tried to simulate my inputs for you. 
#Here the problem manifests itself already at m, n = 1600, 4950.
np.random.seed(0)
m, n= 1600, 4950
A = np.round(np.random.rand(m, n)/(1-1/m**(1/2))/2)*255*(m**(1/2)/1800)**1.1
b = 255*(np.random.randn(m))

x = cp.Variable(n, integer=True)
objective = cp.Minimize(cp.sum_squares(A @ x - b))
constraints = []
constraints.extend([x >= 0.0, x <= 1])
prob = cp.Problem(objective,constraints)
prob.solve(solver=cp.SCIP,verbose=True,scip_params={"limits/time": 50})
print("Status: ", prob.status)
print("The optimal value is", prob.value)
print("A solution x is")
print(x.value)
print("Non-zero ",sum([i for i in x.value]))

我的问题如下：

  1. 有没有更好的选择来解决我的问题？
  2. 也许我应该使用其他求解器？我找不到另一个免费的解决方案来解决我的问题。
  3. 如果 SCIP 是最好的选择，我怎样才能让它适用于大的 m 和 n？我不明白为什么它停止工作...
python mathematical-optimization pyomo cvxpy scip
1个回答
0
投票

我假设“绝对误差之和”足以作为解决方案而不是 SSE 的质量衡量标准。通过这样做，模型可以很容易地线性化，我认为 MIP 求解器比二次求解器更容易处理。请注意，下面的两个示例都assume有一个good（在这种情况下实际上是精确的）解决方案。如果矩阵基本上是“噪声”并且您正在寻找最合适的矩阵，那么求解器几乎肯定会更加困难，因为问题可能更加退化。

在查看解决方案之前，我认为您对上述

A
的构建存在问题。做一个小的 m x n 矩阵，看看它是否是你所期望的。我得到了很多重复的，所有的正值，所以 x 向量 的“全零”完全有可能是 b
 向量的最佳解决方案
，应该平均为 ~0。我用正/负值和结构做了一些改变，以获得你想要的密度。

如果您使用我的 

cvxpy
shell 并将其切换回 SSE，我很想知道您获得的性能与线性性能相比如何。我无法找到合适的求解器来解决这个问题，因为我不常使用
cvxpy
。

cvxpy
对于像这样可以很容易地用矩阵格式表达的问题来说是很自然的。我在 
pyomo
中跳过几个圈以在解决之前到达同一点，但它工作得很好。在我 8 岁的时候。旧机器，
pyomo
需要大约 200 秒来构建最大尺寸的问题（如您在结果中看到的那样）

my 

cvxpy
下面的解决方案适用于较小的实例，但 barfs 对于 20k x 20k 左右的矩阵“不可行”，所以可能有一些我不理解的内部内容，因为它显然是可行的。也许对
cvxpy
有更好了解的人可以修补并看看发生了什么。

CVXPY 用于小矩阵：

import cvxpy as cp
import numpy as np
from time import time

#I tried to simulate my inputs for you. 
#Here the problem manifests itself already at m, n = 1600, 4950.
np.random.seed(0)
m, n= 5_000, 5_000  # 20_000, 20_000  <--- barfs
A_base = (np.random.randint(-1000, 1000, (m, n)))
marks = np.random.binomial(n=1, p=(n*m)**0.5/(m*n), size=(m, n))
A = np.where(marks, A_base, 0)
true_x = np.random.binomial(n=1, p=0.6, size=n)
b = A @ true_x

tic = time()
x = cp.Variable(n, integer=True)
objective = cp.Minimize(cp.sum(cp.abs(A @ x - b)))
constraints = []
constraints.extend([x >= 0.0, x <= 1])
prob = cp.Problem(objective,constraints)
res=prob.solve() #solver=cp.SCIP,verbose=True,scip_params={"limits/time": 50})
toc = time()
print("Status: ", prob.status)
print("The optimal value is", prob.value)
print("A solution x is")
print(x.value)
print("Non-zero ",sum([i for i in x.value if i is not None ]))
print(f'solver time:  {toc-tic : 0.1f} sec.')

CVXPY输出：

Status:  optimal
The optimal value is 0.0
A solution x is
[ 1.  1.  1. ... -0.  1. -0.]
Non-zero  1907.0
solver time:   3.5 sec.
[Finished in 5.5s]

PYOMO 用于全尺寸矩阵

import numpy as np
import pyomo.environ as pyo
from time import time
import sys

np.random.seed(0)
m, n = 28_680, 62_500
A_base = (np.random.randint(-1000, 1000, (m, n)))
marks = np.random.binomial(n=1, p=(n*m)**0.5/(m*n), size=(m, n))
A = np.where(marks, A_base, 0)
true_x = np.random.binomial(n=1, p=0.6, size=n)
b = A @ true_x
# b = np.random.randint(-100, 100, n)
# print(b)

tic = time()
model = pyo.ConcreteModel()

model.N = pyo.Set(initialize=range(n))
model.M = pyo.Set(initialize=range(m))

# initializing x only because some values of x[n] may end up unconstrained and this prevents "None" in result
model.x           = pyo.Var(model.N, within=pyo.Binary, initialize=0)
model.row_err     = pyo.Var(model.M)
model.row_abs_err = pyo.Var(model.M)

# constrain the rowsum/error appropriately
@model.Constraint(model.M)
def rowsum(model, m):
    return model.row_err[m] == sum(A[m,n] * model.x[n] for n in model.N if A[m,n]) - b[m]


# constrain the abs error
@model.Constraint(model.M)
def abs_1(model, m):
    return model.row_abs_err[m] >= model.row_err[m]

@model.Constraint(model.M)
def abs_2(model, m):
    return model.row_abs_err[m] >= -model.row_err[m]

# minimize sum of absolute row errors
model.obj = pyo.Objective(expr=sum(model.row_abs_err[m] for m in model.M))
toc = time()
print(f'starting solver...  at time: {toc-tic: 0.1f} sec.')
solver = pyo.SolverFactory('cbc')
res = solver.solve(model)
toc = time()

print(f'total solve time:  {toc-tic: 0.1f} sec')
print(res)
# model.x.display()

x = np.array([pyo.value(model.x[n]) for n in model.N], dtype=int)

print(f'verifying solution:  {np.sum(A@x - b)}')

PYOMO输出：

starting solver...  at time:  228.2 sec.
solver time:   384.8 sec

Problem: 
- Name: unknown
  Lower bound: 0.0
  Upper bound: 0.0
  Number of objectives: 1
  Number of constraints: 40941
  Number of variables: 50488
  Number of binary variables: 30839
  Number of integer variables: 30839
  Number of nonzeros: 20949
  Sense: minimize
Solver: 
- Status: ok
  User time: -1.0
  System time: 123.18
  Wallclock time: 152.91
  Termination condition: optimal
  Termination message: Model was solved to optimality (subject to tolerances), and an optimal solution is available.
  Statistics: 
    Branch and bound: 
      Number of bounded subproblems: 21
      Number of created subproblems: 21
    Black box: 
      Number of iterations: 2668
  Error rc: 0
  Time: 153.19930219650269
Solution: 
- number of solutions: 0
  number of solutions displayed: 0

verifying solution:  0
最新问题
© www.soinside.com 2019 - 2024. All rights reserved.