Originally posted by: anthzed

Hi,

I am coding a Long Short portfolio where Sum of Longs=1 and Sum of Shorts=-1.

To avoid a simultaneous Long and Short position I use a SOS1 constraint implemented as Big M.

For sample random data, the code woks sometimes and sometimes it fails.

Below is the code (in cvxpy) and sample data where the code does not work.

Help would be appreciated.

 

Regards, Anth

 

 

 

import numpy as  np
import cvxpy as cp

n=20

#ret=np.random.randn(n,1)
ret=(np.array([0.633847,0.0224574,0.787235,-1.45919,-0.858576,-0.90331,0.513264,-0.595323,0.452494,0.475518, \
     -1.02247,0.0798508,-0.120809,0.6326,-0.274057,-0.450731,0.709198,-0.196911,-0.101243,0.178535])).reshape(20,1)

lamb=5

Shorts=cp.Variable(n)
Longs=cp.Variable(n)

constraints=[]

constraints += [sum(Shorts)==-1]
constraints += [sum(Longs)==1]

constraints += [Shorts>=-1]
constraints += [Shorts<=0]

constraints += [Longs>=0]
constraints += [Longs<=1]

# Big M
M=10
I=cp.Variable(n,boolean=True)

constraints += [Longs<=M*I]
constraints += [-Shorts<=M*(1-I)]

obj=cp.Minimize(lamb*cp.sum_squares(Shorts+Longs)-ret.T*Shorts-ret.T*Longs)
prob=cp.Problem(obj,constraints)

obj11=prob.solve(solver=cp.CPLEX)

if prob.status=='optimal':
    ww=np.vstack((Shorts.value,Longs.value)).T
    wt=np.sum(ww,axis=1)
    summ=[prob.solver_stats.solve_time, sum(wt), sum(abs(wt)), obj11]
else:
    print(prob.status)

#CPLEXOptimizers
#DecisionOptimization

Can you please specify what "works" and "does not work" means? Does it crash? Does it produce wrong answers? Does it fail to solve?

Also, since not everybody is familiar with cvxpy it would help if you could export the models as SAV (or LP) and attach them here.

#CPLEXOptimizers
#DecisionOptimization

Originally posted by: anthzed

For the problem above I get:

CPLEX Error 1217:No solutions exist.

However, when I run using ECOS_BB solver, it finds a optimal solution.

 

Further, I randomly generate the vector ret, and sometimes I get a optimal solution with CPLEX.

I will need to learn SAV or LP, and then will resend the problem.

 

#CPLEXOptimizers
#DecisionOptimization

I wonder about this

constraints += [sum(Shorts)==-1]
constraints += [sum(Longs)==1]

That says that the short and long options must each sum to 1.

But then you also have

constraints += [Longs<=M*I]
constraints += [-Shorts<=M*(1-I)]

Which says you cannot pick a short and a long option at the same time. These are contradictory constraints, or am I missing something?

#CPLEXOptimizers
#DecisionOptimization

Originally posted by: anthzed

 

They do not contradict.

A individual asset can't have a long and short option at the same time. But at a portfolio level, the shorts sum to -1 and the longs sum to 1.

The sav file is attached and I generated it by the cvxpy command:

obj11=prob.solve(solver=cp.CPLEX,cplex_filename='Q:\LSops.sav')

 

#CPLEXOptimizers
#DecisionOptimization

Originally posted by: EdKlotz

I downloaded the SAV file and ran it with the CPLEX 12.9 interactive optimizer.   I didn't see any problems:

 

MIP - Integer optimal, tolerance (0.0001/1e-06):  Objective = -1.7131737013e-01
Current MIP best bound = -1.7133442092e-01 (gap = 1.70508e-05, 0.01%)
Solution time =   14.84 sec.  Iterations = 3689877  Nodes = 707149 (29)
Deterministic time = 3976.48 ticks  (267.88 ticks/sec)

 

While I might hope for a bit faster performance given that this is a small model, the optimization did finish.   When you indicate you got:

CPLEX Error 1217:No solutions exist.

 

which routine issues this error?   Is it the actual prob.solve method, or one of the solution querying routines that follow:

 

obj11=prob.solve(solver=cp.CPLEX)

if prob.status=='optimal':
    ww=np.vstack((Shorts.value,Longs.value)).T
    wt=np.sum(ww,axis=1)
    summ=[prob.solver_stats.solve_time, sum(wt), sum(abs(wt)), obj11]

 

 

#CPLEXOptimizers
#DecisionOptimization

Originally posted by: anthzed

So it then appears to be a cvxpy issue

#CPLEXOptimizers
#DecisionOptimization

Originally posted by: EdKlotz

That  probably the most likely  possibility, but I would not conclude that based on the information so far.    I see numerous other possibilities.

 

