In short
The optimisation problem below is declared infeasible when run with Mosek, but is solvable (easily and accurately) using the open-source solver ECOS.
I'm wondering: why is such an advanced commercial solver like Mosek failing to solve this problem?
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()
_
More details
I'm testing out some optimizers, and have been looking at Mosek. I've downloaded a trial licence and am using Mosek v8.1, and cvxpy 0.4.10.
I'm finding that Mosek doesn't seem to stick to constraints very accurately, or fails in cases with many constraints. That's what I'd like help with - why is Mosek imprecise on these constraints, and why does it fail for clearly solvable problems?
In the below script I solve a simple problem with just two constraints (all variables positive, summing to 1), then re-solve it with several added constraints, which I call "sumproduct" constraints.
These added contraints are all of the form "the weights for some subset of variables must sum to less than b_i" for some constant b_i. I pack these constraints into a matrix equation A @ weights <= b.
When I run the script at the bottom of this post using the in-built solver ECOS, it solves the basic problem easily, giving an optimal value of 2.63...:
Objective function value:
2.6338492447784283
Discrepancy: all weights > 0
4.735618828548444e-13
Discrepancy: sum(weights) <= 1
1.3322676295501878e-15
You can see that I'm also calculating what I call the discrepancy for each constraint. This is the amount by which the optimiser is "going over" the constraint in the returned weights. So ECOS here is slightly "breaking the rules" defined by the constraints, but not by much.
I then ask ECOS to solve a much more constrained problem, with 100 added "sumproduct" constraints. These sumproduct constraints are of the form A @ weights <= b, and A has 486 ones, the rest zeros.
number of 1s in A
486
Then I re-solve the problem, and see a revised set of optimal weights. The optimal value is a little less than it was before (because of the added constraints), and ECOS is still "obeying" all the constraints to within a very good accuracy:
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
If I run the same script with Mosek, I find that on the basic problem, Mosek solves it, but is already quite far out on one of the constraints:
Objective function value:
2.633643747862593
Discrepancy: all weights > 0
7.039232392536552e-06
Discrepancy: sum(weights) <= 1
0
That means that we have several weights less than zero summing to -7e-6, which isn't accurate enough for my liking.
Then when it comes to solving the more constrained problem, Mosek is failing completely and declaring it PRIMAL_INFEASIBLE.
Could anyone offer any ideas as to why Mosek is failing like this? I've seen it being extremely inaccurate on constraints in other cases too. I tried to increase the precision using the parameter intpnt_co_tol_pfeas, but whenever I do this the solver starts failing often.
Thanks in advance for any help with this. Here's my example code. Running with solver='ECOS' works, running with solver='MOSEK' fails.
