{
"metadata": {
"name": "",
"signature": "sha256:78d0399fcdd9547fc8d70f4f87a4a351fdc82b8db5387566ce4a2e097d644778"
},
"nbformat": 3,
"nbformat_minor": 0,
"worksheets": [
{
"cells": [
{
"cell_type": "heading",
"level": 2,
"metadata": {},
"source": [
"DCP Errors with Constraints"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Created by David Stonestrom 6/4/2014.\n",
"\n",
"
\n",
"The purpose of this document is to demonstrate debugging CVXPY DCP errors resulting from constraints. \n",
"\n",
"
Let's look at an optimization problem with a simple objective and two constraints. \n",
"$$ \\textbf{maximize } f\\left(x\\right) = \\sum_{i=1}^{n}{x_i}$$\n",
"\n",
"$$\\text{subject to } Ax \\succeq b$$\n",
"\n",
"$$\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\sum_{i = 1}^{n}{e^{x_i}} = c$$\n",
"\n",
"
\n",
"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."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"import cvxpy as cvx\n",
"import numpy\n",
"\n",
"m = 3; # number of inequality constraints\n",
"n = 6; # number of variables\n",
"\n",
"# generate some fake data\n",
"numpy.random.seed(0) # for repeatibility\n",
"A = numpy.random.random([m,n])\n",
"b = numpy.random.randn(m,1) + 0.5\n",
"c = 100 * numpy.random.random(1)\n",
"\n",
"#objective\n",
"x = cvx.Variable(n,1)\n",
"objectiveExpression = cvx.sum(x)\n",
"\n",
"# writing the constraints individually with names will make debugging them easier\n",
"constraint_1 = A*x >= b\n",
"constraint_2 = cvx.sum(cvx.exp(x)) == c\n",
"\n",
"# try the problem as written\n",
"objective = cvx.Maximize(objectiveExpression)\n",
"constraints = [constraint_1, constraint_2]\n",
"problem = cvx.Problem(objective, constraints)\n",
"\n",
"# Try solving the problem. Print error if any.\n",
"try:\n",
" problem.solve()\n",
"except Exception as e:\n",
" print e"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"Problem does not follow DCP rules.\n"
]
}
],
"prompt_number": 1
},
{
"cell_type": "heading",
"level": 3,
"metadata": {},
"source": [
"Diagnosing the DCP error"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"All right, lets look at the three components of the problem."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"print \"Objective: \", objective.is_dcp()\n",
"print \"Constraint 1: \", constraint_1.is_dcp()\n",
"print \"Constraint 2: \", constraint_2.is_dcp()"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"Objective: True\n",
"Constraint 1: True\n",
"Constraint 2: False\n"
]
}
],
"prompt_number": 2
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Okay, we can dig into constraint 2 further by examining the curvature on each side."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"print \"Constraint 2 curvatures:\"\n",
"print constraint_2.lh_exp.curvature\n",
"print constraint_2.rh_exp.curvature"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"Constraint 2 curvatures:\n",
"CONVEX\n",
"CONSTANT\n"
]
}
],
"prompt_number": 3
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Since the second constraint is an equality constraint, it needs each side to be AFFINE or CONSTANT in order to be DCP compliant. "
]
},