Ошибка решателя Mosek при добавлении ограничений к задаче оптимизации (переменная 10000, с использованием Python / cvxpy)

Суммируя

Приведенная ниже проблема оптимизации объявляется невыполнимой при запуске с Mosek, но решается (легко и точно) с помощью решающей программы с открытым исходным кодом ECOS.

Мне интересно: почему такой продвинутый коммерческий решатель, как Mosek, не может решить эту проблему?

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. Я загрузил пробную лицензию и использую Mosek v8.1 и cvxpy 0.4.10.

Я обнаружил, что Mosek, похоже, не очень точно придерживается ограничений или терпит неудачу в случаях со многими ограничениями. Вот в чем я хотел бы помочь - почему Mosek неточно описывает эти ограничения и почему он не работает для четко решаемых проблем?

В приведенном ниже сценарии я решаю простую задачу всего с двумя ограничениями (все переменные положительны, сумма которых равна 1), а затем повторно решаю ее с несколькими добавленными ограничениями, которые я называю ограничениями «суммирования».

Все эти добавленные ограничения имеют форму «веса для некоторого подмножества переменных должны в сумме быть меньше, чем 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 дополнительных ограничений. Эти ограничения суммирования имеют вид 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, что, на мой взгляд, недостаточно точно.

Затем, когда дело доходит до решения более ограниченной проблемы, Mosek полностью терпит неудачу и объявляет об этом PRIMAL_INFEASIBLE.

Может ли кто-нибудь предложить какие-либо идеи относительно того, почему Мосек терпит неудачу? Я видел, что это было крайне неточно в отношении ограничений и в других случаях. Я попытался повысить точность с помощью параметра intpnt_co_tol_pfeas, но всякий раз, когда я это делаю, решающая программа часто дает сбой.

Заранее благодарим за любую помощь с этим. Вот мой пример кода. Запуск с solver='ECOS' работает, запуск с solver='MOSEK' не работает.