CVXOPT solvers.lp faster without the option solver = 'glpk' (??)

[guillaume lecue]

I'm using CVXOPT with python 3 on MAC to solve an \ell_1-minimization procedure (here it is for the Basis Pursuit procedure).

According to <a \href = 'https://scaron.info/blog/linear-programming-in-python-with-cvxopt.html'> Stéphane Caron's webpage </a>, using the solver = 'glpk' option should be much faster. 

So I did the following installation

1. Install glpk via brew :

> brew install glpk

2. Edit the setup.py of CVXOPT:

 First thing set variable with BUILD_GLPK to 1: 

> # Set to 1 if you are installing the glpk module.

>BUILD_GLPK = 1 

Then indicate the path to *libglpk*

># Directory containing libglpk (used only when BUILD_GLPK = 1).

>GLPK_LIB_DIR = '/usr/local/Cellar/glpk/4.52/lib'

Then indicate the path to *glpk.h*

># Directory containing glpk.h (used only when BUILD_GLPK = 1).

>GLPK_INC_DIR = '/usr/local/Cellar/glpk/4.52/include'

3. Go to the CVXOPT's 'setup.py' folder and run

> python setup.py install

So everything went well. I can use solvers.lp(c, G, h, A, b, solver = 'glpk') with the solver = 'glpk' option BUT my problem is that:

*** It is much slower with the solver = 'glpk' option than with no option.

Here is the result  I get:

%timeit solvers.lp(c, G, h, A, b)

1 loops, best of 3: 296 ms per loop

%timeit solvers.lp(c, G, h, A, b, solver = 'glpk')

1 loops, best of 3: 2.08 s per loop

I tried several cases and always solver = 'glpk' was much slower.

Do you have any idea why it is much slower with glpk than with the build-in CVX solver ?

Thanks for your help,
Guillaume.

I'm using the following matrices:

N, m, s = 1000, 200, 30
A, x = measurement_matrix(m, N), signal(N, s) 
y = np.dot(A,x)
c, G, h, A, b = cvx_mat(A, y)

def cvx_mat(A, y):
    '''A, y: numpy array or cvx matrices
    return parameters c, G, h, A, b for CVXOPT solvers.lp method'''
    A = matrix(A)
    y = matrix(y)
    m, N = matrix(A).size
    # matrix c, G, h, A, b
    c = matrix(np.ones(2*N))
    G = -spdiag([1]*2*N)
    h = matrix(np.zeros(2*N))
    A = matrix([[A],[-A]])
    b = y
    return c, G, h, A, b

def measurement_matrix(m, N, rand_type = "gauss", power = 2):
    if rand_type == "student":
        A = np.random.standard_t(power, size=(m, N))
    else:
        A = np.random.randn(m, N)
        A = np.multiply(np.power(np.absolute(A),power), np.sign(A))          
    return A 

def signal(N, s):
    """Construct a signal of lenght n with s gaussian variables localized randomly and zeros everywhere else"""
    x = np.zeros(N)
    sel = np.random.permutation(N)[0:s]
    x[sel] = np.random.randn(s)
    return x