{
"metadata": {
"name": "",
"signature": "sha256:985eecf19ba0222b4ddd21d85cb036f7d9f7a7034e5c728f43b645208a603e6f"
},
"nbformat": 3,
"nbformat_minor": 0,
"worksheets": [
{
"cells": [
{
"cell_type": "heading",
"level": 2,
"metadata": {},
"source": [
"Dealing with DCP Errors"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Created by David Stonestrom 6/4/2014.\n",
"\n",
"
\n",
"The purpose of this document is to walk through the debugging process for DCP errors in CVXPY with a very simple problem. \n",
"\n",
"
Problem Statement:
\n",
"$$\\textbf{minimize } f\\left(x\\right) = \\sqrt{1 + ax^2} - \\frac{x}{4}$$ \n",
"\n",
"with no constraints. Note: the optimization variable $x$ is a scalar and the parameter $a \\geq 0$ is problem data. The first attempt is to just type it in and press go."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"import cvxpy as cvx\n",
"\n",
"# toy problem: minimize [sqrt(1 + a*x^2) - x/4] with parameter a > 0\n",
"x = cvx.Variable()\n",
"a = cvx.Parameter(nonneg=True)\n",
"\n",
"a.value = 0.1; # Eventualy we might want a tradeoff curve for a, but that comes after we can solve at least one instance.\n",
"\n",
"# normaly, I don't bother giving the objective expression its own name until things go wrong, but I'll save us some time here.\n",
"objectiveExpression = cvx.sqrt(1 + a * cvx.square(x)) - x / 4\n",
"\n",
"# instead I would just encode the objective function in this line\n",
"objective = cvx.Minimize(objectiveExpression);\n",
"\n",
"# let's solve the problem:\n",
"\n",
"problem = cvx.Problem(objective)\n",
"\n",
"# Try solving the problem. Print error if any.\n",
"try:\n",
" problem.solve()\n",
"except Exception as e:\n",
" print e\n"
],
"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": [
"Diognosing the problem"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This is only easy to do with small funcitons, but I like to plot the functions which may be causing problems. In this case the only function is the objective."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"import numpy, matplotlib.pyplot\n",
"%matplotlib inline\n",
"x = numpy.linspace(0,8,100)\n",
"matplotlib.pyplot.plot(x, numpy.sqrt(1 + a.value * x**2) - x / 4)"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 2,
"text": [
"[]"
]
},