Created by David Stonestrom 6/4/2014.
The purpose of this document is to demonstrate debugging CVXPY DCP errors resulting from constraints.
Let's look at an optimization problem with a simple objective and two constraints.
$$ \textbf{maximize } f\left(x\right) = \sum_{i=1}^{n}{x_i}$$
where $x \in \mathbb{R}^n$ is the variable and $A \in \mathbb{R}^{n \times m}_{+}$, $b \in \mathbb{R}^m$, and $c \in \mathbb{R}$ are the problem data.
In [1]:
import cvxpy as cvx
import numpy
m = 3; # number of inequality constraints
n = 6; # number of variables
# generate some fake data
numpy.random.seed(0) # for repeatibility
A = numpy.random.random([m,n])
b = numpy.random.randn(m,1) + 0.5
c = 100 * numpy.random.random(1)
#objective
x = cvx.Variable(n,1)
objectiveExpression = cvx.sum(x)
# writing the constraints individually with names will make debugging them easier
constraint_1 = A*x >= b
constraint_2 = cvx.sum(cvx.exp(x)) == c
# try the problem as written
objective = cvx.Maximize(objectiveExpression)
constraints = [constraint_1, constraint_2]
problem = cvx.Problem(objective, constraints)
# Try solving the problem. Print error if any.
try:
problem.solve()
except Exception as e:
print e
All right, lets look at the three components of the problem.
In [2]:
print "Objective: ", objective.is_dcp()
print "Constraint 1: ", constraint_1.is_dcp()
print "Constraint 2: ", constraint_2.is_dcp()
Okay, we can dig into constraint 2 further by examining the curvature on each side.
In [3]:
print "Constraint 2 curvatures:"
print constraint_2.lh_exp.curvature
print constraint_2.rh_exp.curvature
Since the second constraint is an equality constraint, it needs each side to be AFFINE or CONSTANT in order to be DCP compliant.
We can make the equality constraint DCP compliant with a change of variables $y_i = e^{x_i}$. Note that $x_i = ln\left(y_i\right)$. The optimization problem is now:
$$ \textbf{maximize } g\left(y\right) = \sum_{i=1}^{n}{log\left(y\right)}$$$$\text{subject to } A\log\left(y\right) \succeq b$$$$\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\sum_{i = 1}^{n}{y} = c$$Lets look at how this transforms the constraints.
In [4]:
y = cvx.Variable(n,1) # y = exp(x)
constraint_1 = A * cvx.log(y) >= b
constraint_2 = cvx.sum(y) == c
print "constraint 1:"
print "left hand side: ", constraint_1.lh_exp.curvature
print "right hand side: ", constraint_1.rh_exp.curvature
print "DCP: ", constraint_1.is_dcp()
print "\nconstraint 2:"
print "left hand side: ", constraint_2.lh_exp.curvature
print "right hand side: ", constraint_2.rh_exp.curvature
print "DCP: ", constraint_2.is_dcp()
So we've fixed the second constraint, but broken the first one. In this case, the variables are small enough to print out the full expression for the first constraint.
In [5]:
print "Right hand side:\n", constraint_1.rh_exp
print "\nLeft hand side:\n", constraint_1.lh_exp
Notice that the right and left hand sides are switched from how the constraint was entered; CVXPY does this in order to represent all constraints as $\text{LHS } \leq \text{ RHS}$.
The matrix multiplying $\ln\left(y\right)$ is all positive, so the right hand side should be concave. However, CVXPY requires the user to specify when a coefficient has a particular sign if that sign matters to the convexity of the expression. This is done with parameters.
In [6]:
Aparam = cvx.Parameter(m,n,nonneg=True)
Aparam.value = A
constraint_1 = Aparam * cvx.log(y) >= b
print "constraint 1:"
print "left hand side: ", constraint_1.lh_exp.curvature
print "right hand side: ", constraint_1.rh_exp.curvature
print "DCP: ", constraint_1.is_dcp()
Now that both constraints are good, lets check the objective and solve the problem
In [7]:
objectiveExpression = cvx.sum(cvx.log(y))
print "objective Curvature: ", objectiveExpression.curvature
All right, everything should work now.
In [8]:
constraints = [constraint_1, constraint_2]
objective = cvx.Maximize(objectiveExpression)
problem = cvx.Problem(objective, constraints)
print "DCP rules test: ", problem.is_dcp()
try:
print "solving... "
problem.solve()
# print some stuff
print problem.status
print "Optimal value: ", problem.value
print "\nResidual of equality constraint: ", abs(sum(y.value) - c)
print "\nResidual of inequality constraints:\n", Aparam.value.dot(numpy.log(y.value)) - b
print "\ny:\n", y.value
print "\nx:\n", numpy.log(y.value)
except Exception as e:
print(e)