Simple example not converging, while MATLAB does.

[jah]

Hi, I'm having some trouble getting a simple problem to converge.  I have a working version using MATLAB as follows:

    A = [1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.;

         1.,  1.,  0.,  0.,  0.,  0.,  0.,  0.;

         0.,  0.,  1.,  1.,  0.,  0.,  0.,  0.;

         0.,  0.,  0.,  0.,  1.,  1.,  0.,  0.;

         0.,  0.,  0.,  0.,  0.,  0.,  1.,  1.;

         1.,  0.,  1.,  0.,  0.,  0.,  0.,  0.;

         0.,  1.,  0.,  1.,  0.,  0.,  0.,  0.;

         0.,  0.,  0.,  0.,  1.,  0.,  1.,  0.;

         0.,  0.,  0.,  0.,  0.,  1.,  0.,  1.;

         1.,  0.,  0.,  0.,  1.,  0.,  0.,  0.;

         0.,  1.,  0.,  0.,  0.,  1.,  0.,  0.;

         0.,  0.,  1.,  0.,  0.,  0.,  1.,  0.;

         0.,  0.,  0.,  1.,  0.,  0.,  0.,  1.]

    b = [1; 1/4; 1/4; 1/4; 1/4; 1/2; 0; 1/4; 1/4; 1/2; 0; 1/4; 1/4]

    [m,n] = size(A)

    cvx_begin

    variable x(n) nonnegative

    minimize (norm(A*x - b, inf))

    cvx_end

If you run this, you find that:  x = [1/4; 0; 1/4; 0; 1/4; 0; 0; 1/4]

Now I want to implement this using cvxopt. My code is as follows:

import numdifftools

from cvxopt import matrix

from cvxopt.solvers import cp, options

import numpy as np

import scipy.linalg as linalg

def full_rank(A, b):

    """

    From a linear system Ax = b, return Bx = c such that B has full rank.

    """

    U, S, Vh = linalg.svd(A)

    Smat = linalg.diagsvd(S, A.shape[0], A.shape[1])

    # See np.linalg.matrix_rank

    tol = S.max() * max(A.shape) * np.finfo(S.dtype).eps

    rank = np.sum(S > tol)

    B = np.dot(Smat, Vh)[:rank]

    c = np.dot(U.transpose(), b)[:rank]

    return B, c, rank

A = np.array([

    [ 1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.],

    [ 1.,  1.,  0.,  0.,  0.,  0.,  0.,  0.],

    [ 0.,  0.,  1.,  1.,  0.,  0.,  0.,  0.],

    [ 0.,  0.,  0.,  0.,  1.,  1.,  0.,  0.],

    [ 0.,  0.,  0.,  0.,  0.,  0.,  1.,  1.],

    [ 1.,  0.,  1.,  0.,  0.,  0.,  0.,  0.],

    [ 0.,  1.,  0.,  1.,  0.,  0.,  0.,  0.],

    [ 0.,  0.,  0.,  0.,  1.,  0.,  1.,  0.],

    [ 0.,  0.,  0.,  0.,  0.,  1.,  0.,  1.],

    [ 1.,  0.,  0.,  0.,  1.,  0.,  0.,  0.],

    [ 0.,  1.,  0.,  0.,  0.,  1.,  0.,  0.],

    [ 0.,  0.,  1.,  0.,  0.,  0.,  1.,  0.],

    [ 0.,  0.,  0.,  1.,  0.,  0.,  0.,  1.]

])

b = np.array([ 1.  ,  0.25,  0.25,  0.25,  0.25,  0.5 ,  0.,

               0.25,  0.25,  0.5 ,  0.  ,  0.25,  0.25])

# Dimension of optimization variable

n = 8

# Number of nonlinear constraints

m = 0

# Nonnegativity constraint

G = matrix( -1 * np.eye(n) )   # G should have shape: (K,n) = (n,n)

h = matrix( np.zeros((n,1)) )  # h should have shape: (K,1) = (n,1)

# Linear equality constraints (these are not independent constraints)

AA, bb, rank = full_rank(A, b)

A = matrix(AA)

b = matrix(bb)  # now a column vector

def func(x):

    return np.max(np.abs(A * x - b))

grad = numdifftools.Gradient(func)

hess = numdifftools.Hessian(func)

def F(x=None, z=None):

    # x has shape: (n,1)

    # z has shape: (m+1,1) and is the Hessian of f_0

    if x is None and z is None:

        p = np.ones(n) / n

        return (m, matrix(p))

    xarr = np.array(x)[:,0]

    if np.any(xarr > 1) or np.any(xarr < 0):

        return None

    if not np.allclose(np.sum(xarr), 1):

        return None

    f = func(xarr)

    Df = grad(xarr)

    Df = matrix(Df.reshape((1,n)))

    if z is None:

        return (f, Df)

    H = hess(xarr)

    H = matrix(H)

    return (f, Df, H)

