Подходит ли моя проблема для выпуклой оптимизации, и если да, как выразить ее с помощью cvxpy?

Я предлагал показать вам, как представить это в виде линейной программы, так что поехали. Я использую 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)

Я не понял, но просто некоторые хакерские штуки, основанные на предположении:

• вам нужен cvxpy-эквивалент поведения правил широковещательной передачи numpy для массивов (m, 1) + (1, n)

Так что по-умному:

m = 3
n = 4
np.random.seed(1)

a = np.random.randn(m, 1)
b = np.random.randn(1, n)

a
array([[ 1.62434536],
   [-0.61175641],
   [-0.52817175]])

b
array([[-1.07296862,  0.86540763, -2.3015387 ,  1.74481176]])


a + b
array([[ 0.55137674,  2.48975299, -0.67719333,  3.36915713],
   [-1.68472504,  0.25365122, -2.91329511,  1.13305535],
   [-1.60114037,  0.33723588, -2.82971045,  1.21664001]])

Давайте имитируем это с помощью np.kron, в котором эквивалент cvxpy:

aLifted = np.kron(np.ones((1,n)), a)
bLifted = np.kron(np.ones((m,1)), b)

aLifted
array([[ 1.62434536,  1.62434536,  1.62434536,  1.62434536],
   [-0.61175641, -0.61175641, -0.61175641, -0.61175641],
   [-0.52817175, -0.52817175, -0.52817175, -0.52817175]])

bLifted
array([[-1.07296862,  0.86540763, -2.3015387 ,  1.74481176],
   [-1.07296862,  0.86540763, -2.3015387 ,  1.74481176],
   [-1.07296862,  0.86540763, -2.3015387 ,  1.74481176]])

aLifted + bLifted
array([[ 0.55137674,  2.48975299, -0.67719333,  3.36915713],
   [-1.68472504,  0.25365122, -2.91329511,  1.13305535],
   [-1.60114037,  0.33723588, -2.82971045,  1.21664001]])

Давайте полуслепо проверим cvxpy (мы только размеры; слишком ленив, чтобы настроить проблему и исправить переменную для проверки вывода :-D):

import cvxpy as cp
x = cp.Variable((m, 1))
y = cp.Variable((1, n))

cp.kron(np.ones((1,n)), x) + cp.kron(np.ones((m, 1)), y)
# Expression(AFFINE, UNKNOWN, (3, 4))

# looks good!

Теперь несколько предостережений:

• i don't know how efficient cvxpy can reason about this matrix-form internally
  • unclear if more efficient as a simple list-comprehension based form using cp.vstack and co (it probably is)
• this operation itself kills all sparsity
  • (if both vectors are dense; your matrix is dense)
  • cvxpy and more or less all convex-optimization solvers are based on some sparsity assumption
    • scaling this problem up to machine-learning dimensions will not make you happy
• для вашей проблемы, вероятно, существует гораздо более краткая математическая теория, чем использовать (с учетом разреженности) (довольно) общий (DCP, реализованный в cvxpy, является подмножеством) выпуклая оптимизация