cvxpy / cvxpy

A Python-embedded modeling language for convex optimization problems.
https://www.cvxpy.org
Apache License 2.0
5.44k stars 1.07k forks source link

How to construct a SOCP problem and solve #1185

Closed ghost closed 3 years ago

ghost commented 3 years ago

Hi Im looking at this SOCP problem (problem 11) (here is the scihub link in case you can't access), and here is a snippet of the problem:

Screenshot from 2020-12-09 15-10-30

I've looked at this example, but is still stuck and can't get the problem to be solved. I keep getting

Problem does not follow DCP rules

on the objective and constraints.

Here is my attempt:

import cvxpy as cp
import numpy as np

n_assets = 5

covmat_sqrt = cp.Parameter((n_assets, n_assets))

x = cp.Variable(n_assets)

zeta = []
zero = []
for n in range(n_assets):
    zeta.append((covmat_sqrt @ x)[n])
    zero.append(cp.Constant(0))

psi = cp.sqrt(cp.sum_squares(covmat_sqrt @ x)/n_assets)
gamma = cp.sqrt(cp.sum_squares(covmat_sqrt @ x))

objective = cp.Minimize(psi)
constraint_soc1 = [cp.SOC(zeta[n], (covmat_sqrt @ x)[n]) for n in range(n_assets)]
constraint_soc2 = [cp.SOC(x[n] * zeta[n], gamma ** 2) for n in range(n_assets)]
constraint_soc3 = [cp.SOC(zeta[n], zero[n]) for n in range(n_assets)]
constraint_soc4 = [cp.SOC(x[n], zero[n]) for n in range(n_assets)]
constraint = [cp.sum_squares(covmat_sqrt @ x) <= n_assets * (psi ** 2),
              cp.sum(x) == 1, psi >= 0, gamma >= 0]

constraint_all = constraint + constraint_soc1 + constraint_soc2 + constraint_soc3 + constraint_soc4

matrix = np.random.random((n_assets, n_assets))
covmat_sqrt.value = matrix.T @ matrix # to make positive definite

prob = cp.Problem(objective, constraint_all)

prob.solve()
print("status:", prob.status)
print("optimal value", prob.value)
rileyjmurray commented 3 years ago

We usually refer people to Stackoverflow for these kinds of questions. However, your case is easy to identify. The issue is that t**2 <= x*y is a "rotated second-order cone" constraint, which CVXPY does not automatically identify. For this type of constraint, you should use cp.quad_over_lin(t,y) <= x. The behavior of quad_over_lin is described here: https://www.cvxpy.org/tutorial/functions/index.html.

rileyjmurray commented 3 years ago

Follow-up: it seems your constraint x.T @ Q @ x <= N * p**2 will also cause a problem. This can also be reformulated using the quad_over_lin function. If you have more questions, please post them to Stackoverflow with a CVXPY tag.

ghost commented 3 years ago

Hi, I've now only noticed that the sample code I gave does not match the equation symbols, sorry about that.

Anyway I've posted the question on Stack Overflow. I've taken your comments into consideration, but is unfortunately still stuck. I've also added more details about me needing to use cvxpylayers as well. Any help is appreciated, thanks !

ghost commented 3 years ago

hi @rileyjmurray I've updated my question again as I realized its probably a multivariate problem.

rileyjmurray commented 3 years ago

@blenderben2 it looks like you've got some good advice from Erling since you made the comment above. If you remain stuck, you can ping me here (again) and I will chime in on that thread.

ghost commented 3 years ago

hi @rileyjmurray I am very new to convex optimization in general. Ideally, I would like to solve expression 8 (as in the snippet) within cvxpy, however I am a bit confused as to what is my current issue and what next steps to take.

From my point of view, expression 8 is sufficient in terms of the maths to be a SOCP problem, no alteration required (I might be wrong). So the suggestion from Erling about the cookbook looks like just SOCP formulation in terms of math, which I presumably already have (expression 8).

Although I am getting an ECOS unable to be solved issue, I think this is due to how I write the math SOCP part of the constraints (the ones with i = 1, ...., N) in cvxpy that is wrong. Why I say this? Because I am referencing this SOCP example by cvxpy documentation and I see there is a loop section which mainly involves the Parameter constants/data (A, b, c, d) SOCP constraints,

for i in range(m):
    A.append(np.random.randn(n_i, n))
    b.append(np.random.randn(n_i))
    c.append(np.random.randn(n))
    d.append(np.linalg.norm(A[i] @ x0 + b, 2) - c[i].T @ x0)

while my SOCP constraints involve Variables. I have tried to put these Variables in a loop without success.

My question is, is there a way to do everything in cvxpy while using the default solver (ECOS)? I want to rule out that my interpretation of the expression 8's math is at least correctly written in a SOCP way in cvxpy code before considering other solvers.

Thank you for your patience.

rileyjmurray commented 3 years ago

When getting this problem to work with CVXPY, you should not use Expression 8, because Expression 8 is very much NOT allowed in SOCP modeling. You should use the formulation in display (11), and ONLY that formulation. Based on the reformulation pointed out by Erling you can (and it seems you have) written the problem in a way that CVXPY accepts. From here you've found the ECOS failed due to numerical issues, so you need to try another solver. In this case you can use SCS, which is also a default solver in CVXPY.

I don't really have anything else to offer at this point, but best of luck in getting this to work for your situation.

ghost commented 3 years ago

Hi @rileyjmurray I have made a big mistake in my previous comment. By expression 8 I mean 11 which is what is written in EDIT 2 and 3 of the stack overflow question. I'm very sorry for this confusion.

I will take your word that my formulation 11 written in cvxpy is indeed correct. I will try different solvers as you suggestion. Thank you for your time.

ghost commented 3 years ago

Based on the reformulation pointed out by Erling you can (and it seems you have) written the problem in a way that CVXPY accepts

One question is that how does cvxpy know that my problem is SOCP without specifying the cp.SOC

rileyjmurray commented 3 years ago

The purpose of CVXPY is to handle converting appropriate nonlinear programs (such as those involving quad_over_lin) into suitable standard forms (such as SOCP's).

ghost commented 3 years ago

I was able to solve using XPRESS. Thank you for your time.