• cvxpy presumably calls CPLEX's Python API, so it remains possible that the issue is in the underlying CPLEX Python API.
• The model is pretty simple, and your code to build it seems pretty simple as well.   But it remains possible that you have something in your code that creates a problem that manifests itself with the no solution error you report.   I would suggest you create a really, really simple model with the cpxpy program.   Start with just one constraint and a handful of integer variables and a linear objective.   Then add a few convex quadratic terms to the objective.   Can you reproduce it there?  
• Take the SAV file you exported, and run it through the mipex2.py example that comes with your CPLEX distribution.   See what happens there.   At least then you are running a python program rather than interactive CPLEX, so the differences with your cvxpy program are fewer.
• It's always possible there's some issue with the specific Python version you are using.   Please provide the version and the Python developer on the machine in question

Let me also reiterate my question from my previous post: which function call in your code results in the no solution error?

Do these tests and post the results, plus any other ideas you have for tests that might help narrow this down.   If you are also familiar with the DoCplex Python API, you could probably build and solve your model pretty quickly with that, which would be a useful test.   If not, and the above tests don't reveal any useful information, post again to the web and one of us here will try to build your model with the DoCplex API and see if we can reproduce it that way.   If we can, then cvxpy is off the hook.   If we cannot, then, along with the other tests, it does look like cvxpy is involved in the problem.

 

#CPLEXOptimizers
#DecisionOptimization

Originally posted by: anthzed

The problem happens after prob.solve   (and not the post querying).

Python is 3.6

I ran it with mipex2.py and it worked.

Seems like cvxpy is the culprit.

I am not familiar with DoCplex so if you could code it I would really appreciate it (I kept the original problem as thin as possible).

#CPLEXOptimizers
#DecisionOptimization

Originally posted by: EdKlotz

Python 3.6 should be fine.    I did write a DoCplex version of your model, which is attached.  It did not result in the 1217 error, and cvxpy looks like the culprit.  I'll get into the program below, but first let me get back to the 1217 error you reported with your code.  I asked a colleague familiar with the cvxpy github repositories to look into this, and he concluded the 1217 errors is superfluous and harmless.   When you solve a MILP or MIQP, there is no associated LP, basis matrix, dual values, or various other info that pertains only to LPs.   The cvxpy API code appears to be calling an LP only solution query routine, which results in the 1217 error.   But this does not stop anything from proceeding; it just generates a confusing but superfluous error message.   You can check this by verifying that the following segment of your code still executes:

if prob.status=='optimal':
    ww=np.vstack((Shorts.value,Longs.value)).T
    wt=np.sum(ww,axis=1)
    summ=[prob.solver_stats.solve_time, sum(wt), sum(abs(wt)), obj11]

If that executes, you are getting solution values, which shows that a solution actually is available.  And if those code does not execute, what status do you get in the code just below it that does?

else:
    print(prob.status)

 

If we are talking about the same behavior, I bet you don't see any status here with a 1217 value.

 

Returning to the DoCplex version I have attached, it's a bit more general; you specify the return data in a file on the command line, and after that on the command line specify the multiplier for the quadratic objective term.    I ran it with your return data, and no 1217 error occurs.   I have attached both the program and your sample data; you would run it with

python longshort.py 5.0

to reproduce your data.   As a disclaimer, I am not an extremely Pythonic programmer, so if you see my code uses some form of for loop where you used a more vector oriented operation, do not assume that means DoCplex lacks the more vector oriented approach you used. 

 

I did learn one other thing that might be of interest regarding faster performance.   Your code builds the sum of squares of each long + short term for an asset as follows:

cp.sum_squares(Shorts+Longs)

