I'm new to SOCP and want to try to get familiar with the format and how to solve it with cvxopt in python. However, for a simple toy example I'm struggling to get the right input format.

The problem I want to solve has the form

$$ \max_x c^T x$$ subject to $$ \| A_ix + b_i\|_2\le c_i^Tx + d_i, i=1,\dots, m$$ $$Fx = g $$

The example I came up with has the following parameters: $c=(0.02, 0.06)$. Additionally I'm given a symmetric positive semi-definite matrix $$\Sigma=\begin{bmatrix} 0.000025 & 0.000046\\ 0.000046& 0.024025 \end{bmatrix} $$ and have the constraints: $$ \sqrt{x^T\Sigma x} \le d_{max}=0.05$$ $$ 0\le x_i \le 1, i = 1, 2$$ $$ \sum_{i=1}^2 x_i = 1$$

I then first transformed the constrained into the proper SOCP format. In a next step I want to make them cvxopt acceptable.

1. Transform to SOCP The constraint $$ \sum_{i=1}^2 x_i = 1$$ is simple by taking $F=\begin{bmatrix} 1 & 1 \end{bmatrix}$ and $ g = 1$. Next, for the constraint $ \sqrt{x^T\Sigma x} \le d_{max}$ we use that for a positive semi-definite matrix $ \sqrt{x^T\Sigma x} = \|\Sigma^{\frac{1}{2}}x\|_2$ where $\Sigma^{\frac{1}{2}}$ denotes the Cholesky decomposition. This means, we have $A_1 = \Sigma^{\frac{1}{2}}$, $b_1=0$, $c_1 = 0$ and $d_1 = d_{max}$

The last constraints of the type $ 0\le x_i \le 1, i = 1, 2$ can be converted into a proper SOCP format by taking $A_2 = 0, b_2 =0, c_2=(-1, 0)$ and $d_2 = 1$. This gives $ x_1 \le 1$. To get $0 \le x_1$ we take $A_3 = 0, b_3 =0, c_3=(1, 0)$ and $d_3 = 0$. We get two additional constraints for $x_2$ in the same way with changing $c_4 = (0, -1)$ and $c_5=(0, 1)$.

2. Transform SCOP to cvxopt format

The cvxopt has slightly different format. The objective function and equality constraint is the same. However, the second order constraints have the form

$$ G_kx + s_k = h_k$$ $$s_{k0}\ge \|s_{k1}\|_2$$ where $s_k=(s_{k0}, s_{k1})$

The solver needs as input the list of the matrices $G_k$ and vectors $h_k$. If I'm not wrong I transfer each of the constraints $\| A_ix + b_i\|_2\le c_i^Tx + d_i$ to this format by taking

$$ G_i =\begin{bmatrix} -c_i \\ -A_i \end{bmatrix}$$ and $$ h_i = \begin{bmatrix} d_i \\ b_i \end{bmatrix}$$

Python toy model I've wanted to test this in python. Unfortunately, it doesn't work as expected and I don't see exactly what is wrong. I'm not sure if it is the implementation or one of the transformations above. Here is the python toy example

In [201]: import numpy as np

In [202]: import cvxopt as cpt

In [203]: Sigma = np.asarray([[0.000025, 0.000046],[0.00046, 0.024025]])

In [205]: np.linalg.eig(Sigma)
Out[205]: 
(array([2.41183657e-05, 2.40258816e-02]), array([[-0.99981638, -0.00191659],
        [ 0.01916244, -0.99999816]]))

As you can see the the matrix is really positive semi-definite.

In [37]: c = cpt.matrix(np.asarray([[0.02, 0.06]]).T)

In [38]: F = cpt.matrix(np.asarray([[1.0, 1.]]), tc='d')

In [39]: g = cpt.matrix(np.asarray([[1.0]]), tc='d')

In [40]: A0 = cpt.matrix(np.linalg.cholesky(Sigma))

In [41]: G1 = cpt.matrix(np.vstack(([[0, 0]],-A0)), tc='d')

In [42]: h1 = cpt.matrix(np.asarray([[0.05],[0.0], [0.0]]), tc='d')

In [43]: G2 = cpt.matrix(np.asarray([[-1, 0],[0, 0]]), tc='d')

In [44]: h2 = cpt.matrix(np.asarray([[1.0],[0.0]]), tc='d')

In [45]: G3 = cpt.matrix(np.asarray([[1, 0],[0, 0]]), tc='d')

In [46]: h3 = cpt.matrix(np.asarray([[0.0],[0.0]]), tc='d')

In [47]: G4 = cpt.matrix(np.asarray([[0, -1],[0, 0]]), tc='d')

In [48]: h4 = cpt.matrix(np.asarray([[1.0],[0.0]]), tc='d')

In [49]: G5 = cpt.matrix(np.asarray([[0, 1],[0, 0]]), tc='d')

In [50]: h5 = cpt.matrix(np.asarray([[0.0],[0.0]]), tc='d')

In [51]: Gq = [G1]

In [52]: [Gq.append(Gi) for Gi in [G2, G3, G4, G5]]

In [53]: hq = [h1]

In [54]: [hq.append(hi) for hi in [h2, h3, h4, h5]]

Gq and hq are lists which is required by cvxopt.solvers.scop. If I know call the solver I get:

