Как решить проблему целочисленной оптимизации с помощью Python?

Я исходил из предположения, что «сумма абсолютных ошибок» достаточно хороша в качестве меры качества решения, а не SSE. Таким образом, модель может быть легко линеаризована, и я думаю, что гораздо более податливым является решатель MIP, а не квадратичный решатель. Обратите внимание, что в обоих приведенных ниже примерах предполагается, что существует хорошее (фактически точное в данном случае) решение. Если матрицы в основном «шумовые», и вы ищете наилучшее соответствие, то решатели почти наверняка будут испытывать гораздо большие трудности, поскольку проблема, вероятно, будет более вырожденной.

Прежде чем взглянуть на решение, я думаю, что есть проблемы с вашей конструкцией A выше. Создайте небольшую матрицу m x n и посмотрите, соответствует ли она вашим ожиданиям. Я получаю много дубликатов, все положительные значения, поэтому вполне возможно, что «все нули» для вектора x являются оптимальным решением для вектора b, который должен в среднем ~ 0. Я немного изменил значения положительных/отрицательных значений и конструкцию, чтобы получить желаемую плотность.

Если вы возьмете мою оболочку cvxpy и переключите ее обратно на SSE, мне будет любопытно посмотреть, какую производительность вы получите по сравнению с линейной. Я не мог заставить подходящий решатель работать для этого, так как я не так часто использую cvxpy.

cvxpy является естественным для таких задач, которые легко могут быть выражены в матричном формате. Я прыгаю через несколько обручей в pyomo, чтобы добраться до той же точки перед решением, но это работает довольно хорошо. В мои 8 лет. старая машина, pyomo требуется около 200 секунд, чтобы построить задачу максимального размера (как вы видите в результатах)

мое 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)}')

Выход ПИОМО:

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