简而言之
下面的优化问题在使用 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'失败。