result = cp(F=F, G=G, h=h, dims={'l':n, 'q':[], 's':[]}, A=A, b=b)

out = np.asarray(result['x'])

print out

Let me know if I've implemented this incorrectly, but I think it's correct (if a bit verbose). However when this runs, it does not converge. The feasible set consists of a single point which is the answer that MATLAB gives. Here is the log:

     pcost       dcost       gap    pres   dres

 0:  0.0000e+00  1.5752e-01  9e+00  1e+00  1e+00

 1:  6.1745e-01  2.9734e-01  5e+00  6e-01  4e-01

 2:  9.7355e-01  4.3862e-01  3e+00  3e-01  3e-01

 3:  1.0115e+00  4.5308e-01  3e+00  3e-01  2e-01

 4:  1.0677e+00  4.6902e-01  3e+00  2e-01  2e-01

 5:  1.0957e+00  4.6693e-01  3e+00  2e-01  2e-01

 6:  1.1030e+00  4.6411e-01  3e+00  2e-01  2e-01

 7:  1.1049e+00  4.6324e-01  3e+00  2e-01  2e-01

 8:  1.1054e+00  4.6301e-01  3e+00  2e-01  2e-01

 9:  1.1054e+00  4.6300e-01  3e+00  2e-01  2e-01

10:  1.1054e+00  4.6299e-01  3e+00  2e-01  2e-01

11:  1.1054e+00  4.6299e-01  3e+00  2e-01  2e-01

What should I be doing to get this working with cvxopt?

[jah]

On Tuesday, March 17, 2015 at 12:03:59 AM UTC-5, jah wrote:

 

result = cp(F=F, G=G, h=h, dims={'l':n, 'q':[], 's':[]}, A=A, b=b)

Sorry this line should be:

result = cp(F=F, G=G, h=h, dims={'l':n, 'q':[], 's':[]})

 

without any linear equality constraints. Even so, it still fails to converge:

     pcost       dcost       gap    pres   dres
 0:  0.0000e+00  1.5752e-01  9e+00  1e+00  1e+00

 1: -4.4086e-06  1.5751e-01  9e+00  1e+00  1e+00
 2: -5.5092e-06  1.5751e-01  9e+00  1e+00  1e+00
 3: -5.7843e-06  1.5751e-01  9e+00  1e+00  1e+00
 4: -5.8530e-06  1.5751e-01  9e+00  1e+00  1e+00
 5: -5.8874e-06  1.5751e-01  9e+00  1e+00  1e+00
 6: -5.8877e-06  1.5751e-01  9e+00  1e+00  1e+00
 7: -5.8878e-06  1.5751e-01  9e+00  1e+00  1e+00
 8: -5.8879e-06  1.5751e-01  9e+00  1e+00  1e+00
 9: -4.1058e-15  1.5752e-01  9e+00  1e+00  1e+00
10: -6.1587e-15  1.5752e-01  9e+00  1e+00  1e+00
11: -6.1908e-15  1.5752e-01  9e+00  1e+00  1e+00
12: -6.2228e-15  1.5752e-01  9e+00  1e+00  1e+00
13: -6.2549e-15  1.5752e-01  9e+00  1e+00  1e+00
14: -6.2870e-15  1.5752e-01  9e+00  1e+00  1e+00

15: -6.3191e-15  1.5752e-01  9e+00  1e+00  1e+00

[Martin]

Hi,

