How to transform this SOCP to the format required by cvxopt

Question

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.

Answer 1

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]