From CVX to CVXPY or CVXOPT

silgon picture silgon · Jun 4, 2015 · Viewed 9.9k times · Source

I've been trying to pass some code from Matlab to Python. I have the same convex optimization problem working on Matlab but I'm having problems passing it to either CVXPY or CVXOPT.

n = 1000;
i = 20;
y = rand(n,1);
A = rand(n,i);
cvx_begin
variable x(n);
variable lambda(i);
minimize(sum_square(x-y));
subject to
    x == A*lambda;
    lambda >= zeros(i,1);
    lambda'*ones(i,1) == 1;
cvx_end

This is what I tried with Python and CVXPY.

import numpy as np
from cvxpy import *

# Problem data.
n = 100
i = 20
np.random.seed(1)
y = np.random.randn(n)
A = np.random.randn(n, i)

# Construct the problem.
x = Variable(n)
lmbd = Variable(i)
objective = Minimize(sum_squares(x - y))
constraints = [x == np.dot(A, lmbd),
               lmbd <= np.zeros(itr),
               np.sum(lmbd) == 1]

prob = Problem(objective, constraints)

print("status:", prob.status)
print("optimal value", prob.value)

Nonetheless, it's not working. Does any of you have any idea how to make it work? I'm pretty sure my problem is in the constraints. And also it would be nice to have it with CVXOPT.

Answer

silgon picture silgon · Jun 4, 2015

I think I got it, I had one of the constraints wrong =), I added a random seed number in order to compare the results and check that are in fact the same in both languages. I leave the data here so maybe this is useful for somebody someday ;)

Matlab

rand('twister', 0);
n = 100;
i = 20;
y = rand(n,1);
A = rand(n,i);
cvx_begin
variable x(n);
variable lmbd(i);
minimize(sum_square(x-y));
subject to
    x == A*lmbd;
    lmbd >= zeros(i,1);
    lmbd'*ones(i,1) == 1;
cvx_end

CVXPY

import numpy as np
from cvxpy import *

# random seed
np.random.seed(0)

# Problem data.
n = 100
i = 20
y = np.random.rand(n)
# A = np.random.rand(n, i)  # normal
A = np.random.rand(i, n).T  # in this order to test random numbers

# Construct the problem.
x = Variable(n)
lmbd = Variable(i)
objective = Minimize(sum_squares(x - y))
constraints = [x == A*lmbd,
               lmbd >= np.zeros(i),
               sum_entries(lmbd) == 1]

prob = Problem(objective, constraints)
result = prob.solve(verbose=True)

CVXOPT is pending.....