python - 将约束添加到优化问题时 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）解决了一个简单的问题，然后用几个添加的约束重新解决它，我称之为“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”约束。这些 sumproduct 约束的形式为A @ weights <= b，并且A有 486 个 1，其余为 0。

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。

任何人都可以提供任何关于为什么 Mosek 会失败的想法吗？我也看到它在其他情况下的约束也非常不准确。我尝试使用参数来提高精度intpnt_co_tol_pfeas，但是每当我这样做时，求解器就会经常失败。

提前感谢您对此的任何帮助。这是我的示例代码。运行solver='ECOS'正常，运行solver='MOSEK'失败。

标签： pythonconvex-optimizationcvxpymosek

python - 将约束添加到优化问题时 Mosek 求解器失败（10000 个变量，使用 Python/cvxpy）

问题描述

简而言之

_

更多细节

解决方案

推荐阅读