Я предложил показать вам, как представить это как линейную программу, так что вот так. Я использую Pyomo, так как я более знаком с этим, но вы можете сделать что-то подобное в PuLP.
Чтобы запустить это, вам нужно сначала установить Pyomo и линейный программный решатель, такой как glpk. glpk должен работать для задач разумного размера, но если вы обнаружите, что решение занимает слишком много времени, вы можете попробовать (гораздо более быстрый) коммерческий решатель, такой как CPLEX или Gurobi.
Вы можете установить Pyomo через pip install pyomo или conda install -c conda-forge pyomo. Вы можете установить glpk из https://www.gnu.org/software/glpk/ или через conda install glpk. (Я думаю, что PuLP поставляется с встроенной версией glpk, так что это может сэкономить вам шаг.)
Вот сценарий. Обратите внимание, что это вычисляет абсолютную ошибку как линейное выражение, определяя одну переменную для положительного компонента ошибки и другую для отрицательной части. Затем он стремится минимизировать сумму обоих. В этом случае решатель всегда будет устанавливать один в ноль, так как это простой способ уменьшить ошибку, и тогда другой будет равен абсолютной ошибке.
import random
import pyomo.environ as po
random.seed(1)
# ~50% sparse data set, big enough to populate every row and column
m = 10 # number of rows
n = 10 # number of cols
data = {
(r, c): random.random()
for r in range(m)
for c in range(n)
if random.random() >= 0.5
}
# define a linear program to find vectors
# x in R^m, y in R^n, such that x[r] + y[c] is close to data[r, c]
# create an optimization model object
model = po.ConcreteModel()
# create indexes for the rows and columns
model.ROWS = po.Set(initialize=range(m))
model.COLS = po.Set(initialize=range(n))
# create indexes for the dataset
model.DATAPOINTS = po.Set(dimen=2, initialize=data.keys())
# data values
model.data = po.Param(model.DATAPOINTS, initialize=data)
# create the x and y vectors
model.X = po.Var(model.ROWS, within=po.NonNegativeReals)
model.Y = po.Var(model.COLS, within=po.NonNegativeReals)
# create dummy variables to represent errors
model.ErrUp = po.Var(model.DATAPOINTS, within=po.NonNegativeReals)
model.ErrDown = po.Var(model.DATAPOINTS, within=po.NonNegativeReals)
# Force the error variables to match the error
def Calculate_Error_rule(model, r, c):
pred = model.X[r] + model.Y[c]
err = model.ErrUp[r, c] - model.ErrDown[r, c]
return (model.data[r, c] + err == pred)
model.Calculate_Error = po.Constraint(
model.DATAPOINTS, rule=Calculate_Error_rule
)
# Minimize the total error
def ClosestMatch_rule(model):
return sum(
model.ErrUp[r, c] + model.ErrDown[r, c]
for (r, c) in model.DATAPOINTS
)
model.ClosestMatch = po.Objective(
rule=ClosestMatch_rule, sense=po.minimize
)
# Solve the model
# get a solver object
opt = po.SolverFactory("glpk")
# solve the model
# turn off "tee" if you want less verbose output
results = opt.solve(model, tee=True)
# show solution status
print(results)
# show verbose description of the model
model.pprint()
# show X and Y values in the solution
for r in model.ROWS:
print('X[{}]: {}'.format(r, po.value(model.X[r])))
for c in model.COLS:
print('Y[{}]: {}'.format(c, po.value(model.Y[c])))
Просто чтобы завершить рассказ, вот решение, которое ближе к вашему первоначальному примеру. Он использует cvxpy, но с моим подходом разреженных данных.
Я не знаю "официального" способа выполнения поэлементных вычислений с помощью cvxpy, но, похоже, все в порядке, если использовать стандартную функцию Python sum с большим количеством отдельных cp.abs(...) вычислений.
Это дает решение, которое немного хуже, чем линейная программа, но вы можете исправить это, отрегулировав допуск решения.
import cvxpy as cp
import random
random.seed(1)
# Problem data.
# ~50% sparse data set
m = 10 # number of rows
n = 10 # number of cols
data = {
(i, j): random.random()
for i in range(m)
for j in range(n)
if random.random() >= 0.5
}
# Construct the problem.
x = cp.Variable(m)
y = cp.Variable(n)
objective = cp.Minimize(
sum(
cp.abs(x[i] + y[j] + data[i, j])
for (i, j) in data.keys()
)
)
prob = cp.Problem(objective)
result = prob.solve()
print(x.value)
print(y.value)