You can express this problem as a linear program (see. eg. Boyd & Vandenberghe's "Convex Optimization", Section 6.1), and hence you can use solvers.conelp() or solvers.lp(). Note also that there is a Python modeling package called CVXPY that is similar to CVX for MATLAB.

Martin

[jah]

On Tuesday, March 17, 2015 at 3:18:08 AM UTC-5, Martin wrote:

Hi,

You can express this problem as a linear program (see. eg. Boyd & Vandenberghe's "Convex Optimization", Section 6.1), and hence you can use solvers.conelp() or solvers.lp(). Note also that there is a Python modeling package called CVXPY that is similar to CVX for MATLAB.

Hi Martin, 

Thanks for the tip. If I use lp directly (or use cvx.modeling), cvxopt complains about a rank issue. So I'm still not able to get "equivalent" code to work with cvxopt. Here is what I did:

    x = variable(n)

    t = variable()

    AA = matrix(A)

    bb = matrix(b)

    c1 = (-t <= AA * x - bb)

    c2 = ( AA * x - bb <= t)

    op(t, [c1, c2]).solve()

Any other help you can provide would be greatly appreciated.  Here is the complete code demonstrating two different ways of formulating the problem, both of which fail:

    http://codepad.org/X5gP7bBe

I've also pasted the complete code at the end of this email.

------------------

"""
Equivalent MATLAB code:

    A = [1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.;
         1.,  1.,  0.,  0.,  0.,  0.,  0.,  0.;
         0.,  0.,  1.,  1.,  0.,  0.,  0.,  0.;
         0.,  0.,  0.,  0.,  1.,  1.,  0.,  0.;
         0.,  0.,  0.,  0.,  0.,  0.,  1.,  1.;
         1.,  0.,  1.,  0.,  0.,  0.,  0.,  0.;
         0.,  1.,  0.,  1.,  0.,  0.,  0.,  0.;
         0.,  0.,  0.,  0.,  1.,  0.,  1.,  0.;
         0.,  0.,  0.,  0.,  0.,  1.,  0.,  1.;
         1.,  0.,  0.,  0.,  1.,  0.,  0.,  0.;
         0.,  1.,  0.,  0.,  0.,  1.,  0.,  0.;
         0.,  0.,  1.,  0.,  0.,  0.,  1.,  0.;
         0.,  0.,  0.,  1.,  0.,  0.,  0.,  1.]
    b = [1; 1/4; 1/4; 1/4; 1/4; 1/2; 0; 1/4; 1/4; 1/2; 0; 1/4; 1/4]
    [m,n] = size(A)
    cvx_begin
    variable x(n) nonnegative
    minimize (norm(A*x - b, inf))
    cvx_end

Answer:

    x = [1/4; 0; 1/4; 0; 1/4; 0; 0; 1/4]

"""

from cvxopt import matrix
from cvxopt.solvers import lp
from cvxopt.modeling import variable, op

def method1(fullrank=True):
    # ValueError: Rank(A) < p or Rank([G; A]) < n

    x = variable(n)
    t = variable()

    if fullrank:

        AA, bb, rank = full_rank(A, b)

    else:
        AA, bb = A, b

    AA = matrix(A)
    bb = matrix(b)

    c1 = (-t <= AA * x - bb)
    c2 = ( AA * x - bb <= t)

    op(t, [c1, c2]).solve()


def method2():
    # ValueError: Rank(A) < p or Rank([G; A]) < n

    # min c^T z     :   G z <= h

    # z = [x; t]

    # c = [0; 1]
    c = np.zeros(n + 1)
    c[-1] = 1
    c = matrix(c)

    # G = [A, -1; -A -1]
    col = (A.shape[0], 1)
    negone = -np.ones(col)
    Gone = np.hstack([A, negone])
    Gtwo = np.hstack([-A, negone])
    G = np.vstack([Gone, Gtwo])
    G = matrix(G)

    # h = [b; -b]
    h = np.hstack([b, -b])
    h = matrix(h)

    result = lp(c, G, h)

    out = np.asarray(result['x'])
    print out

method1()

[Martin]

I obtain the same solution if I add the constraint x >= 0 to your "method 1".

    c1 = (-t <= AA * x - bb)
    c2 = ( AA * x - bb <= t)

    c3 = ( x >= 0 )

    op(t, [c1, c2, c3]).solve()

[jah]

Sigh...clearly, I'm still learning.

Much appreciated.

[jah]

On Tuesday, March 17, 2015 at 1:31:33 PM UTC-5, Martin wrote:

I obtain the same solution if I add the constraint x >= 0 to your "method 1".

    c1 = (-t <= AA * x - bb)
    c2 = ( AA * x - bb <= t)
    c3 = ( x >= 0 )

    op(t, [c1, c2, c3]).solve()

I just wanted to mention that if you do the following in MATLAB:

    cvx_begin

    variable x(n)

    minimize (norm(A*x - b, inf))

    cvx_end

(that is, without the nonnegative constraint), then MATLAB gives: 

0.357660907871685

-0.107660907871686

0.142339092128315

0.107660907871686

0.142339092128315

0.107660907871686

0.107660907871685

0.142339092128315

whereas CVXOPT raises the exception I mentioned previously. Not sure if that is expected or not.

[Martin]

You're not doing the rank reduction correctly. The matrix A in your example does not have full column rank, and hence adding a vector from null(A) to x does not change the objective value. This basically means that you can reduce the number of variables in your problem. 

To see this, let x = V1*y + V2*z where V = [V1 V2] is the right singular vectors of A. It immediately follows that A*x = U1*S1*y where the matrix U = [U1 U2] denotes the left singular vectors and S1 is a diagonal matrix with the nonzero singular values. Thus, we can replace A*x with U1*S1*y in the original problem which reduces the number of variables to the dimension of y which is equal to the (column) rank of A.

After solving for y, we obtain the solutions to the original problem parameterized by z as x = V1*y + V2*z. The solution returned by CVX in MATLAB is one of many solutions.

Hope this helps.