
{
"metadata": {},
"output_type": "display_data",
"png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD/CAYAAADllv3BAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xt4VNW5x/FvCPfoIAaLFi8lN/CGGoRggpghCsEKCAhF\nUUvEcCJaq6bp4dhaUVpRTo9p6w0FDVqF2odLbFCKEYnQBCJBJIqWQLgIgoggJiFCArPOH0umXCcJ\nTLJnJr/P8+zHTPbsmTcI76y8e613hRljDCIiEjJaOB2AiIj4lxK7iEiIUWIXEQkxSuwiIiFGiV1E\nJMQosYuIhJh6Jfbi4mLcbvdx38/Ly6N3794kJiYyY8YMADweDxkZGSQmJuJ2uykvL/dvxCIi4lPL\nup4wdepUXn/9dc4444yjvl9bW8tDDz1ESUkJ7du3JykpiSFDhvCvf/2LAwcOUFRURHFxMZmZmeTm\n5jbaDyAiIkerc8QeExPDvHnzOHYd0+eff05MTAwdOnSgVatW9O3bl6VLl1JYWMigQYMASEhIoKSk\npHEiFxGRE6ozsQ8fPpyWLY8f2FdUVNChQwfv4zPPPJPvvvuOiooKXC6X9/vh4eF4PB4/hSsiInU5\n5ZunHTp0oLKy0vu4srKSs846C5fLddT3PR4PLVroHq2ISFOps8Z+Mt27d2f9+vV8++23REREsHTp\nUrKysggLCyMvL4+RI0eyYsUKevToccLr27SJoaZGN1ZFRBoiOjqaDRs2+H6SqYdNmzaZa665xhhj\nzKxZs8xLL71kjDEmLy/P9OrVy/Ts2dM8//zzxhhjPB6PycjIMImJiSYxMdGsW7fuhK8JmB076vPu\nznr00UedDqFeFKf/BEOMxihOfwuWOOuTtus1Yv/JT35CUVERALfeeqv3+zfddBM33XTTUc8NCwvj\nhRdeqNcnz9tvw7hx9XqqiIjUk6PF73/8w8l3FxEJTY4m9iVLoLrayQjqlpyc7HQI9aI4/ScYYgTF\n6W/BEmd9hP1Qs2n6Nw4Lw+02PPggDB7sRAQiIsEnLCzsuHVFx3J0xD5kiMoxIiL+5uiIfcMGQ1IS\nbN8OmuouIlK3gB+xR0dDp06wcqWTUYiIhBbHx8kqx4iI+FdAJPa33nI6ChGR0OF4Yu/dG3bvhrpW\nyIqISP04nthbtIChQ2H+fKcjEREJDY4ndoBhw5TYRUT8xdHpjoffuqYGOneGzz6D885zIhoRkeAQ\n8NMdD2vdGgYN0uwYERF/CIjEDirHiIj4S0CUYgAqK6FLF9i6FY7YcU9ERI4QNKUYgDPPhH79bI92\nERE5dQGT2MGWY3JznY5CRCS4BUwpBmDXLoiNha++grZtnYhKRCSwBVUpBuCcc+CKKyA/3+lIRESC\nV0AldoARI2DuXKejEBEJXj4Tu8fjISMjg8TERNxuN+Xl5Uednz17NvHx8SQmJpKdne39fnx8PG63\nG7fbzbgG7lY9fDjk5dlFSyIi0nAtfZ3Mzc2lpqaGoqIiiouLyczMJPeHu5u7d+/m4YcfZvXq1XTo\n0AG3201ycjIXX3wxAEuWLDmlgM4/H+Li7H6oAwee0kuIiDRrPkfshYWFpKamApCQkEBJSYn3XHl5\nOVdccQVnnXUWYWFh9OnTh6VLl1JaWkp1dTUDBw4kJSWF4uLiBgd1yy0qx4iInCqfib2iogKXy+V9\nHB4ejsfjASA2Npa1a9fy9ddfU11dzeLFi6murqZ9+/ZkZWWxaNEipk2bxpgxY7zX1Nfw4Xba48GD\np/ATiYg0cz5LMS6Xi8rKSu9jj8dDix82J+3YsSPZ2dmMGDGCyMhI4uPj6dSpE3FxccTExAA2+UdG\nRrJjxw66dOly3OtPmjTJ+3VycjLJyckAdO0KF14Iy5aB2326P6KISPAqKCigoKCgYRcZH+bOnWvG\njh1rjDFm+fLl5sYbb/Seq62tNY8++qgxxpj9+/ebhIQEU15ebqZNm2YmTJhgjDHmyy+/NN27dzeH\nDh067rXreGvzxBPG/PAyIiLyg7pypzHG+FygZIxhwoQJlJaWApCTk8OqVauoqqoiPT2dyZMnk5ub\nS3h4OBkZGdx1110cPHiQtLQ0tmzZAsDUqVPp06fPca9d1yT7sjJIToZt2+xmHCIiUr8FSgG18vRY\nPXrA889D375NFJSISIALupWnx7rlFpgzx+koRESCS0CP2D/7DAYMgC++UDlGRARCYMR+ySXQsSMU\nFTkdiYhI8AjoxA7ws5/Bm286HYWISPAI6FIM2Nkx111nZ8eEhzdBYCIiASzoSzFg+8ace65drCQi\nInUL+MQOthzz9787HYWISHAI+FIMwMaNcM018OWX0NJnEwQRkdAWEqUYgKgo2zvmgw+cjkREJPAF\nRWIHGDVKs2NEROojKEoxAFu2wNVXw/bt0KpVIwYmIhLAQqYUA3DRRRAbq42uRUTqEjSJHeC222D2\nbKejEBEJbEFTigHYuRO6dbPlmPbtGykwEZEAFlKlGIDOnSEhAfLynI5ERCRwBVViB1uOmTXL6ShE\nRAJXUJViACoq4IILYPNm2/lRRKQ5CblSDIDLZXu0z53rdCQiIoEp6BI7wK23qhwjInIyQVeKAdi/\nH378Y/jkE+jSxc+BiYgEsNMuxXg8HjIyMkhMTMTtdlNeXn7U+dmzZxMfH09iYiLZ2dn1usYf2raF\nm2+Gv/3N7y8tIhL0fCb23NxcampqKCoq4sknnyQzM9N7bvfu3Tz88MO8//77FBYW8tZbb7F69Wpy\nc3M5cODACa/xpzvugL/+tVFeWkQkqPlM7IWFhaSmpgKQkJBASUmJ91x5eTlXXHEFZ511FmFhYfTp\n04elS5dSWFjIoEGDTniNP113HezZY8sxIiLyHz4Te0VFBS6Xy/s4PDwcj8cDQGxsLGvXruXrr7+m\nurqaxYsXs2/fPp/X+DXwFjBmjEbtIiLH8rlthcvlorKy0vvY4/HQooX9LOjYsSPZ2dmMGDGCyMhI\n4uPj6dSpE7t37z7pNceaNGmS9+vk5GSSk5MbFPwdd8ANN8CUKdoPVURCU0FBAQUFBQ27yPgwd+5c\nM3bsWGOMMcuXLzc33nij91xtba159NFHjTHG7N+/3yQkJJjy8nKf1xypjreut549jXn3Xb+8lIhI\nwKtP7vQ5Yh82bBj5+fkkJSUBkJOTw+zZs6mqqiI9PZ3w8HB69uxJeHg4GRkZREVF0bVr1+OuaUyH\nb6LecEOjvo2ISNAIynnsR/r6a4iLg23b4Iwz/BCYiEgAC8mWAsf60Y8gKQnmz3c6EhGRwBD0iR3g\nzjvhtdecjkJEJDAEfSkG4Pvv4fzzYfVquPBCv7ykiEhAahalGIB27WDUKI3aRUQgREbsAB9+aDfh\nWL8ewsL89rIiIgGl2YzYAXr1gjZt4F//cjoSERFnhUxiDwuDtDRo5GnzIiIBL2RKMQBffQUXXwxb\nt2pOu4iEpmZVigE491zo2xfmzHE6EhER54RUYgeVY0REQqoUA1BTY+e0FxVBTIzfX15ExFHNrhQD\n0Lo13H47vPKK05GIiDgj5EbsAJ9/DikpsGULtGrVKG8hIuKIZjliBzszJioK3n7b6UhERJpeSCZ2\ngPR0mD7d6ShERJpeSJZiAPbtgwsugDVr7H9FREJBsy3FAEREwOjRmvooIs1PyI7Ywbbxvflm2LhR\nm12LSGho1iN2gKuugk6dID/f6UhERJpOSCd2gPHj4cUXnY5CRKTp+EzsHo+HjIwMEhMTcbvdlJeX\nH3V+/vz59OrVi969ezNt2jTv9+Pj43G73bjdbsaNG9c4kdfTbbfBBx/Yza5FRJoDnzX2efPmsWDB\nAl555RWKi4uZMmUKubm53vNdu3Zl9erVREREcMkll1BSUkKbNm1ITEzko48+8v3GTVBjP+zee21J\n5rHHmuTtREQazWnX2AsLC0lNTQUgISGBkpKSo863atWKvXv3Ul1djTGGsLAw1qxZQ3V1NQMHDiQl\nJYXi4uLT/DFO3z33wIwZUFvrdCQiIo3PZ2KvqKjA5XJ5H4eHh+PxeLyPMzMz6dmzJ5dffjmDBw/G\n5XIRERFBVlYWixYtYtq0aYwZM+aoa5xw2WV2Jeo//uFoGCIiTaKlr5Mul4vKykrvY4/HQ4sW9rPg\niy++4Nlnn2XLli20b9+e22+/nTlz5jBkyBBifmirGBsbS2RkJDt27KBLly7Hvf6kSZO8XycnJ5Oc\nnOyHH+nE7rkHXngBRoxotLcQEfG7goICCgoKGnRNnTX2vLw8cnJyWLFiBZMnT+btHxqwlJWVMWrU\nKFauXEmrVq144IEHuPTSS/F4PJSWlvLcc8+xfft2UlJSWLt2rfcDwfvGTVhjBzhwAC68EJYtg7i4\nJntbERG/qk/u9JnYjTFMmDCB0tJSAHJycli1ahVVVVWkp6eTnZ3NrFmzaNu2LTExMUz/oTlLWloa\nW7ZsAWDq1Kn06dPnlILzt4kTbb/2p59u0rcVEfGb007sjcmJxL5pE/TqBV98Ae3bN+lbi4j4RbNf\neXqsrl0hMRHeeMPpSEREGk+zSuwAv/gFPPMMOPN7iohI42t2if366+HgQbsaVUQkFDW7xB4WBvfd\nZ0ftIiKhqFndPD2sqgouugg++sj+V0QkWOjm6UmccQbceaddsCQiEmqa5YgdYMMGuOYaO/WxXTvH\nwhARaRCN2H2IiYGEBHj9dacjERHxr2ab2AEefBCys8HhHmUiIn7VrBN7//7QujUsWuR0JCIi/tOs\nE3tYGDz0kHrHiEhoabY3Tw+rqbGtBhYuhB49nI5GRMQ33Tyth9at7YIljdpFJFQ0+xE7wJ49EB0N\nn30G553ndDQiIienEXs9nX023HYbPPus05GIiJw+jdh/cHjB0saNcOaZTkcjInJiGrE3QEwMuN3w\nwyZQIiJBSyP2I6xaBTffDOXl9qaqiEig0Yi9gXr2hO7dtcOSiAQ3jdiPsXixnf64di200MeeiASY\n0x6xezweMjIySExMxO12U15eftT5+fPn06tXL3r37s20adPqdU2g698fIiIgL8/pSERETk1LXydz\nc3OpqamhqKiI4uJiMjMzyc3N9Z5/6KGHWL16NREREVxyySWMHj2a999/nwMHDpz0mkAXFgYTJ8KT\nT8KQIfaxiEgw8TliLywsJDU1FYCEhARKSkqOOt+qVSv27t1LdXU1xhjCwsIoLCxk0KBBJ70mGAwb\nBrt3Q0GB05GIiDScz8ReUVGBy+XyPg4PD8dzRI/bzMxMevbsyeWXX87gwYPp0KFDndcEg/BwePhh\nmDzZ6UhERBrOZynG5XJRWVnpfezxeGjxwx3FL774gmeffZYtW7bQvn17br/9dubMmePzmmNNmjTJ\n+3VycjLJycmn8aP415gx8PjjUFgISUlORyMizVVBQQEFDSwf+EzsSUlJ5OXlMXLkSFasWEGPI9of\n7t+/n/DwcNq0aUOLFi340Y9+xN69e31ec6wjE3ugadXK1tp//3vb+VFExAnHDnofe+yxOq/xOd3R\nGMOECRMoLS0FICcnh1WrVlFVVUV6ejrZ2dnMmjWLtm3bEhMTw/Tp0wkPDz/umri4uOPfOECnOx7p\nwAG7InXePOjVy+loRETqlzs1j70OzzwD770Hb73ldCQiIkrsfvH99xAVBf/8J1xxhdPRiEhzp5YC\nftCuHWRl2RupIiLBQCP2eqiuhthYWLAArrrK6WhEpDnTiN1P2re3M2R+9zunIxERqZtG7PW0f78d\ntc+ZAwkJTkcjIs2VRux+1LYt/OY38OijTkciIuKbEnsD3HUXrFtnV6OKiAQqJfYGaN0aHnkEfvtb\nCKIqkog0M0rsDXTnnbB9O7z7rtORiIicmBJ7A7VsCX/4g50lE2RNK0WkmVBiPwUjRtgmYW++6XQk\nIiLH03THU/T++5CeDp9/bmvvIiJNQdMdG1H//rbz4/TpTkciInI0jdhPw+rVcOONsH49nHGG09GI\nSHOgEXsju+oqO3L/3/91OhIRkf/QiP00ffGFTfBr1sD55zsdjYiEOvVjbyIPPwxffgmvvup0JCIS\n6pTYm0hFBXTrZtv69uzpdDQiEspUY28iLhc89hhkZqrVgIg4T4ndT+66C3bv1t6oIuI8n6UYj8fD\nhAkTKC0tpU2bNsyYMYPo6GgAdu7cyejRo73P/fjjj3nqqacYP3488fHxdOjQAYCoqChefvnl4984\nhEoxh+XnQ0YGrF1r2/yKiPjbadfY582bx4IFC3jllVcoLi5mypQp5ObmHve85cuX88gjj5Cfn8+B\nAwdITEzko48+Ou3ggtGwYXD11bZ3u4iIv512jb2wsJDU1FQAEhISKCkpOe45xhjuv/9+XnjhBcLC\nwlizZg3V1dUMHDiQlJQUiouLT+NHCD5PPw3Z2bB1q9ORiEhz5TOxV1RU4HK5vI/Dw8PxHNPSMC8v\nj8suu4zY2FgAIiIiyMrKYtGiRUybNo0xY8Ycd00o69oV7r0XfvUrpyMRkeaqpa+TLpeLyspK72OP\nx0OLFkd/Frzxxhs88MAD3sdxcXHExMQAEBsbS2RkJDt27KBLly7Hvf6kSZO8XycnJ5OcnHwqP0PA\n+e//hksugSVLwO12OhoRCWYFBQUUFBQ06Jo6a+x5eXnk5OSwYsUKJk+ezNtvv33Uc6KjoykvL/c+\nfvHFFyktLeW5555j+/btpKSksHbt2uM+EEK1xn7Y3LkwaRJ89JFt8Ssi4g+nffPUGOOdFQOQk5PD\nqlWrqKqqIj09nV27djFw4MCjbpQePHiQtLQ0tmzZAsDUqVPp06fPKQUXzIyB1FS4/nrIynI6GhEJ\nFVp56rDyckhIgFWr4KKLnI5GREKBVp46LDoaHnwQ7rtPK1JFpOkosTeyrCw7cp8/3+lIRKS5UCmm\nCSxdCmPGwGefwZlnOh2NiAQz1dgDyLhx0K4dPPus05GISDBTYg8g334Ll10Gf/sbXHut09GISLDS\nzdMA0rGjHa2PGwfff+90NCISyjRib2KjRtm2A0895XQkIhKMVIoJQDt3Qo8e8PbbtgukiEhDqBQT\ngDp3hv/7Pxg7FvbvdzoaEQlFGrE7wBgYMQJiYmDqVKejEZFgolJMANu1C664At58U7NkRKT+VIoJ\nYOecA9Omwc9/Dkd0RhYROan6joWV2B00ZIjt1/7QQ05HIiKB7tNPISWlfs9VYndYdjYsXgwn2EpW\nRIRvv4X774f+/WH48Ppdo8TuMJcL3ngD/uu/YNs2p6MRkUBx6BC8+CJcfDEcPGh7Td13X/2u1c3T\nAPGHP0B+vh29h4c7HY2IOOmDD+CXv4SzzoI//QmuvPI/5zQrJogcOmR3W7r+evjNb5yORkScsHkz\n/PrX8OGH8Mc/2mnRYWFHP0ezYoJIeDi8/jo88wwUFTkdjYg0paoqO6C7+mq7Mv3zz+GWW45P6vWl\nxB5AunSB6dNh9Gg7z11EQpvHAzk50K0bbN0Ka9bAb39rW3yfDpViAtDEifDRR7BwoertIqFq6VK7\ndWabNraO3rt3/a477VKMx+MhIyODxMRE3G435eXl3nM7d+7E7XZ7j44dO/LSSy9hjDnpNVI/v/89\n1NTA5MlORyIi/rZhg62d33GH3TqzsLD+Sb3ejA9z5841aWlpxhhjVqxYYYYOHXrC5xUVFZmUlBTj\n8XjM3LlzzdixY+u8po63bvZ27DDmxz82ZuFCpyMREX/Ys8eYhx4yJjLSmCeeMKa6+tRepz650+eI\nvbCwkNTUVAASEhIoKSk50QcD999/Py+88AJhYWEUFhYyaNAgn9dI3c49F2bPtl0gN250OhoROVU1\nNfDnP9s6emUlrF0L//M/p19H98VnYq+oqMDlcnkfh4eH4/F4jnpOXl4el112GbGxsfW+RuqnXz97\np/zmm+1dcxEJHsbAvHlw6aWwaBG8/z689JJt3d3YWvo66XK5qDyiQ5XH46FFi6M/C9544w0eeOCB\nBl1z2KRJk7xfJycnk5yc3JDYm4X77oOPP7Yj97//HU7yRykiAaSoyNbPq6rg+efhhhtO/bUKCgoo\nKCho2EW+6jRH1suXL19ubrzxxuOeExUV1eBr6lsnEmv/fmP69DHm8cedjkREfCkrM2bECGPOP9+Y\nmTONOXjQ/+9Rn9zpc8Q+bNgw8vPzSUpKAiAnJ4fZs2dTVVVFeno6u3btokOHDnVeI6enTRv7K13v\n3vbXuvo2AhKRprFzJzz+uN1f4aGH4K9/bdwael00jz2IrFoFqal2v1S/T48SkQarqrJbXf7lL3Dn\nnfaeWKdOjfueaikQYnr2hJdftjdTN292OhqR5qumBp57DmJjoawMVq60LbgbO6nXl89SjASeIUNg\n0yb46U/twoazznI6IpHmw+Ox5ZZHHrFJfeHCozsvBgqVYoLU/ffb+bDvvGNr8CLSeIyxSfzhh6Ft\nW5gyxe5+5gS17Q1hhw7Bz35mu7/97W/qKSPSWJYtswl9zx67b8LQoafeddEfVGMPYeHhduelPXvg\nF7+o/ya3IlI/q1bBoEF2w/n0dCgttfe3nEzq9aXEHsTatIH586G4GB57zOloRELDp5/aJl1Dhtjj\n3/+2M16C6bdiJfYg53LZOvusWfD0005HIxK8yspgzBi7i1lSku3CeM890Lq105E1nBJ7COjc2e6V\n+uyzdgqWiNRfeblt2ZGUBJdcAuvX20VGTi4wOl2a7hgiLrjAJvfrrrMlmrvvdjoikcC2aZO9GZqb\na+9TbdgAxyykD1pK7CGka1eb3N1uaNnSjkJE5GibNsETT9g2HRMm2BLM2Wc7HZV/KbGHmNhYeO89\n202upgbGj3c6IpHAUF5u55/Pn28T+vr1oZfQD1NiD0Hdu0NBAaSkwP79djGTSHNVVmZH6AsW2DbY\nGzZAx45OR9W4lNhDVHQ0fPAB9O9vk/uvf+10RCJN65NPbEJfvPg/Cb25tODQrJgQdtFFNrnPnGkT\nuxYxSXOwciUMGwYDBkB8vC3B/O53zSepgxJ7yDv/fLsketkySEuD2lqnIxLxP2NgyRJ7b2nECDuB\nYONGu4vRmWc6HV3TU6+YZmLfPhg1yi6HfvNNiIhwOiKR0+fxwFtvwZNPwt69MHGiXWQUjIuK6ktN\nwOQotbV2lkxpKeTlwY9/7HREIqfmwAHbK2nqVLv6euJE25wrmJb9nyo1AZOjtGoFr7xif1Xt08du\nki0STPbutaPzrl3t5u7PP297JQ0f3jySen0psTczYWG2Bekf/2jrkf/4h9MRidRt82Z48EGIirL7\nECxcCP/8p531FQzdFpuaz+mOHo+HCRMmUFpaSps2bZgxYwbR0dHe8ytXriQzMxNjDF26dOG1116j\ndevWxMfHeze5joqK4uWXX27cn0IabNQoO2tmxAj46CM7a6CFPuYlwHz4oW1ul58Pd90Fa9bY9hni\nm8/EnpubS01NDUVFRRQXF5OZmUlubi4AxhjGjx/P3LlziYqKYvr06WzatImLLroIgCVLljR+9HJa\nEhKgpARGjrS9p19/PXR6ZUjwOnjQrg7905/gyy/hl7+El16ytXSpH59jtMLCQlJTUwFISEigpKTE\ne66srIzIyEiefvppkpOT2bt3L926dWPNmjVUV1czcOBAUlJSKC4ubtyfQE7LuefaBRwXXQS9eqnu\nLs7Zs8feDI2Ohr/8BTIz7aKiBx9UUm8on4m9oqIC1xF/ouHh4Xg8HgC++eYbioqK+MUvfsF7773H\n4sWLWbJkCREREWRlZbFo0SKmTZvGmDFjvNdIYGrd2rb8nTTJ1t2ff16LmaTplJba2VrR0fDZZ7Y5\n17Jl9oZoS62NPyU+E7vL5aKystL72OPx0OKHQmxkZCQxMTF069aNli1bkpqaSklJCXFxcYwZMwaA\n2NhYIiMj2bFjRyP+COIvt90GRUUwfbotz3z7rdMRSaiqrbWzWvr1s9vPnX++3alo5kzo2dPp6IKf\nz8/DpKQk8vLyGDlyJCtWrKBHjx7ec1FRUVRVVVFeXk50dDTLli3j7rvvJicnh9LSUp577jm2b99O\nRUUF55133glff9KkSd6vk5OTSU5O9ssPJacuNhaWL7ctCK64wk6PvP56p6OSULF1q62Xv/yy/bt2\n//12/nmrVk5HFrgKCgooKCho0DU+FygZY7yzYgBycnJYtWoVVVVVpKens2TJEiZOnIgxhqSkJLKz\nszl48CBpaWls2bIFgKlTp9KnT5/j31gLlALeu+/CuHG278aTT0L79k5HJMHo0CFYtAhefNGWWMaM\ngYwMuPRSpyMLTlp5Kqft22/h3nvt7Jnp0+0OTSL1sW2b/Y1vxgx7k378eBg9Gs44w+nIgpsSu/jN\nW2/Z1qeDBtmZC82pU57UX02N7Xs+YwasWAE/+5lN6Fdd5XRkoUMtBcRvhg6FTz+1sxQuvdT26dDn\nshxWWmo3gL7gAjv//NZb7Yj9hReU1J2gEbs02PLldvTevj088wxceaXTEYkTvv4aZs+G116DXbvg\n5z+3R0yM05GFNpVipNEcOmRnNjzyiL25+uijcJLJTxJCqqttf6HXX4d//QuGDIE77rA9W9SEq2mo\nFCONJjzc1k4//9xuZHDZZbbfzBHLHiRE1Nbaplt33GFbPc+caWvn27bZ0foNNyipBxqN2MUvNm+2\no/d334Vf/cruAq/NPILXwYN2Q/Q337R9W+LibN181Cjo3Nnp6Jo3lWKkyX36KTz+OCxdahN8Roam\ntwWL2lq7vdycOZCbCxdeaKcnjhxpewlJYFBiF8d88glMngzvv29LNvffb+cyS2DZt8/2Nc/NhXfe\nsatBR4607Zx/8hOno5MTUWIXx5WX237as2bZm6z33qteIE7bts1ujbhggV0J2qeP/X8zZAh06eJ0\ndFIXJXYJGLt22UUr06bZ2TMTJsAtt6hNQVOorbVTVBcutMfWrXah2eDBMHCgFpsFGyV2CTiHDsHb\nb9uFK8XF9lf+tDS45hptceYvxtg+5vn59mZ2QYFtiTtokD0SEtQON5gpsUtA+/JLOx965ky76/yo\nUfa46iol+YbassUm8PfftzdADx2y0xAHDICUFM1kCSVK7BIUjLF7Wf7973Z6Hdh67+DBcO21aul6\nrEOH7IYFhfMVAAAHtUlEQVTOy5fbGvnSpfaDsV8/m8T797c3QfXhGJqU2CXoHE7yeXn2WL8e3G7b\nE/7665tnwvryS1i50h4ffmiPc8+15au+fW1Cb45/Ls2VErsEva++suWF996zNWOPB5KS7JGYCD16\nQJs2TkfpH7W1tjb+ySd279mPP4bVq+33e/WyR+/edhZLp05ORytOUWKXkGIMbNoEhYW2T0lRkZ1O\n2b27nULZo4ftPHnppfCjHwXmCNYY2+N+40YoK4N16+yWcP/+t318/vm2PcOVV/7nuPDCwPxZxBlK\n7BLyqqtty9hVq+yq17Vr7eHx2Jkgh48LLrBztLt0sWWMTp38P9KvrbVJe+dOe3z1lS2jbN1qjy++\nsB9MxkBUlF2m362bPbp3h4sv1vRPqZsSuzRLxsDu3XY0f/jYts0m2W3bbLvZb76xiT0yElwu28js\nzDOhbVv7/TZt7JTAwyNlj8cm7poae6OyuhqqquzKzYoK2LPHfq9jRzsDpXNn+wFy3nn2Q+XwERUF\nZ5+tEbicOiV2kZMwxibk3bttR8rKSvv4wAF71NTYRlhHat36P0f79rYHzhln2A+Gs8+2HwxK2NLY\nlNhFRELMafdj93g8ZGRkkJiYiNvtpry8/KjzK1eupF+/flx77bWMHj2ampqaOq8REZHG5TOx5+bm\nUlNTQ1FREU8++SSZmZnec8YYxo8fz8yZM1m2bBkpKSls2rSJ3NxcDhw4cMJrglFBQYHTIdSL4vSf\nYIgRFKe/BUuc9eEzsRcWFpKamgpAQkICJSUl3nNlZWVERkby9NNPk5yczN69e+nWrRuFhYUMGjTo\nhNcEo2D5n604/ScYYgTF6W/BEmd9+GwFVFFRgcvl8j4ODw/H4/HQokULvvnmG4qKinjuueeIjo7m\npptu4uqrr/Z5jYiIND6fid3lclF5xCaWRyboyMhIYmJi6NatGwCpqamUlJT4vEZERJqA8WHu3Llm\n7Nixxhhjli9fbm688UbvuQMHDpiuXbuaDRs2GGOMGT58uHnnnXd8XnOk6OhoA+jQoUOHjgYc0dHR\nvtK2McYYn9MdjTFMmDCB0tJSAHJycli1ahVVVVWkp6ezZMkSJk6ciDGGpKQksrOzT3hNXFzcyd5C\nRET8zLF57CIi0jhU/BYRCTFNntiDbQFTcXExbrfb6TBOqLa2ljvuuIN+/fqRkJBAXl6e0yGd0KFD\nh7jrrrvo27cv1157LWvXrnU6JJ++/vprLrjgAsrKypwO5aTi4+Nxu9243W7GjRvndDgnNWXKFBIT\nE+nVqxevvvqq0+Gc0Kuvvur9s+zTpw/t2rWjoqLC6bCO4vF4vP+G+vXrx7p163xfUGcV3s/mzp1r\n0tLSjDHGrFixwgwdOrSpQ6i3p556ylx++eXmmmuucTqUE8rJyTEPPvigMcaYPXv2mAsvvNDhiE4s\nNzfXjBs3zhhjTEFBQUD/P6+pqTE333yz6datm1m3bp3T4ZzQ999/b6666iqnw6jTkiVLzODBg40x\nxlRVVZnf/e53DkdUt3vvvddMnz7d6TCOs3DhQjNq1ChjjDH5+flmxIgRPp/f5CN2X4ueAk1MTAzz\n5s0L2J42I0eO5PHHHwfsJ3rLAN2heOjQobz44osAbN68mY4dOzoc0cllZWVxzz33cN555zkdykmt\nWbOG6upqBg4cSEpKCsXFxU6HdELvvvsul19+OTfffDODBw9myJAhTofkU0lJCWvXruXuu+92OpTj\ntGvXju+++w5jDN999x2tW7f2+fwmzwTBtIBp+PDhbN682ekwTioiIgKAyspKRo4cyR/+8AeHIzq5\n8PBwxo4dy/z585kzZ47T4ZzQzJkzOeeccxgwYABTpkwJ2A/0iIgIsrKyGDduHOvXr2fQoEGUlZUF\n3L+hXbt2sXXrVhYsWMDGjRsZMmQI//73v50O66SeeOIJJk2a5HQYJ5SUlMT+/fvp3r07u3fvrrPs\n2uR/E7SAyb+2bt1K//79ufPOOxk9erTT4fg0c+ZMysrKSE9P5/vvv3c6nOPk5OSQn5+P2+3m448/\n5uc//zk7d+50OqzjxMXFMWbMGABiY2OJjIxkx44dDkd1vE6dOjFgwABatmxJXFwcbdu25ZtvvnE6\nrBPau3cvZWVlXHfddU6HckJTp04lKSmJdevWef9u1tTUnPT5TZ5Rk5KSeOeddwBYsWIFPXr0aOoQ\nQsbOnTsZMGAAU6dOZezYsU6Hc1J//etfmTJlCmB/pWzRokVAfph/8MEHFBQUsGTJEq688kpee+01\nOnfu7HRYx8nJyfE219u+fTsVFRUBWTrq27cv//znPwEb5759+4iMjHQ4qhNbunQpKSkpTodxUvv2\n7fNWOjp27EhtbS2HDh066fObvBQzbNgw8vPzSUpKAuxf0kAXFqC7JzzxxBN89913PP74495a+8KF\nC2nbtq3DkR3tlltuYezYsVx33XXU1tby5z//mTahsgO1A8aNG0daWhr9+vUD7L+hQPyg/OlPf8rS\npUvp3bs3Ho+H559/PmD/LZWVlREdHe10GCeVlZVFWloa1157LbW1tUyZMoV27dqd9PlaoCQiEmIC\n72NeREROixK7iEiIUWIXEQkxSuwiIiFGiV1EJMQosYuIhBgldhGREKPELiISYv4fbfP7NEY37rgA\nAAAASUVORK5CYII=\n",
"text": [
""
]
}
],
"prompt_number": 2
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Definitely looks convex, let's dig into the expression:"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"print(objective.is_dcp())\n",
"print(objectiveExpression.curvature)"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"False\n",
"UNKNOWN\n"
]
}
],
"prompt_number": 3
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Okay... why? First, lets check that the parameters have signs specified. Compare the declaration and results for a (defined above) and b."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"# a = cvx.Parameter(nonneg=True) # copied from above, sign is set here\n",
"# a.value = 0.1; \n",
"print \"a value: \", a.value, \"\\t a sign: \", a.sign\n",
"\n",
"# b is only included so you can see a common problem for DCP rules\n",
"b = cvx.Parameter() # notice that sign is not set here\n",
"b.value = 3\n",
"print \"b value: \", b.value, \"\\t b sign: \", b.sign\n"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"a value: 0.1 \t a sign: POSITIVE\n",
"b value: 3 \t b sign: UNKNOWN\n"
]
}
],
"prompt_number": 4
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"If `a` could be negative (sign NEGATIVE or UNKNOWN), the expression could be non-convex (plot it with a = -1 if you want to check this). \n",
"\n",
"\n",
"\n",
"Since that's not the problem, let's look at the objective expression next. It first breaks apart at the '-'"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"x = cvx.Variable() # remember we changed x to linspace when plotting? now we want it back as a cvxpy variable\n",
"\n",
"# objectiveExpression = cvx.sqrt(1 + a * cvx.square(x)) - x / 4 # copied from above\n",
"expression1 = cvx.sqrt(1 + a * cvx.square(x))\n",
"expression2 = - x / 4\n",
"print expression1.curvature\n",
"print expression2.curvature\n"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"UNKNOWN\n",
"AFFINE\n"
]
}
],
"prompt_number": 5
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"So, expression1 is the problem. The outermost function there is sqrt. Its argument is $1 + a x^2$. Note that $1 + ax^2$ is the function $+$ applied to arguments $1$ and $a x^2$, so its curvature depends on its monotonicity and the curvature of cvx.square(x)."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"print(\"cvx.sqrt:\")\n",
"print \"curvature: \", cvx.sqrt(x).curvature # note that curvature is a property \n",
"print \"monotonicity: \", cvx.sqrt(x).monotonicity() # but monotonicity() is a function\n",
"\n",
"expression3 = (1 + a * cvx.square(x))\n",
"print \"\\nexpression3: \", expression3\n",
"print \"curvature: \", expression3.curvature\n",
"print \"monotonicity: \", expression3.monotonicity()\n",
"\n",
"print(\"\\na * cvx.square:\")\n",
"print (a * cvx.square(x)).curvature"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"cvx.sqrt:\n",
"curvature: CONCAVE\n",
"monotonicity: ['INCREASING']\n",
"\n",
"expression3: 1 + param1 * square(var3)\n",
"curvature: CONVEX\n",
"monotonicity: ['INCREASING', 'INCREASING']\n",
"\n",
"a * cvx.square:\n",
"CONVEX\n"
]
}
],
"prompt_number": 6
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Starting from the inside: expression3 is convex since it is an affine function and it is increasing on the constant argument '1' and increasing on the convex argument 'cvx.square.'\n",
"\n",
"\n",
"Next we consider cvx.sqrt: It is an increasing function and the input is convex, so this seems good for convexity. The problem is that sqrt is a concave function. \n",
"\n",
"\n",
"\n"
]
},