In other words, you are taking the sum of squares of some linear expressions, which causes the model building code in the cvxpy API (and probably the DoCplex API as well, although I didn't test it) to have to create extra variables and constraints to model these linear expression.   It is unaware that in fact only one of any short or long activity for a particular asset can take on a nonzero value (and you should not expect that level of preprocessing by any modeling language or modeling API).   So when I look at the LP file of the model you attached to the forum, I see these constraints and variables, where x1,...,x20 appear as squared terms in the objective:

 

 + [ 10 x1 ^2
      + 10 x2 ^2 + 10 x3 ^2 + 10 x4 ^2 + 10 x5 ^2 + 10 x6 ^2 + 10 x7 ^2
      + 10 x8 ^2 + 10 x9 ^2 + 10 x10 ^2 + 10 x11 ^2 + 10 x12 ^2 + 10 x13 ^2
      + 10 x14 ^2 + 10 x15 ^2 + 10 x16 ^2 + 10 x17 ^2 + 10 x18 ^2 + 10 x19 ^2
      + 10 x20 ^2 ] / 2
Subject To
 c1:   - x1 + x21 + x22  = 0
 c2:   - x2 + x23 + x24  = 0
 c3:   - x3 + x25 + x26  = 0
 c4:   - x4 + x27 + x28  = 0
 c5:   - x5 + x29 + x30  = 0
 c6:   - x6 + x31 + x32  = 0
 c7:   - x7 + x33 + x34  = 0
 c8:   - x8 + x35 + x36  = 0
 c9:   - x9 + x37 + x38  = 0
 c10:  - x10 + x39 + x40  = 0
 c11:  - x11 + x41 + x42  = 0
 c12:  - x12 + x43 + x44  = 0
 c13:  - x13 + x45 + x46  = 0
 c14:  - x14 + x47 + x48  = 0
 c15:  - x15 + x49 + x50  = 0
 c16:  - x16 + x51 + x52  = 0
 c17:  - x17 + x53 + x54  = 0
 c18:  - x18 + x55 + x56  = 0
 c19:  - x19 + x57 + x58  = 0
 c20:  - x20 + x59 + x60  = 0

 

My formulation doesn't have this because I exploit the complementarity condition on the long/short activities when building the objective:

#                                                                                                                                                               
#   Constraints are done; build the objective                                                                                                                   
#                                                                                                                                                               
    objexpr = model.quad_expr()
    objexpr += lambd*(model.sumsq(Long) + model.sumsq(Short))
    objexpr -= (model.scal_prod(Long, retcoeffs) + model.scal_prod(Short, retcoeffs))
    model.set_objective("min", objexpr)

 

On other words for any asset [i], if you know that Long[i]*Short[i] = 0, then you know (Long[i] + Short[i])^2 = Long[i]^2 + Short[i]^2.   This permits the creation of sum of squares of individual variables rather than linear expressions, and avoids the need for the auxiliary variables and associated constraints in your model.    I mention this because the performance difference was surprisingly different.   With CPLEX 12.9 and 12 threads, your version of the model took 14-15 seconds and 700000 nodes, whereas my formulation took about .1 of a second:

 

MIP - Integer optimal solution:  Objective = -1.7131737013e-01
Solution time =    0.10 sec.  Iterations = 406  Nodes = 146
Deterministic time = 3.33 ticks  (35.00 ticks/sec)

 

I don't think I made a mistake in the formulation, but you can easily test this by modifying the way you build the sum of squares term to try this out.   This is independent of the Python modeling API used.    I'm still a bit puzzled why this would make such a huge difference (particularly regarding node count), and we'll look into that and see if your presolve could pick up this automatically.   But in the meantime, try this change in your code and see if it helps performance.

 

 

 

 

#CPLEXOptimizers
#DecisionOptimization

The issue you ran into with the 1217 error was reported here in the cvxpy github issue tracker and has already been addressed; it should be fixed in the next cvxpy release (whenever that is).

#CPLEXOptimizers
#DecisionOptimization

Originally posted by: anthzed

Top job guys.

Cvxpy works now and the speed improvements from the objective adjustment are fantastic.

#CPLEXOptimizers
#DecisionOptimization

Originally posted by: EdKlotz

I'm glad to hear you're all good now.   We'll have a look and see if CPLEX can pick up that speed improvement automatically, but in the meantime you know what to do.

#CPLEXOptimizers
#DecisionOptimization

Originally posted by: EdKlotz

I had a closer look at the two different formulations, and I now think I see why the alternate formulation in the attached DoCplex program (which could just as easily have been generated using the cvxpy API) it tighter.   I'll illustrate with the simplest  version of this model:  a single asset with long and short positions available, returns of 0 and risk coefficients of 1:

 

original model:

min (x +y)^2  s.t.
x*y = 0                                                // formulated via big M, SOS, indicators, whatever.
0 <= x <= 1                                        // long
-1 <= y <= 0                                      // short; note bounds
x - y = 1                          

longshort.py:

min x^2 + y^2  s.t.                               // only the objective differs.
x*y = 0
0 <= x <= 1                                        // long
-1 <= y <= 0                                      // short; note bounds
x - y = 1 

 

The initial relaxation removes the x*y = 0 constraint.  In the relaxation, we can set x = 1/2, y = -1/2.    In the first formulation, we get an objective of 0, which is
clearly the minimum possible objective.   By contrast, in the second formulation x = 1/2, y = -1/2. gives us an objective of 1/2, much tighter.   The difference is
that in the first objective the 2xy term in the relaxation subtracts off 1/2, taking the objective down from 1/2 to 0.   But if we remove that x*y term, we get the
tighter relaxation value.     I bet there are some other MIQPs out there that can be tightened in this way.  

#CPLEXOptimizers
#DecisionOptimization