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):
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.
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.