Search code examples
portfoliocvxpyconvex-optimizationcvxopt

How to construct a SOCP problem and solve using cvxpy and cvxpylayers


I'm trying to solve a SOCP problem using cvxpy and integrating it to cvxpylayers. I'm 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 (note min (p-t) comes from an adaptation of problem 4 using expression 8 in the link above):

enter image description here

Go to ---> EDIT 2

OLD

I've looked at this example, but is still stuck and can't get the problem to be solved. Here is a sample code at my attempt:

import cvxpy as cp
import numpy as np

N = 5

Q_sqrt = cp.Parameter((N, N))
x = cp.Variable(N)

z = []
zero = []
for n in range(N):
    z.append((Q_sqrt @ x)[n])
    zero.append(cp.Constant(0))
    
p = cp.sqrt(cp.sum_squares(Q_sqrt @ x)/N) # based on expression 8
t = cp.sqrt(cp.min(x * (Q_sqrt @ x)))     # based on expression 8

objective = cp.Minimize(p - t)            # based on expression 8

constraint_soc1 = [cp.SOC(z[n], (Q_sqrt @ x)[n]) for n in range(N)]
constraint_soc2 = [cp.SOC(cp.quad_over_lin(t, z[n]), x[n]) for n in range(N)]
constraint_soc3 = [cp.SOC(z[n], zero[n]) for n in range(N)]
constraint_soc4 = [cp.SOC(x[n], zero[n]) for n in range(N)]

constraint_other = [cp.sum_squares(Q_sqrt @ x) <= N * (p ** 2),
                    cp.sum(x) == 1, p >= 0, t >= 0]

constraint_all = constraint_other + constraint_soc1 + constraint_soc2 + constraint_soc3 + constraint_soc4

matrix = np.random.random((N, N))
Q_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)

Note that I am using cp.sum_squares(Q_sqrt @ x) syntax to ensure it is not only DCP-compliant but DPP-compliant as well (refer Section 4.1).

Firstly, I noticed my p = cp.sqrt(cp.sum_squares(Q_sqrt @ x)/N) part of the objective (first term of expression 8) is not DCP and is not sure how to fix it.

Also, I was told that this x.T @ Q @ x <= N * p**2 constraint was required to be reformulated using cp.quad_over_lin, however, as stated I require it to follow DPP rules (for cvxpylayers) as well, hence I am unsure how to change cp.sum_squares(Q_sqrt @ x) <= N * (p ** 2) to follow cp.quad_over_lin.

Lastly, I am sure I have other parts of my code wrongly implemented as well. Any help on this problem is really appreciated!

EDIT 1: note that the final solution does not have to call cvxpylayers, but just be DPP-compliant so that it can be integrated with cvxpylayers after.

__________________________________________________________________________

NEW

EDIT 2: I have now realized that this is probably a multivariate problem, to minimize p - t given variables x,z,p,t. I don't know why this didn't cross my mind earlier.

Given the snippet problem above, my latest attempt is as follows:

import cvxpy as cp
import numpy as np

N = 5

Q_sqrt = cp.Parameter((N, N))
Q = cp.Parameter((N, N))

x = cp.Variable(N)
z = cp.Variable(N)
p = cp.Variable()
t = cp.Variable()

objective = cp.Minimize(p - t)

constraint_soc = [z == Q @ x, x.value * z >= t ** 2, z >= 0, x >= 0]

constraint_other = [cp.quad_over_lin(Q_sqrt @ x, N) <= p ** 2, cp.sum(x) == 1, p >= 0, t >= 0]

constraint_all = constraint_other + constraint_soc

matrix = np.random.random((N, N))
a_matrix = matrix.T @ matrix
Q.value = a_matrix
Q_sqrt.value = np.sqrt(a_matrix)

prob = cp.Problem(objective, constraint_all)

prob.solve(verbose=True)
print("status:", prob.status)
print("optimal value", prob.value)

However, I am still getting the error:

DCPError: Problem does not follow DCP rules. Specifically: The following constraints are not DCP: quad_over_lin(param7788 @ var7790, 5.0) <= power(var7792, 2.0) , because the following subexpressions are not: |-- quad_over_lin(param7788 @ var7790, 5.0) <= power(var7792, 2.0)

for this x.T @ Q @ x <= N * p**2 constraint (my attempt: cp.quad_over_lin(Q_sqrt @ x, N) <= p ** 2).

Besides that the problem is also SOCP as mentioned above, so my constraints in constraint_soc is probably wrong, but I am not sure how to loop/assign variables based on this.

EDIT 3: after ErlingMOSEK's suggestion on the constraint, I've changed my constraint from cp.quad_over_lin(Q_sqrt @ x, N) <= p ** 2 to cp.quad_over_lin(Q_sqrt @ x, p) <= N * p and it now runs, however,

A new error appeared

Solver 'ECOS' failed. Try another solver, or solve with verbose=True for more information.

which I assume is due to the problem being SOCP (my second issue) and I am unsure how to construct the SOCP just based on this example.

EDIT 4: using EDIT 2/3 and Erling's answer along with a different solver as Erling suggested, the problem was solved. Note I used XPRESS solver instead.


Solution

  • You should change

    x.T @ Q @ x <= N * p**2

    to

    (x.T @ Q @ x)/p <= N * p

    assuming p>=0.

    Btw if you want to know more about how to formulate this as a SOCP, then consult the Mosek modelling cookbok.