{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## Isotonic regression\n", "\n", "Problem data generation. (This code is in a Jupyter Notebook.)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "import cvxpy as cp\n", "import numpy as np\n", "from matplotlib import pyplot as plt\n", "%matplotlib inline\n", "\n", "np.random.seed(0)\n", "A = np.random.randn(20,20)\n", "x_true = np.cumsum(np.random.rand(20))\n", "b = A@x_true + 2*np.random.randn(20)\n", "\n", "\n", "def huber_loss(x) :\n", " return np.sum((x**2)*(np.abs(x)<=1)+(2*np.sign(x)*x-1)*(np.abs(x)>1))\n", "\n", "def huber_grad(x) :\n", " return 2*x*(np.abs(x)<=1) + 2*np.sign(x)*(np.abs(x)>1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "CVXPY solution" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "24.733210916870362\n", "[-0.28 -0.28 0.92 2.14 3.39 3.39 3.39 3.73 3.73 3.73 5.29 5.52 5.52 6.08\n", " 6.61 7.35 7.35 7.35 7.79 7.79]\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWoAAAD4CAYAAADFAawfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAXvElEQVR4nO3df4zkdX3H8ed7+aWrVtRbBcHbBWtVaoLCLFFpzY0cBq3xSv8YuSxiWs3GDVhpaq6Ey5yy2zPx2hrbSqbZqlXLFBwVKjFq5XrfizEpsHsICB4o4u15Ht4tNgJ2g4j37h+fmbvdYWb3O/ud73y/O/t6JJPvzPf7/ey+73vDi+99vp/v52vujoiI5NdA1gWIiMjyFNQiIjmnoBYRyTkFtYhIzimoRURy7uQ0fuiGDRt8ZGQkjR8tItKX9u3b97i7D7XalkpQj4yMMDs7m8aPFhHpS2Y2125brK4PM/srM3vQzB4ws5vN7HndK09ERJazYlCb2VnAXwIFd38DcBJwRdqFiYhIEPdi4snA883sZGAQOJxeSSIistiKQe3uPwf+HjgIPAY84e7fad7PzMbNbNbMZufn57tfqYjIOhWn6+MlwBbgHOCVwAvM7Mrm/dx92t0L7l4YGmp54VJERFYhTtfHZuCn7j7v7r8FbgXemm5ZIiK9sWsXRNHSdVEU1veifRxxgvog8GYzGzQzAy4B9nevBBGR7IyOQql0ImyjKHweHe1N+zji9FHfBXwVuAf4Qb3NdPdKEBHJTrEItVoI1x07wrJWC+t70T6OWKM+3P1j7v46d3+Du7/P3X/TvRJERLJVLMLEBExNhWWnIZu0/Uo014eIrHtRBJUKlMth2dznnHb7lSioRWRda/Qp12owOXmiGyNu2EYRlLY8Te2krUz+7QC1k7ZS2vJ0V8NaQS0i69rMzNI+5Uaf88xMzPbT36f2zOUUj9wC7hSP3ELtmcuZmf5+12q0NJ6ZWCgUXJMyici6MDICcy3mUxoehgMHYv8YM9vn7oVW23RGLSJSrYbAHRgIy2o1ftuDBztbvwoKahFZ36pVGB8PZ8XuYTk+Hj+sN27sbP0qKKhFZH3bvh0WFpauW1gI6+PYuRMGB5euGxwM67tEQS0iiWR9C3biW7iTdl2MjcH0dOiTNgvL6emwvksU1CKSSNa3YCe+hbsbXRdjY+HC4bFjYdnFkAbA3bv+uvDCC11E1o89e9w3bHAvl8Nyz54et7/+Dt8w8LiXmfQNA4/7nuvviN/4ppvcBwfdQw91eA0OhvU9BMx6m0xVUItIV5TLIVHK5R63rwdtmRtCe27oPGhvusl9eNjdLCx7HNLuCmoRSVmmZ9TDw76HTb6Bo17mBt/AUd/DphC4a4iCWkRS0wjZRrg2f069PcUT4QzHQ3sPxc7/MBlaLqh1MVFEEkl8C3bS9qdfSo0SRfaG9uylRomZ0y/t7A+SY7qFXETWtsYNK4vHQg8Odn2IXNoS3UJuZq81s3sXvZ40s2u7X6aIyCr0YBxz1uI84eVhd3+ju78RuBBYAG5LvTIRWT+SzLUB6Y9jztjJHe5/CfATd28xVZSIyCo0d1005tqAvgvc1er0YuIVwM1pFCIi61TSuTbWgdhBbWanAu8BvtJm+7iZzZrZ7Pz8fLfqE5F+14NpQte6Ts6o3wnc4+5HWm1092l3L7h7YWhoqDvViUj/68E0oWtdJ0G9FXV7iEi39WCa0LUuVlCb2SBwKXBruuWIyLqzDobXJRVr1Ie7LwAvS7kWEVmvxsYUzMvQLeQiIjmnoBZZ5xI/IUVSp6AWWecSPyFFUtfpnYki0mcas9WVSjAxAZXK0tnsJHs6oxYRisUQ0lNTYamQzhcFtYgQReFMulwOy+Y+a8mWglpknWv0SddqMDl5ohtEYZ0fCmqRdS7pE1aA5NOUyrL0hBcRSaZPnrCStURPeBERWZamKU2dglpEktE0palTUItIMpqmNHUKahFJRtOUpk5BLSLJaJrS1CmoRda4XEyq1OdPAc+aglpkjevKpEoaB51rcZ/wcrqZfdXMHjKz/Wb2lrQLE5F4ikWoje+mtPmX7LApSpt/SW18d/z5OhrjoOfmwD0sx8cV1jkS94z6H4Fvu/vrgPOB/emVJCIdqVYpfnoLE8duZIoyE8dupPjpLfGDVuOgc2/FOxPN7PeA+4BzPeZtjLozUaSHRkaI5s6hRI0JKlSYoEaJ4vBPQ3/xSgYGwpl0M7PQ5yw9kfTOxHOBeeDfzOz7ZvZZM3tBi18ybmazZjY7Pz+fsGQRiSuaO5cSNWqUmORj1ChRokY0d268H6Bx0LkXJ6hPBi4AKu7+JuD/gOuad3L3aXcvuHthaGioy2WKSDszp18azqDZC0CRvdQoMXP6pfF+gMZB516coD4EHHL3u+qfv0oIbhHJgW2f2Uhx8O4l64qDd7PtMzHPiDUOOvdWfBSXu//CzH5mZq9194eBS4Afpl+aiMTSCNTt28P8Ghs3hrPhToJ2bEzBnGNxn5n4YaBqZqcCjwJ/nl5JItIxBW1fixXU7n4v0PJqpIiIpEt3JoqI5JyCWiQPdAu3LCNuH7WIpKX5UVaNW7hB/c4C6IxaJLFduyDavnvJGXG0fXfs2et2XXOQaOGiJeuihYvYdY2ekCKBglokodEndlP6xPlEc+eAe7id+xPnM/rE7njtf3VHuJOQTQBEbKJEjdFf3ZFi1bKWKKhFEipWP3j8tu0d3HD8du5i9YPx2g8/2rr98KMpVy5rhYJaJKmDBymylwkqTLGDCSrhdu64D3fduZPi4N1L2w/erVu45TgFtUhSGzcSsYkKE5SZpMJE6MaIO6nR2BjRtV+nMnA1ZaaoDFxNdO3XdSFRjlNQy7qX9FFW0dhnW89eN/bZeO0jKE1vprb7ZUx6mdrul1Ga3vycmmT9UlDLupf0UVYzL95M7fr7wvzPZhSHf0rt+vuYefHmeO1noFbj+BNZisXweWZmFX8Y6UsrPjhgNfTgAFlrGuE8MQGVytLgFOmFpA8OEOl7xWII6ampsFRIS54oqEUIZ9SVCpTLYan+YckTBbWse41uj1oNJifDcnGftUjWFNSy7ulinuRdrIuJZnYAeAr4HfBsuw7vBl1MFBHpzHIXEzuZPa/o7o93qSYREYlJXR8iIjkXN6gd+I6Z7TOz8TQLEhGRpeJ2fVzs7ofN7OXAHWb2kLt/d/EO9QAfB9gYd44DERFZUawzanc/XF8eBW4DLmqxz7S7F9y9MDQ01N0qRUTWsRWD2sxeYGYvarwH3gE8kHZhIiISxOn6eAVwm5k19v8Pd/92qlWJiMhxKwa1uz8KnN+DWkREpAUNzxMRyTkFtYhIzimoRQCqVRgZgYGBsKxWs65I5LhObiEX6U/VKoyPw8JC+Dw3Fz6DnlsouaAzapHt20+EdMPCQlgvkgMKapGDBztbL9JjCmqRdlMeaCoEyQkFtcjOnTA4uHTd4GBYL5IDCmqRsTGYnobhYTALy+lpXUiU3NCoDxEIoaxglpzSGbWISM4pqEVEck5BLSKScwpqEZGcU1CLiOScglr6gyZVkj4We3iemZ0EzAI/d/d3p1eSSIc0qZL0uU7OqD8C7E+rEJFV06RK0udiBbWZnQ38CfDZdMsRWQVNqiR9Lu4Z9aeBbcCxdjuY2biZzZrZ7Pz8fFeKE4lFkypJn1sxqM3s3cBRd9+33H7uPu3uBXcvDA0Nda1AkRVpUiXpc3HOqC8G3mNmB4BbgLeb2U2pViXSCU2qJH3O3D3+zmabgI+uNOqjUCj47OxswtJERNYPM9vn7oVW2zSOWhLbtQuiaOm6KArr10J7kbzrKKjdfa/GUEuz0VEobXma6IytMDBAdMZWSlueZnS0h+1LJ8I6isLnuO1F8k7zUUtixcNVas/cROmpLzHB66kcmaB22uUUD18JrNxPnLh9EWq1EM4TE1CphM/FYvI/m0guuHvXXxdeeKHLOjI87A5e5gYH9zI3uENY34v2deVyaFYud1i/SA4As94mU9VHLckdPEjEJipMUGaSChNEbIp/w0nS9oTujkoFyuWwbO6zFlnT2iV4kpfOqNeXPa+4wjdw1PewyR18D5vC51dc0Zv2e9w3bAjLVp9F1gJ0Ri1pmiluo3baVRTZC0CRvdROu4qZ4rbetJ9Z2ifd6LOemen0TyKSTx2No45L46jXoWo1TIJ08GC4dXvnzs5uOEnaXmSN0zjqnMt6HHFXxiGPjcGBA3DsWFh2GrJJ24v0MQV1DiQdB5x1exFJWbvO6yQvXUzsXOMCWLm8ugthWbcXkWTQxcT8KxbDzRpTU2HZ6c0aWbcXkfQoqHMi6TjgrNuLSIranWoneanrozNJxwFn3V5EkkNdH/mWdBxw1u1FJF0aRy0ikgMaRy0isobFeWbi88zsbjO7z8weNLMbelGYiIgEceaj/g3wdnf/tZmdAnzPzL7l7nemXJuIiBDjjLp+QfLX9Y+n1F/d79iWZKpVGBmBgYGwrFZ7215EUhPrCS9mdhKwD/h94EZ3vyvVqqQz1SqMj8PCQvg8Nxc+Q7w5M5K2F5FUdfoU8tOB24APu/sDTdvGgXGAjRs3Xjg3N9fNOmU5IyMhXJsND4cJjtJuLyKJdW3Uh7v/CtgLXNZi27S7F9y9MDQ0tKpCZZXaPQmlgyesJGovIqmKM+pjqH4mjZk9H9gMPJR2YdKBjRs7W9/t9iKSqjhn1GcCkZndD8wAd7j7N9ItSzqycycMDi5dNzgY1veivYikasWLie5+P/CmHtQiq9W44LfaJ6QkbS8iqdIt5HmhR1GJrGvLXUyMNTxPUqbhcSKyDM31kQfbt58I6YaFhbBeRNY9BXUeaHiciCxDQZ0HGh4nIstQUOeBhseJyDIU1HkwNgbT0+GWbbOwnJ7WhUQRATTqIz/GxhTMItKSzqhFRHJOQS0iknMKahGRnFNQi4jknIJaRCTnFNRdsGsXRNHSdVEU1ouIJKWg7oLRUShteZrojK0wMEB0xlZKW55mdDTrykSkH2gcdRcUD1epPXMTpae+xASvp3Jkgtppl1M8fCWgsdEikkycR3G9yswiM9tvZg+a2Ud6Udiasn07xd98mwkqTLGDCSoUf/NtzX4nIl0Rp+vjWeCv3f31wJuBq83svHTLWmMOHiRiExUmKDNJhQkiNmn2OxHpihWD2t0fc/d76u+fAvYDZ6Vd2FoSvfy9lKhRo8QkH6NGiRI1ope/N+vSRKQPdHQx0cxGCM9PvKvFtnEzmzWz2fn5+e5Ut0bMFLdRO+0qiuwFoMheaqddxUxxW7aFiUhfiB3UZvZC4GvAte7+ZPN2d59294K7F4aGhrpZY+5tu/lNFD935ZLZ74qfu5JtN+uZwCKSXKxRH2Z2CiGkq+5+a7olrVGa/U5EUhJn1IcBnwP2u/un0i9JREQWi9P1cTHwPuDtZnZv/fWulOsSEZG6Fbs+3P17gPWgFhERaUG3kIuI5JyCWkQk5xTUIiI5p6AWEck5BbWISM4pqEVEck5BLSKScwpqEZGcU1CLiOScglpEJOcU1A3VKoyMwMBAWFarWVckIgLo4bZBtQrj47CwED7PzYXPoKlLRSRzOqOG8BDaRkg3LCzo4bQikgsKamj/EFo9nFZEckBBDbBxY2frRUR6KM4TXj5vZkfN7IFeFJSJnTthcHDpusHBsF5EJGNxzqi/AFyWch3ZGhuD6eklD6dleloXEkUkF+I84eW7ZjaSfikZ08NpRSSn1EcN7NoFUbR0XRSF9SIiWetaUJvZuJnNmtns/Px8t35sT4yOQql0IqyjKHweHc22LhER6GJQu/u0uxfcvTA0NNStH9sTxSLUaiGcd+wIy1otrBcRyZq6PuqKRZiYgKmpsFRIi0hexBmedzPwP8BrzeyQmX0g/bJ6L4qgUoFyOSyb+6xFRLISZ9TH1l4UkqVGn3Sju6NYVPeHiOSHuj6AmZmlodzos56ZybYuEREAc/eu/9BCoeCzs7Nd/7kiIv3KzPa5e6HVNp1Ri4jkXP8EtSb+F5E+1R8PDtDE/yLSx/rjjFoT/4tIH+uPoNbE/yLSx/ojqDXxv4j0sf4Iak38LyJ9rD+CWhP/i0gf649RH6CJ/0Wkb+XnjFrjoEVEWsrHGbXGQYuItJWPM2qNgxYRaSsfQa1x0CIibeUjqDUOWkSkrVhBbWaXmdnDZvaImV3X9So0DlpEpK04j+I6CbgReCdwHrDVzM7rahUaBy0i0lacUR8XAY+4+6MAZnYLsAX4YVcr0ThoEZGW4nR9nAX8bNHnQ/V1S5jZuJnNmtns/Px8t+oTEVn34gS1tVj3nOd3ufu0uxfcvTA0NJS8MhERAeIF9SHgVYs+nw0cTqccERFpFieoZ4DXmNk5ZnYqcAVwe7pliYhIw4pB7e7PAtcA/wXsB2ru/mA3i9i1C6Jo6booCut70V5EJM9ijaN292+6+x+4+6vdveuDm0dHoVQ6EbZRFD6PjvamvYhInuViUqZiEWq1EK4TE1CphM/FYm/ai4jkWT5uISeE6sQETE2FZachm7S9iEhe5SaooyicCZfLYdnc55x2exGRvMpFUDf6lGs1mJw80Y0RN2yTthcRybNcBPXMzNI+5Uaf88xMb9qLiOSZuT/nJsPECoWCz87Odv3nioj0KzPb5+6FVttycUYtIiLtKahFRHJOQS0iknMKahGRnFNQi4jkXCqjPsxsHphbZfMNwONdLKfbVF8yqi8Z1ZdMnusbdveWk/mnEtRJmNlsuyEqeaD6klF9yai+ZPJeXzvq+hARyTkFtYhIzuUxqKezLmAFqi8Z1ZeM6ksm7/W1lLs+ahERWSqPZ9QiIrKIglpEJOcyCWozu8zMHjazR8zsuhbbzcz+qb79fjO7oMf1vcrMIjPbb2YPmtlHWuyzycyeMLN7668dPa7xgJn9oP67nzNVYZbH0Mxeu+i43GtmT5rZtU379PT4mdnnzeyomT2waN1LzewOM/txffmSNm2X/b6mWN/fmdlD9b+/28zs9DZtl/0upFjfx83s54v+Dt/Vpm1Wx+/Li2o7YGb3tmmb+vFLzN17+gJOAn4CnAucCtwHnNe0z7uAbwEGvBm4q8c1nglcUH//IuBHLWrcBHyj18dv0e8/AGxYZnumx7Dp7/sXhMH8mR0/4G3ABcADi9btAq6rv78O+GSb+pf9vqZY3zuAk+vvP9mqvjjfhRTr+zjw0Rh//5kcv6bt/wDsyOr4JX1lcUZ9EfCIuz/q7s8AtwBbmvbZAnzJgzuB083szF4V6O6Pufs99fdPAfuBs3r1+7sk02O4yCXAT9x9tXeqdoW7fxf436bVW4Av1t9/EfjTFk3jfF9Tqc/dv+Puz9Y/3gmc3e3fG1eb4xdHZsevwcwMKAE3d/v39koWQX0W8LNFnw/x3BCMs09PmNkI8Cbgrhab32Jm95nZt8zsD3taGDjwHTPbZ2bjLbbn5RheQfv/QLI8fgCvcPfHIPzPGXh5i33ychz/gvAvpFZW+i6k6Zp618zn23Qd5eH4/TFwxN1/3GZ7lscvliyC2lqsax4jGGef1JnZC4GvAde6+5NNm+8h/HP+fOCfgf/scXkXu/sFwDuBq83sbU3bMz+GZnYq8B7gKy02Z3384srDcdwOPAtU2+yy0nchLRXg1cAbgccI3QvNMj9+wFaWP5vO6vjFlkVQHwJetejz2cDhVeyTKjM7hRDSVXe/tXm7uz/p7r+uv/8mcIqZbehVfe5+uL48CtxG+CfmYpkfQ8IX/x53P9K8IevjV3ek0R1UXx5tsU+mx9HM3g+8GxjzeodqsxjfhVS4+xF3/527HwP+tc3vzfr4nQz8GfDldvtkdfw6kUVQzwCvMbNz6mdcVwC3N+1zO3BVfeTCm4EnGv9E7YV6n9bngP3u/qk2+5xR3w8zu4hwLH/Zo/peYGYvarwnXHR6oGm3TI9hXdszmSyP3yK3A++vv38/8PUW+8T5vqbCzC4D/gZ4j7svtNknznchrfoWX/O4vM3vzez41W0GHnL3Q602Znn8OpLFFUzCiIQfEa4Gb6+v+xDwofp7A26sb/8BUOhxfX9E+OfZ/cC99de7mmq8BniQcBX7TuCtPazv3Prvva9eQx6P4SAheF+8aF1mx4/wP4zHgN8SzvI+ALwM+G/gx/XlS+v7vhL45nLf1x7V9wihf7fxHfyX5vrafRd6VN+/179b9xPC98w8Hb/6+i80vnOL9u358Uv60i3kIiI5pzsTRURyTkEtIpJzCmoRkZxTUIuI5JyCWkQk5xTUIiI5p6AWEcm5/wdc2M8+f5nbLgAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "x = cp.Variable(20)\n", "\n", "objective = cp.Minimize(sum(cp.huber(A@x-b,1)))\n", "\n", "constraints = [x[1:]>=x[:-1]]\n", "\n", "p_star = cp.Problem(objective, constraints).solve(solver='ECOS')\n", "\n", "np.set_printoptions(formatter={'float': lambda x: \"{0:0.2f}\".format(x)})\n", "print(p_star)\n", "print(x.value)\n", "\n", "plt.plot(np.arange(20),x_true, color=\"red\", linestyle=\"\", marker='o')\n", "plt.plot(np.arange(20),x.value, color=\"blue\", linestyle=\"\", marker='x')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "DYS iteration" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Optimal value from CVXPY: 24.733210916870362\n", "Optimal value from DYS: 24.733210919782255\n", "\n", "x^* from CVXPY: [-0.28 -0.28 0.92 2.14 3.39 3.39 3.39 3.73 3.73 3.73 5.29 5.52 5.52 6.08\n", " 6.61 7.35 7.35 7.35 7.79 7.79]\n", "x^k from DYS: [-0.28 -0.28 0.92 2.14 3.39 3.39 3.39 3.73 3.73 3.73 5.29 5.52 5.52 6.08\n", " 6.61 7.35 7.35 7.35 7.79 7.79]\n" ] } ], "source": [ "# Computes projection onto the set {x_{i+1}-x_i >= 0, i even}\n", "def even_proj(y):\n", " proj = y.copy()\n", " for i in range(len(y)//2):\n", " if proj[2*i]>proj[2*i+1]:\n", " proj[2*i] = (proj[2*i]+proj[2*i+1])/2\n", " proj[2*i+1] = proj[2*i]\n", " return proj\n", "\n", "\n", "# Computes projection onto the set {x_{i+1}-x_i >= 0, i odd}\n", "def odd_proj(y):\n", " proj = y.copy()\n", " for i in range((len(y)-1)//2):\n", " if proj[2*i+1]>proj[2*i+2]:\n", " proj[2*i+1] = (proj[2*i+1]+proj[2*i+2])/2\n", " proj[2*i+2] = proj[2*i+1]\n", " return proj\n", "\n", "\n", "\n", "zk = np.zeros(20)\n", "alpha = 0.01\n", "\n", "# DYS\n", "for k in range(5000):\n", " xkhalf = odd_proj(zk)\n", " xk1 = even_proj(2*xkhalf-zk-alpha*A.T@huber_grad(A@xkhalf-b))\n", " zk = zk+xk1-xkhalf\n", " \n", "\n", "print('Optimal value from CVXPY:', p_star)\n", "print('Optimal value from DYS: ', huber_loss(A@xk1-b))\n", "\n", "print('\\nx^* from CVXPY: ',x.value)\n", "print('x^k from DYS: ',xk1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "PAPC iteration" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Optimal value from CVXPY: 24.733210916870362\n", "Optimal value from PAPC: 24.733210919782262\n", "\n", "x^* from CVXPY: [-0.28 -0.28 0.92 2.14 3.39 3.39 3.39 3.73 3.73 3.73 5.29 5.52 5.52 6.08\n", " 6.61 7.35 7.35 7.35 7.79 7.79]\n", "x^k from PAPC: [-0.28 -0.28 0.92 2.14 3.39 3.39 3.39 3.73 3.73 3.73 5.29 5.52 5.52 6.08\n", " 6.61 7.35 7.35 7.35 7.79 7.79]\n" ] } ], "source": [ "# Define the matrix D\n", "from scipy.linalg import circulant\n", "\n", "D = circulant(np.concatenate(([1,-1],np.zeros(18))))[1:,:]\n", "\n", "uk = np.zeros(19)\n", "xk = np.zeros(20)\n", "alpha, beta = 0.01, 10\n", "\n", "# PAPC\n", "for k in range(5000):\n", " uk = np.minimum(uk+beta*D@(xk-alpha*D.T@uk-alpha*A.T@huber_grad(A@xk-b)),0)\n", " xk = xk - alpha*(np.matmul(D.T,uk) + A.T@huber_grad(A@xk-b))\n", " \n", "print('Optimal value from CVXPY:', p_star)\n", "print('Optimal value from PAPC: ', huber_loss(A@xk-b))\n", "\n", "print('\\nx^* from CVXPY: ',x.value)\n", "print('x^k from PAPC: ',xk)" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.3" } }, "nbformat": 4, "nbformat_minor": 4 }