当约束添加到优化问题时，Mosek解算器失败（10000变量，使用Python/cvxpy）

Question

简言之

下面的优化问题被宣布为不可行时，与Mosek运行，但可以解决（容易和准确）使用开源求解器生态系统。在

我想知道：为什么像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。我下载了一个试用版的许可证，使用的是MosekV8.1和CVXPY0.4.10。在

我发现Mosek似乎没有很准确地遵守约束，或者在有很多约束的情况下失败。这就是我想要的帮助-为什么Mosek在这些约束条件上不精确，为什么它在清晰可解的问题上失败？在

在下面的脚本中，我用两个约束（所有变量都是正的，求和为1）来解决一个简单的问题，然后用几个添加的约束重新解决它，我称之为“sumproduct”约束。在

这些附加的约束对于某些常数bˉi来说都是“某些变量子集的权重之和必须小于bˉi”的形式，我将这些约束打包成一个矩阵方程A @ weights <= b。在

当我使用内置的求解器ECOS运行这篇文章底部的脚本时，它很容易地解决了基本问题，给出了2.63的最佳值…：

^{pr2}$

您可以看到，我还在计算每个约束的差异。这是优化程序“越过”返回权重中的约束的量。因此，这里的ECOS有点“违反规则”所定义的约束，但不会太多。在

然后我要求ECOS解决一个更受约束的问题，添加100个“sumproduct”约束。这些sumproduct约束的形式是A @ weights <= b，A有486个1，其余为零。在

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

这意味着我们有几个权重小于0，加起来是-7e-6，这对我来说不够准确。在

在解决更具约束性的问题时，Mosek完全失败了，并声明它PRIMAL_INFEASIBLE。在

有人能解释一下莫塞克为什么会这样失败吗？在其他情况下，我也看到它在约束条件上极不准确。我试图使用参数intpnt_co_tol_pfeas来提高精度，但每当我这样做时，解算器常常会失败。在

在此之前感谢您的任何帮助。下面是我的示例代码。用solver='ECOS'运行有效，用solver='MOSEK'运行失败。在

Answer 1

可能有很多原因导致这两个代码使用不同的公差。比如说问题

x=0.x；x=0

可行吗？只有当你允许x是无穷大的时候才是。在其他的世界里，你的问题可能很严重。在

总的来说，我建议阅读

https://docs.mosek.com/9.0/pythonapi/debugging-log.html

尤其是关于解决方案的总结。如果这没有帮助，那么把你的问题和解决方案摘要放在Mosek google群组里。在