Короче говоря
Приведенная ниже проблема оптимизации объявлена не решаемой при запуске с Mosek, но разрешима (легко и точно) с помощью решателя с открытым исходным кодом ECOS.
Мне интересно: почему такой продвинутый коммерческий решатель, как Мозек, не может решить эту проблему ?
import cvxpy as cvx
import numpy as np
print('cvxpy version:')
print(cvx.__version__)
print('')
np.random.seed(0)
SOLVER = 'ECOS_BB' # Works fine, sticks to constraint thresholds very precisely
# SOLVER = 'MOSEK' # Fails when many "sumproduct" constraints are added
def get_objective_function_and_weights(n, means, std_devs):
weights = cvx.Variable(n)
# Markowitz-style objective "expectation minus variance" objective function
objective = cvx.Maximize(
means * weights
- cvx.sum_entries(cvx.mul_elemwise(std_devs, weights) ** 2)
)
return objective, weights
def calculate_objective_value(weights, means, std_devs):
expec = weights.T @ means
cov = np.power(np.diag(std_devs), 2)
var = weights.T @ cov @ weights
return float(expec - var)
def get_total_discrepancy(weights, A, b):
# We want A @ weights <= b
# Discrepancy is sum of any elements where this is not the case
values = A @ weights
assert values.shape == b.shape
discrepancy_idx = values > b
discrepancies = values[discrepancy_idx] - b[discrepancy_idx]
return discrepancies.sum()
def get_weights_gt_0_discrepancy(weights):
discrepancy_idx = weights < 0
discrepancies = np.abs(weights[discrepancy_idx])
return discrepancies.sum()
def get_sum_weights_le_1_discrepancy(weights):
discrepancy = max(0, weights.sum() - 1)
return discrepancy
def main():
n = 10000
means = np.random.normal(size=n)
std_devs = np.random.chisquare(df=5, size=n)
print()
print(' === BASIC PROBLEM (only slightly constrained) ===')
print()
objective, weights = get_objective_function_and_weights(n, means, std_devs)
constraints = [
weights >= 0,
cvx.sum_entries(weights) == 1,
]
# Set up problem and solve
basic_prob = cvx.Problem(objective, constraints)
basic_prob.solve(solver=SOLVER, verbose=True)
basic_weights = np.asarray(weights.value)
print('Optimal weights')
print(basic_weights.round(6))
print('Objective function value:')
basic_obj_value = calculate_objective_value(basic_weights, means, std_devs)
print(basic_obj_value)
print('Discrepancy: all weights > 0')
print(get_weights_gt_0_discrepancy(basic_weights))
print('Discrepancy: sum(weights) <= 1')
print(get_sum_weights_le_1_discrepancy(basic_weights))
print()
print()
print()
print(' === CONSTRAINED PROBLEM (many added "sumproduct" constraints) ===')
print()
objective, weights = get_objective_function_and_weights(n, means, std_devs)
# We will place all our sumproduct constraints into a single large matrix `A`
# We want `A` to be a fairly sparse matrix with only a few 1s, mostly 0s
m = 100 # number of "sumproduct" style constraints
p = 5 # on average, number of 1s per row in `A`
A = 1 * (np.random.uniform(size=(m, n)) < p/n)
# We look at the observed values of A @ weights from the basic
# problem, and base our constraint on that
observed_values = (A @ basic_weights).reshape((-1, 1))
b = observed_values * np.random.uniform(low=0.90, high=1.00, size=(m, 1))
print('number of 1s in A')
print(A.sum())
new_constraints = [
weights >= 0,
cvx.sum_entries(weights) == 1,
A * weights <= b,
]
# Set up problem and solve
prob = cvx.Problem(objective, new_constraints)
prob.solve(solver=SOLVER, verbose=True)
final_weights = np.asarray(weights.value)
print('Optimal weights')
print(final_weights.round(6))
print('Objective function value:')
constrained_obj_value = calculate_objective_value(final_weights, means, std_devs)
print(constrained_obj_value)
print('Difference vs basic')
print(constrained_obj_value - basic_obj_value)
# Now calculate the discrepancy - the amount by which the returned
# optimal weights go beyond the required threshold
print('Discrepancy: all weights > 0')
print(get_weights_gt_0_discrepancy(final_weights))
print('Discrepancy: sum(weights) <= 1')
print(get_sum_weights_le_1_discrepancy(final_weights))
print('Discrepancy: sumproduct threshold:')
print(get_total_discrepancy(final_weights, A, b))
main()
_
Подробнее
Я тестируюнекоторые оптимизаторы, и смотрели на Мосека.Я скачал пробную лицензию и использую Mosek v8.1 и cvxpy 0.4.10.
Я обнаружил, что Mosek, похоже, не очень точно придерживается ограничений или не работает в случаях со многимиограничения.Вот с чем мне бы хотелось помочь - , почему Mosek неточен в этих ограничениях и почему он не работает для четко решаемых задач ?
В приведенном ниже сценарии я решаю простую проблему с помощью толькодва ограничения (все переменные положительные, суммирующие до 1), затем повторно решите его с несколькими добавленными ограничениями, которые я называю «sumproduct» ограничениями.
Эти добавленные ограничения имеют вид "веса для некоторого подмножества переменных должны быть меньше чем b_i" для некоторой константы b_i.Я упаковываю эти ограничения в матричное уравнение A @ weights <= b.
Когда я запускаю скрипт внизу этого поста, используя встроенный решатель ECOS, он легко решает основную проблему, давая оптимальное значение 2,63....:
Objective function value:
2.6338492447784283
Discrepancy: all weights > 0
4.735618828548444e-13
Discrepancy: sum(weights) <= 1
1.3322676295501878e-15
Вы видите, что я также вычисляю то, что я называю несоответствием для каждого ограничения.Это величина, на которую оптимизатор «преодолевает» ограничение в возвращаемых весах.Таким образом, ECOS здесь немного «нарушает правила», определенные ограничениями, но ненамного.
Затем я прошу ECOS решить гораздо более ограниченную проблему с добавлением 100 дополнительных «sumproduct» ограничений.Эти ограничения на конечный продукт имеют вид A @ weights <= b, а A имеет 486, остальные нули.
number of 1s in A
486
Затем я заново решаю проблему и вижу пересмотренный набор оптимальных весов.Оптимальное значение немного меньше, чем было раньше (из-за добавленных ограничений), и ECOS все еще «соблюдает» все ограничения с очень хорошей точностью:
Objective function value:
2.6338405110044203
Difference vs basic
-8.733774008007344e-06
Discrepancy: all weights > 0
5.963041247103521e-12
Discrepancy: sum(weights) <= 1
9.103828801926284e-15
Discrepancy: sumproduct threshold:
0.0
Если я запускаю то же самоеСценарий с Mosek, я нахожу, что по основной проблеме, Mosek решает ее, но уже довольно далеко от одного из ограничений:
Objective function value:
2.633643747862593
Discrepancy: all weights > 0
7.039232392536552e-06
Discrepancy: sum(weights) <= 1
0
Это означает, что у нас есть несколько весов меньше нуля Суммирование в -7e-6, что недостаточно для моего вкуса.
Затем, когда дело доходит до решения более сложной проблемы, Мосек полностью терпит неудачу и объявляет ее PRIMAL_INFEASIBLE.
Может кто-нибудь предложить какие-либо идеи относительно того, почему Мосек терпит неудачу, как это?Я видел, что это было слишком неточно в отношении ограничений и в других случаях.Я пытался повысить точность, используя параметр intpnt_co_tol_pfeas, но всякий раз, когда я делаю это, решатель начинает часто выходить из строя.
Заранее спасибо за любую помощь в этом.Вот мой пример кода.Работа с solver='ECOS' работает, работа с solver='MOSEK' завершается неудачей.