In [55]: cpt.solvers.socp(-c, Gq=Gq, hq=hq, A=F, b=g)
     pcost       dcost       gap    pres   dres   k/t
 0: -3.9965e-02 -2.1105e+00  1e+01  2e+00  3e-16  1e+00
 1: -3.9369e-02 -1.3687e-02  5e+00  8e-01  1e-15  1e+00
 2: -3.3610e-02  1.0744e+02  1e+03  2e+00  2e-13  1e+02
 3: -3.3610e-02  1.0935e+04  1e+05  2e+00  4e-12  1e+04
 4: -3.3610e-02  1.0936e+06  1e+07  2e+00  1e-09  1e+06
Certificate of primal infeasibility found.
Out[55]: 
{'dual infeasibility': None,
 'dual objective': 1.0,
 'dual slack': 0.9204046569064762,
 'gap': None,
 'iterations': 4,
 'primal infeasibility': None,
 'primal objective': None,
 'primal slack': None,
 'relative gap': None,
 'residual as dual infeasibility certificate': None,
 'residual as primal infeasibility certificate': 5.782986313941418e-08,
 'sl': None,
 'sq': None,
 'status': 'primal infeasible',
 'x': None,
 'y': <1x1 matrix, tc='d'>,
 'zl': <0x1 matrix, tc='d'>,
 'zq': [<3x1 matrix, tc='d'>,
  <2x1 matrix, tc='d'>,
  <2x1 matrix, tc='d'>,
  <2x1 matrix, tc='d'>,
  <2x1 matrix, tc='d'>]}

I don't understand why. Even if I use $d_{max}=0.2$ we should get a solution of $x_1 = 0, x_2=1$. Any help would be much appreciated. Please note that I want to use cvxopt and I'm thankful but not interested in answers / comments suggesting to use another module like cvxpy. I have a restriction for real cases I want to work on, after finishing some of the basic examples, to use cvxopt.

up vote 5 down vote accepted
+50

Instead of direclty answering your question, I will propose to use a, so called, modeling language for optimization problems, which allows to formulate your problem in a natural way, and have the language translate it to the appropriate format required by solvers such as CVXOPT. The code below solves your toy problem using the cvxpy package, which is compatible with many solvers, including CVXOPT.

import numpy as np
import cvxpy
from scipy.linalg import sqrtm  # for finding the squared root of Sigma

x = cvxpy.Variable(2) # optimization vector variable
c = cvxpy.Parameter(2) # placeholder for vector c
A = cvxpy.Parameter((2,2))  # placeholder for sqrtm(Sigma)
dmax = cvxpy.Parameter() # placeholder for parameter dmax

obj = cvxpy.Maximize(cvxpy.sum(x*c))  #define objective function 
cons = [cvxpy.norm(A@x,2)<=dmax,      # define set of constraints
        cvxpy.sum(x)==1, 
        x[0]>=0, x[1]>=0,
        x[0]<=1, x[1]<=1]

prob = cvxpy.Problem(obj, cons)  #setup the problem

# instantiate the problem for a specific value of parameters
A.value = sqrtm(np.asarray([[0.000025, 0.000046], [0.000046, 0.024025]]))
c.value = np.asarray([0.02, 0.06])
dmax.value = 0.05

prob.solve()  # solve the problem

x_opt = x.value # the optimal variable

print(x_opt, np.linalg.norm(A.value @ x_opt), np.sum(x_opt) )

[0.6795 0.3205] 0.049999998463030305 0.9999999999999976

Note how simple and natural is the problem definition. No need to know the technicalities of the solver(s). Here, cvxopt chooses automatically the "best" solver (out of those available in your system and supported by cvxpy) for the problem. You may explicitly indicate the solver (e.g., CVXOPT), if needed.

EDIT: Here's how you can automatically generate the domain constraints instead of manually writing them (helpful when there are many of them and/or you don't know dimensions beforehand)

# define the `cons` python list with cvxpy constraints you can specify before runtime (it may as well be cons=[])
cons = [cvxpy.norm(A@x,2)<=dmax, cvxpy.sum(x)==1]

# generate domain constraints on the fly:
# N is the dimension, possibly available only at runtime; for the toy example, N = 2
for i in range(N):
     cons = cons + [x[i]>=0, x[i]<=1]
  • many thanks for your answer. I do see the simplicity of the solver. Two comments 1: As I stated at the very end. I have some restrictions on using cvxopt directly. Without knowing cvxpy in detail it seems much more cumbersome to write the problem if the dimension of your variables are not know until runtime. 2: My background is pure Mathematics. I have rather less experience in more applied math. That is why I like to understand the detail of problem, what kind of problem it is, what methods you need to solve etc :) – math Aug 10 at 6:11
  • @math I am sorry, I did not see your cvxopt restriction. Knowing the dimensions only at runtime is no big of a problem for cvxpy; you just define variables and/or parameters on the fly. – Stelios Aug 10 at 7:41
  • No worries....I've upvoted your answer but I wait to see if someone can help me with cvxopt. Do you have a reference for the dimension at runtime? Because in a sense you need to create expressions of the form x[0]>= 0 and so on for a not known dimension. I do see that this can be done, but creating a varying matrix in cvxopt is at least IMHO much easier. I would be interested to see an easy way. If this is off topic, I'm sorry – math Aug 10 at 8:24
  • @math Thanks for the upvote :) I have edited my answer with an example on how you would generate automatically constraint specifications at runtime. – Stelios Aug 10 at 10:10
  • many thanks! This is indeed much simpler than what I've expected. – math Aug 10 at 10:12

Your Answer

 

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.