
{
"cell_type": "heading",
"level": 3,
"metadata": {},
"source": [
"Fixing the DCP error"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Since we're pretty confident that the overall function is convex, we should start looking for different phrasings of the problematic part of the objective that will be recognized as convex by the DCP ruleset. \n",
"\n",
"$$\\sqrt{1 + a*x^2} = \\sqrt{1^2 + \\left(\\sqrt{a}*x\\right)^2} = \\left(1^2 + \\left(\\sqrt{a}*x\\right)^2\\right)^\\frac{1}{2} = \\left(\\left[\\begin{array}{ll} 1 & \\sqrt{a}*x\\end{array}\\right] \\left[\\begin{array}{c} 1 \\\\ \\sqrt{a}*x\\end{array}\\right]\\right)^\\frac{1}{2}$$\n",
"\n",
"\n",
"\n",
"Hopefully, one of these popped out to you as $\\left|\\left|\\left[\\begin{array}{c} 1 \\\\ \\sqrt{a}*x\\end{array}\\right]\\right|\\right|_{2}$. We can rewrite the objective to use this cvxpy defined norm, and fix the DCP issues."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"x = cvx.Variable(2) # note x is now a vector, x[0] will be constrained to 1 and x[1] will be x from the earlier attempt\n",
"a = cvx.Parameter(nonneg=True)\n",
"\n",
"a.value = 0.1;\n",
"\n",
"#The following is one way of making the objective expression correct with the new vector x:\n",
"import math\n",
"b = cvx.Parameter(2,2,nonneg=True)\n",
"b.value = numpy.array([[1, 0], [0, math.sqrt(a.value)]]) # for a == 1, b is the 2x2 identity\n",
"\n",
"# the 2-norm replaces the first big expression and x[1] replaces the scalar x\n",
"objectiveExpression = cvx.norm(b*x, 2) - x[1] / 4 \n",
"\n",
"# this constraint forces x[0] to be 1, making the new problem match the origional.\n",
"constraint = [x[0] == 1]\n",
"\n",
"\n",
"# I prefer this version where x[0] is constrained to be 1/sqrt(a), and the whole x vector is multiplied \n",
"# by sqrt(a) before taking the norm. Note: the constant sqrt(a) can be pulled out of the norm expression.\n",
"objectiveExpression = cvx.sqrt(a)*cvx.norm(x, 2) - x[1] / 4 # no need for the parameter b\n",
"constraint = [x[0] == 1/cvx.sqrt(a)]\n",
"\n",
"# form objective\n",
"objective = cvx.Minimize(objectiveExpression);\n",
"problem = cvx.Problem(objective, constraint)\n",
"print \"DCP rules test: \", problem.is_dcp()\n",
"\n",
"# and solve it\n",
"try: \n",
" print \"solving... \"\n",
" problem.solve()\n",
"except Exception as e:\n",
" print(e)\n",
"\n",
"# Print results.\n",
"print problem.status\n",
"print \"Optimal value: \", problem.value\n",
"print \"Full variable x: \"\n",
"print x.value\n",
"print \"Origional optimization variable: x =\", x[1].value"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"DCP rules test: True\n",
"solving... \n",
"optimal\n",
"Optimal value: 0.612372434091\n",
"Full variable x: \n",
"[[ 3.16227766]\n",
" [ 4.08256597]]\n",
"Origional optimization variable: x = 4.08256597129\n"
]
}
],
"prompt_number": 7
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This is the happy outcome for DCP errors; the objective only needed to be restated with a cvxpy built-in function to make things work. The unhappy outcome is that the problem is nonconvex and requires either relaxation or or reformulation to make cvxpy work on it. We'll dive into this in later problems."
]
}
],
"metadata": {}
}
]
}