
{
"cell_type": "heading",
"level": 3,
"metadata": {},
"source": [
"Change of variables"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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:\n",
"\n",
"$$ \\textbf{maximize } g\\left(y\\right) = \\sum_{i=1}^{n}{log\\left(y\\right)}$$\n",
"\n",
"$$\\text{subject to } A\\log\\left(y\\right) \\succeq b$$\n",
"\n",
"$$\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\sum_{i = 1}^{n}{y} = c$$\n",
"\n",
"\n",
"Lets look at how this transforms the constraints. "
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"y = cvx.Variable(n,1) # y = exp(x)\n",
"constraint_1 = A * cvx.log(y) >= b\n",
"constraint_2 = cvx.sum(y) == c\n",
"\n",
"print \"constraint 1:\"\n",
"print \"left hand side: \", constraint_1.lh_exp.curvature\n",
"print \"right hand side: \", constraint_1.rh_exp.curvature\n",
"print \"DCP: \", constraint_1.is_dcp()\n",
"\n",
"print \"\\nconstraint 2:\"\n",
"print \"left hand side: \", constraint_2.lh_exp.curvature\n",
"print \"right hand side: \", constraint_2.rh_exp.curvature\n",
"print \"DCP: \", constraint_2.is_dcp()\n"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"constraint 1:\n",
"left hand side: CONSTANT\n",
"right hand side: UNKNOWN\n",
"DCP: False\n",
"\n",
"constraint 2:\n",
"left hand side: AFFINE\n",
"right hand side: CONSTANT\n",
"DCP: True\n"
]
}
],
"prompt_number": 4
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"print \"Right hand side:\\n\", constraint_1.rh_exp\n",
"print \"\\nLeft hand side:\\n\", constraint_1.lh_exp"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"Right hand side:\n",
"[[ 0.5488135 0.71518937 0.60276338 0.54488318 0.4236548 0.64589411]\n",
" [ 0.43758721 0.891773 0.96366276 0.38344152 0.79172504 0.52889492]\n",
" [ 0.56804456 0.92559664 0.07103606 0.0871293 0.0202184 0.83261985]] * log(var3)\n",
"\n",
"Left hand side:\n",
"[[ 0.94386323]\n",
" [ 0.83367433]\n",
" [ 1.99407907]]\n"
]
}
],
"prompt_number": 5
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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}$. \n",
"\n",
"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."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"Aparam = cvx.Parameter(m,n,nonneg=True)\n",
"Aparam.value = A\n",
"\n",
"constraint_1 = Aparam * cvx.log(y) >= b\n",
"\n",
"print \"constraint 1:\"\n",
"print \"left hand side: \", constraint_1.lh_exp.curvature\n",
"print \"right hand side: \", constraint_1.rh_exp.curvature\n",
"print \"DCP: \", constraint_1.is_dcp()"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"constraint 1:\n",
"left hand side: CONSTANT\n",
"right hand side: CONCAVE\n",
"DCP: True\n"
]
}
],
"prompt_number": 6
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now that both constraints are good, lets check the objective and solve the problem"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"objectiveExpression = cvx.sum(cvx.log(y))\n",
"print \"objective Curvature: \", objectiveExpression.curvature"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"objective Curvature: CONCAVE\n"
]
}
],
"prompt_number": 7
},
{
"cell_type": "heading",
"level": 3,
"metadata": {},
"source": [
"Final solution"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"All right, everything should work now."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"constraints = [constraint_1, constraint_2]\n",
"objective = cvx.Maximize(objectiveExpression)\n",
"problem = cvx.Problem(objective, constraints)\n",
"\n",
"print \"DCP rules test: \", problem.is_dcp()\n",
"try: \n",
" print \"solving... \"\n",
" problem.solve()\n",
"\n",
" # print some stuff\n",
" print problem.status\n",
" print \"Optimal value: \", problem.value\n",
" print \"\\nResidual of equality constraint: \", abs(sum(y.value) - c)\n",
" print \"\\nResidual of inequality constraints:\\n\", Aparam.value.dot(numpy.log(y.value)) - b\n",
" print \"\\ny:\\n\", y.value\n",
" print \"\\nx:\\n\", numpy.log(y.value) \n",
" \n",
"except Exception as e:\n",
" print(e)\n"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"DCP rules test: True\n",
"solving... \n",
"optimal"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n",
"Optimal value: 4.01095525084\n",
"\n",
"Residual of equality constraint: [[ 0.00134525]]\n",
"\n",
"Residual of inequality constraints:\n",
"[[ 1.44312355e+00]\n",
" [ 1.82027513e+00]\n",
" [ 9.80742961e-04]]\n",
"\n",
"y:\n",
"[[ 2.08613192]\n",
" [ 2.3558799 ]\n",
" [ 1.71047705]\n",
" [ 1.72276186]\n",
" [ 1.66771014]\n",
" [ 2.28582697]]\n",
"\n",
"x:\n",
"[[ 0.73531159]\n",
" [ 0.85691429]\n",
" [ 0.53677231]\n",
" [ 0.54392874]\n",
" [ 0.51145151]\n",
" [ 0.82672787]]\n"
]
}
],
"prompt_number": 8
}
],
"metadata": {}
}
]
}