pset 5 solutions
This commit is contained in:
Steven G. Johnson committed
1 parent
da5befa30b
commit
514644194f
2 files changed
+783
-1
No files matched your search
@@ -332,7 +332,7 @@ http://dx.doi.org/10.1137/S1052623499362822) — I used the "linear and separabl
|
||||
## Lecture 29 (Apr 21)
|
||||
|
||||
* [(Discrete) convolutions](https://en.wikipedia.org/wiki/Convolution#Discrete_convolution), translation-invariance, [circulant matrices](https://en.wikipedia.org/wiki/Circulant_matrix), and convolutional neural networks (CNNs)
|
||||
* pset 5 solutions:
|
||||
* [pset 5 solutions](psets/pset5sol.ipynb)
|
||||
* pset 6: coming soon, due Friday May 5.
|
||||
|
||||
**Further reading:** Strang textbook sections IV.2 (circulant matrices) and VII.2 (CNNs), and [OCW lecture 32](https://ocw.mit.edu/courses/18-065-matrix-methods-in-data-analysis-signal-processing-and-machine-learning-spring-2018/resources/lecture-32-imagenet-is-a-cnn-the-convolution-rule/). See also these [Stanford lecture slides](http://cs231n.stanford.edu/slides/2016/winter1516_lecture7.pdf) and [MIT lecture slides](https://mit6874.github.io/assets/sp2020/slides/L03_CNNs_MK2.pdf).
|
||||
@@ -0,0 +1,782 @@
|
||||
{
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "701e5077",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"# 18.065 Problem Set 5 Solutions\n",
|
||||
"\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "74f502b8",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Problem 1 (5+6 points)\n",
|
||||
"\n",
|
||||
"Consider the following optimization problem:\n",
|
||||
"$$\n",
|
||||
"\\min_{x \\in \\mathbb{R}^2} x_1 \\\\\n",
|
||||
"\\mbox{ subject to } x_2 \\le x_1^3 \\mbox{ and } x_2 \\ge 0 \\, .\n",
|
||||
"$$\n",
|
||||
"\n",
|
||||
"**(a)** Draw a sketch of the feasible set in the $x_1, x_2$ plane and indicate the optimum $x_*$.\n",
|
||||
"\n",
|
||||
"**(b)** Show that the optimum $x_*$ does *not* satisfy the KKT conditions, but explain why this is possible because the LICQ conditions are violated (see the last slide of lecture 22).\n",
|
||||
"\n",
|
||||
"(Most problems have local minima that satisfy KKT, but you can see from the picture in (a) that this is a weird case!)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "0075a58b",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"### Solution:\n",
|
||||
"\n",
|
||||
"**(a)** We will use PyPlot to draw our feasible set as a light-blue shaded region, which lies below the cubic $x_2 = x_1^3$ and above the $x_2=0$ axis, and hence is a **cusp-like** region that **only includes** $x_1 \\ge 0$.\n",
|
||||
"\n",
|
||||
"We also the primal optimum $x_*$ as a red star, which lies at $\\boxed{x_1 = 0, x_2 = 0}$, since this is the smallest $x_1$ point in the feasible region."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 1,
|
||||
"id": "f147b166",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"ename": "LoadError",
|
||||
"evalue": "UndefVarError: plot not defined",
|
||||
"output_type": "error",
|
||||
"traceback": [
|
||||
"UndefVarError: plot not defined",
|
||||
"",
|
||||
"Stacktrace:",
|
||||
" [1] top-level scope",
|
||||
" @ In[1]:3"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"# plot constraints\n",
|
||||
"x₁ = range(-2,2,length=200)\n",
|
||||
"plot(x₁, x₁.^3, \"b-\", linewidth=3)\n",
|
||||
"text(1, 1.5, \"x₂ ≤ x₁³\", color=\"b\", rotation=45)\n",
|
||||
"plot(x₁, zero(x₁), \"c-\", linewidth=3)\n",
|
||||
"text(1.5, 0.4, \"x₂ ≥ 0\", color=\"c\")\n",
|
||||
"\n",
|
||||
"# plot feasible set\n",
|
||||
"x₁⁺ = range(0,2,length=100)\n",
|
||||
"fill_between(x₁⁺, zero(x₁⁺), x₁⁺.^3, color=(0.8,0.8,1))\n",
|
||||
"text(1.4, 1.8, \"feasible\")\n",
|
||||
"\n",
|
||||
"# plot optimum\n",
|
||||
"plot([0],[0],\"r*\",markersize=15)\n",
|
||||
"text(0.1,0.4,L\"x_*\",color=\"r\")\n",
|
||||
"\n",
|
||||
"title(\"problem 1(a): feasible set and primal optimum\")\n",
|
||||
"\n",
|
||||
"# center the axes:\n",
|
||||
"ax = gca()\n",
|
||||
"ax.spines[\"left\"].set_position(\"zero\")\n",
|
||||
"ax.spines[\"right\"].set_color(\"none\")\n",
|
||||
"ax.spines[\"bottom\"].set_position(\"zero\")\n",
|
||||
"ax.spines[\"top\"].set_color(\"none\")"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "71a039ef",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"**(b)** Our objective is $f_0(x) = x_1$, and our two constraints $f_i(x) \\le 0$ are $f_1(x) = x_2 - x_1^3$ and $f_2(x) = -x_2$. Therefore, at the optimum $x_* = (0,0)$, their gradients are:\n",
|
||||
"$$\n",
|
||||
"\\left. \\nabla f_0 \\right|_{x_*} = (1,0), \\qquad\n",
|
||||
"\\left. \\nabla f_1 \\right|_{x_*} = (0,1), \\qquad\n",
|
||||
"\\left. \\nabla f_2 \\right|_{x_*} = (0,-1) \\, .\n",
|
||||
"$$\n",
|
||||
"These are depicted graphically on the diagram below.\n",
|
||||
"\n",
|
||||
"But these **cannot possibly** satisfy the KKT condition $\\boxed{\\nabla L = 0}$ (for the Lagrangian $L = f_0 + \\lamba_1 f_1 + \\lambda_2 f_2$), for **any** Lagrange multipliers $\\lambda_1, \\lambda_2$, because $\\nabla f_0$ is **linearly independent of** (*orthogonal* to!) $\\nabla f_1$ and $\\nabla f_2$!\n",
|
||||
"\n",
|
||||
"It is possible to violate the KKT conditions at $x_*$ in this problem because the optimum **violates LICQ**: the constraint gradients $\\nabla f_1$ and $\\nabla f_2$ are linearly dependent (anti-parallel!) at $x_*$.\n",
|
||||
"\n",
|
||||
"This happens because the feasible set **has a cusp** (where two constraints become parallel) right at the primal optimum. Obviously, this is an odd case that doesn't arise too often in practice — most optimization problems typically satisfy LICQ and KKT conditions at any local optimum."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 2,
|
||||
"id": "344ade48",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"ename": "LoadError",
|
||||
"evalue": "UndefVarError: plot not defined",
|
||||
"output_type": "error",
|
||||
"traceback": [
|
||||
"UndefVarError: plot not defined",
|
||||
"",
|
||||
"Stacktrace:",
|
||||
" [1] top-level scope",
|
||||
" @ In[2]:3"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"# plot constraints\n",
|
||||
"x₁ = range(-2,2,length=200)\n",
|
||||
"plot(x₁, x₁.^3, \"b-\", linewidth=3)\n",
|
||||
"text(1, 1.5, \"x₂ ≤ x₁³\", color=\"b\", rotation=45)\n",
|
||||
"plot(x₁, zero(x₁), \"c-\", linewidth=3)\n",
|
||||
"text(1.5, 0.4, \"x₂ ≥ 0\", color=\"c\")\n",
|
||||
"\n",
|
||||
"# plot feasible set\n",
|
||||
"x₁⁺ = range(0,2,length=100)\n",
|
||||
"fill_between(x₁⁺, zero(x₁⁺), x₁⁺.^3, color=(0.8,0.8,1))\n",
|
||||
"text(1.4, 1.8, \"feasible\")\n",
|
||||
"\n",
|
||||
"# plot optimum\n",
|
||||
"plot([0],[0],\"r*\",markersize=15)\n",
|
||||
"text(0.1,0.4,L\"x_*\",color=\"r\")\n",
|
||||
"\n",
|
||||
"title(L\"problem 1(b): gradients $\\nabla f_0, \\nabla f_1, \\nabla f_2$ at $x_*$\")\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"# gradients\n",
|
||||
"arrow(0,0,0,5,color=\"b\", width=0.04, head_width=0.2, head_length=0.8)\n",
|
||||
"text(0.1, 5, L\"\\nabla f_1 = \\nabla (x_2 - x_1^3)\", color=\"b\", fontsize=14)\n",
|
||||
"\n",
|
||||
"arrow(0,0,0,-5,color=\"c\", width=0.04, head_width=0.2, head_length=0.8)\n",
|
||||
"text(0.1, -5, L\"\\nabla f_2 = \\nabla (-x_2)\", color=\"c\", fontsize=14)\n",
|
||||
"\n",
|
||||
"arrow(0,0,1,0,color=\"r\", width=0.2, head_width=0.8, head_length=0.2, zorder=10)\n",
|
||||
"text(1, -2, L\"\\nabla f_0 = \\nabla (x_1)\", color=\"r\", fontsize=14)\n",
|
||||
"\n",
|
||||
"# center the axes:\n",
|
||||
"ax = gca()\n",
|
||||
"ax.spines[\"left\"].set_position(\"zero\")\n",
|
||||
"ax.spines[\"right\"].set_color(\"none\")\n",
|
||||
"ax.spines[\"bottom\"].set_position(\"zero\")\n",
|
||||
"ax.spines[\"top\"].set_color(\"none\")"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "3a49c4b1",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Problem 2 (6+6+6 points)\n",
|
||||
"\n",
|
||||
"Consider the convex problem:\n",
|
||||
"$$\n",
|
||||
"\\min_{x\\in \\mathbb{R}^n} \\Vert b - Ax \\Vert_2^2 \\\\\n",
|
||||
"\\mbox{ subject to } \\Vert x \\Vert_2^2 \\le r^2\n",
|
||||
"$$\n",
|
||||
"for some $r > 0$, $m \\times n$ matrix $A$ (of rank $n$), and $b \\in \\mathbb{R}^m$ — that is, least-squares optimization with the solution constrained to lie inside a sphere of radius $r$.\n",
|
||||
"\n",
|
||||
"**(a)** What is the Lagrange dual function $g(\\lambda)$? (You can give a closed-form expression. Hint: review Tikhonov-regularized least-squares.) Define a corresponding Julia function `g(λ; r=1.0)` for the sample parameters given below (this syntax defines an optional keyword argument `r` that defaults to $r=1$). Make a plot of $g(\\lambda)$ for $r=1$ and $r=0.5$ for $\\lambda \\ge 0$ to verify that it looks concave with a single maximum.\n",
|
||||
"\n",
|
||||
"**(b)** If the unconstrained least-square solution $\\hat{x} = (A^T A)^{-1} A^T b$ satisfies $\\Vert \\hat{x} \\Vert_2 < r$, then what must be true of the derivative $g'(0)$? What if $\\Vert \\hat{x} \\Vert_2 > r$?\n",
|
||||
"\n",
|
||||
"Check in Julia that $g'(0)$ matches your expectations by computing the derivative using automatic differentiation:\n",
|
||||
"```jl\n",
|
||||
"using ForwardDiff\n",
|
||||
"dgdλ(λ; r=1.0) = ForwardDiff.derivative(λ -> g(λ; r=r), λ)\n",
|
||||
"```\n",
|
||||
"and evaluating it at `dgdλ(0; r=???)` for two values of `r`: one $r > \\Vert \\hat{x} \\Vert_2$ (so that the constraint is inactive) and one $r < \\Vert \\hat{x} \\Vert_2$ (so that the constraint is active).\n",
|
||||
"\n",
|
||||
"(In principle, you could take this derivative by hand using matrix calculus, but it's pretty error-prone.)\n",
|
||||
"\n",
|
||||
"**(c)** You can take the *second* derivative of $g(\\lambda)$ via AD by:\n",
|
||||
"```jl\n",
|
||||
"d²gdλ²(λ; r=0.5) = ForwardDiff.derivative(λ -> dgdλ(λ; r=0.5), λ)\n",
|
||||
"```\n",
|
||||
"Use this to implement a Newton iteration to maximize $g(\\lambda)$ (for $\\lambda \\ge 0$) by finding a root of $g'(\\lambda)$, starting with an initial guess of $\\lambda=0$, for $r = 0.5$. (It should converge in only a few iterations. The solution should have $\\lambda > 0$ in this case because ...?) To at least 8 significant digits, give the resulting dual optimum $\\lambda_*$ and the primal optimum $x_*$ (strong duality holds in this convex problem!), and check that $x_*$ is feasible."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "15c6ac1c",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"### Solution:\n",
|
||||
"\n",
|
||||
"**(a)** The dual function $g(\\lambda)$ is the found by minimizing the Lagrangian over $x$, i.e.\n",
|
||||
"$$\n",
|
||||
"g(\\lambda) = \\min_x \\underbrace{\\left[ \\Vert b - Ax \\Vert_2^2 + \\lambda (\\Vert x \\Vert_2^2 - r^2) \\right]}_{L(x,\\lambda)}\n",
|
||||
"$$\n",
|
||||
"where $\\Vert x \\Vert_2^2 - r^2 \\le 0$ is our inequality constraint. But this is exactly in the form of Tikhonov-regularized least-squares (plus a constant $-\\lambda r^2$), so we know from class that the solution is\n",
|
||||
"$$\n",
|
||||
"\\hat{x}_\\lambda = (A^T A + \\lambda I)^{-1} (A^T b)\n",
|
||||
"$$\n",
|
||||
"and hence $g(\\lambda) = \\boxed{L(\\hat{x}_\\lambda, \\lambda)}$ in terms of the formulas above. We can implement this as follows:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 3,
|
||||
"id": "1a89466c",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"using LinearAlgebra\n",
|
||||
"\n",
|
||||
"m = 5\n",
|
||||
"n = 4\n",
|
||||
"A = [ -9 2 -2 3\n",
|
||||
" -5 -3 9 3\n",
|
||||
" -1 -6 9 -2\n",
|
||||
" -3 -4 5 4\n",
|
||||
" -8 9 -6 4 ]\n",
|
||||
"b = [1,2,3,4,5];"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 4,
|
||||
"id": "e5a615e1",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"g (generic function with 1 method)"
|
||||
]
|
||||
},
|
||||
"execution_count": 4,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"function g(λ; r=1)\n",
|
||||
" x̂_λ = (A'A + λ*I) \\ (A'b)\n",
|
||||
" return norm(b - A*x̂_λ)^2 + λ*(norm(x̂_λ)^2 - r^2)\n",
|
||||
"end"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "4001e6f4",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Plotting it, we obtain:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 5,
|
||||
"id": "521ec11c",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAjkAAAHICAYAAACszftMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAA9hAAAPYQGoP6dpAABSBElEQVR4nO3de3zO9f/H8cdlJzbbMDuY44ivs5wlcigikSRFRA7fr6LIt/NREtW3lFIqHSRERCGp5SzHZHIokTmEOdscx+zz++P927jaxjbbPtfheb/dPrdd1+f67LpefcKeex8dlmVZiIiIiHiYQnYXICIiIpIfFHJERETEIynkiIiIiEdSyBERERGPpJAjIiIiHkkhR0RERDySQo6IiIh4JIUcERER8UgKOSIiIuKRFHJERETEIynkiIiIiEdSyBFxMxMnTsThcKQfvr6+lClThgceeIB9+/bl6WcNHz4ch8PBkSNHslXTrl278vTzr8WiRYvo27cvVatWJSgoiNKlS3PHHXewfv36TK8fMWIE1atXJzU1Ncef9cknn1C6dGlOnz59rWXnSGY1nzt3Dl9fX4KDg3n44Ycz/T676hUpaAo5Im7qs88+Y9WqVcTGxjJgwAC+/PJLmjdvrh9c/2/8+PHs2rWLIUOGMH/+fMaOHcuhQ4do0qQJixYtcrp2//79vP7664wYMYJChXL+z2Lv3r0JCgri9ddfz6vyryqrmh0OB4sWLaJRo0aMGzeOP//80yXqFbGFJSJu5bPPPrMAa926dU7nn3/+eQuwJk+enOX3nj59Okef9eKLL1qAdfjw4WzVFB8fn6P3z08HDx7McO7kyZNWZGSkdfPNNzudf+KJJ6zSpUtbFy9ezPXnvfHGG1ZoaGiO73FuXa3mBQsWWIA1derUTF8v6HpF7KCWHBEP0aRJEwB2794NXOpq+vXXX+natSvFixenUqVKAKxYsYKbb76Z4OBgAgMDadq0Kd99912W77137166dOlCSEgIoaGh9OzZk8OHD1+1pu3bt9OjRw8iIiIICAigWrVqvPfee07XpNX522+/cffddxMaGkqJEiUYNmwYKSkpbNu2jXbt2hEcHEyFChWy3foQERGR4VzRokWpXr06e/fuTT93/vx5PvnkE3r06JGhFWfHjh088MADVK5cmcDAQEqXLk3Hjh3ZtGlThve+7777SEpKYtq0admqLzOWZTFu3DiqV69O4cKFqVmzJvPnz6dHjx7UqFEjWzWnqVq1KgBxcXGZvp4X9Yq4OoUcEQ+xY8cOAMLDw53Od+nSheuuu44ZM2bwwQcfsHTpUlq3bk1iYiKffPIJX375JcHBwXTs2JHp06dn+t533nkn1113HTNnzmT48OF888033HrrrVy4cCHLerZu3UrDhg3ZvHkzb775JvPmzaNDhw488sgjvPTSSxmu79atG3Xq1OHrr79mwIABvPXWWzz66KN07tyZDh06MHv2bFq3bs2TTz7JrFmzcnWPEhMT+fXXX50Cw5o1azh69CitWrXKcP3+/fsJCwvj1VdfZcGCBbz33nv4+vrSuHFjtm3b5nRtVFQUVatWzRAWHQ4HLVu2zFZ9Dz74IMOGDaNTp07MmTOHhx56iD59+rBs2TIaNmyYrZrTvPbaa0DWISerekU8it1NSSKSM2ldQ6tXr7YuXLhgnTx50po3b54VHh5uBQcHWwkJCZZlXepqeuGFF5y+v0mTJlZERIR18uTJ9HMpKSlWzZo1rTJlylipqanp59Pe49FHH3V6jylTpjh1jWXWXXXrrbdaZcqUsRITE52+d/DgwVbhwoWtY8eOOX3Gm2++6XTd9ddfbwHWrFmz0s9duHDBCg8Pt7p06ZLT22ZZlmXdd999lq+vr/XLL7+kn3vttdcsIP2+XUlKSop1/vx5q3LlyhnuSdr7R0ZGOp3z8fGxWrdufdX3njFjhgVYn3/+udP54cOHW4D17rvvZrvmH374wQKs4sWLW+Hh4Vl+Zmb1ingSteSIuKkmTZrg5+dHcHAwt99+O1FRUXz//fdERkY6XXfXXXelPz59+jRr1qyha9euFC1aNP28j48PvXr14u+//87QQgGma+Ny3bp1w9fXl8WLF2da27lz51i4cCF33nkngYGBpKSkpB+33XYb586dY/Xq1U7fc/vttzs9r1atGg6Hg/bt26ef8/X15brrrkvvksuJ559/nilTpvDWW29Rv3799PP79+/H4XBQsmTJDN+TkpLCqFGjqF69Ov7+/vj6+uLv78/27dv5/fffM1wfERHBoUOHSElJcXqPhQsXXrW+V155hRtvvJH777/f6XyVKlUAnFpyrlTziRMn6Nu3L3fccQcPPfQQhw8fZv/+/Zl+Zmb1ingSX7sLEJHcmTRpEtWqVcPX15fIyEhKlSqV6XWXnz9+/DiWZWV6bXR0NABHjx7N8FpUVJTTc19fX8LCwjK9Nu09UlJSePfdd3n33Xczveaf09JLlCjh9Nzf35/AwEAKFy6c4XxSUlKm75mVl156iZEjR/LKK68wePBgp9fOnj2Ln58fPj4+Gb5v2LBhvPfeezz55JO0aNGC4sWLU6hQIfr378/Zs2czXF+4cGEsy+LcuXNOIfJqEhISiIuL4+23387w2t69e/H19aVOnTrZqnnQoEFcuHCBCRMmsGzZMsB0WaX9/82LekXchUKOiJuqVq0aDRo0uOp1Docj/XHaD+kDBw5kuC7tt/3MWgcSEhIoXbp0+vOUlBSOHj1KWFhYpp9ZvHjx9NahQYMGZXpNTEzMVWvPCy+99BLDhw9n+PDhPPPMMxleL1myJOfPn+f06dMEBQU5vTZ58mTuv/9+Ro0a5XT+yJEjFCtWLMN7HTt2jICAgBwHhr///hsg0/D5448/UqtWLaewl1XNM2fOZOrUqcydO5fw8HDq1asHmJBz22235Vm9Iu5C3VUiXiQoKIjGjRsza9Ysp5aI1NRUJk+eTJkyZdK7Ry43ZcoUp+dfffUVKSkpWQ6oDQwMpFWrVmzYsIHatWvToEGDDEdWASkvvfzyywwfPpznnnuOF198MdNr0mYh/fXXXxleczgcBAQEOJ377rvvslx0cefOnVSvXj3HdabNBNu6davT+e+//56FCxdmCLOZ1ZyQkMDAgQP597//nd71FxMTQ7FixdiwYUOe1iviLtSSI+JlRo8eTZs2bWjVqhWPPfYY/v7+vP/++2zevJkvv/zSqeUnzaxZs/D19aVNmzZs2bKF559/njp16tCtW7csP2fs2LE0a9aM5s2b8+CDD1KhQgVOnjzJjh07mDt3boYF+fLam2++yQsvvEC7du3o0KFDhjFAaVPu04La6tWrqV27ttM1t99+OxMnTqRq1arUrl2b9evX87///Y8yZcpk+LzU1FTWrl1Lv379nM77+vrSokWLK47LKVu2LI0aNWLMmDGEh4dTo0YNfv75ZyZMmABAo0aNnK7PrOYBAwZQvHhxxowZ43Rt3bp1M51hlVW9Ih7F5oHPIpJDWS0G+E9XWshv+fLlVuvWra2goCCrSJEiVpMmTay5c+dm+R7r16+3OnbsaBUtWtQKDg62unfv7rTYXlaLAcbHx1t9+/a1Spcubfn5+Vnh4eFW06ZNrZEjR161zt69e1tBQUEZamrRooVVo0aNK/63p10HZHlcrnnz5tZtt92W4T2OHz9u9evXz4qIiLACAwOtZs2aWcuXL7datGhhtWjRwunahQsXpt+rywEZrs1MfHy81a5dO6to0aJWeHi41adPH+vTTz+1AGvnzp0Zrr+85gkTJlg+Pj7WqlWrMlw3bNgwy+FwWElJSdmqV8STOCzLsgo+WomIuI6vv/6ae+65h927dzuNPcqJXr16sXPnTn7++ec8qens2bPcfPPNhISEsGDBggyvX2vNeV2viCtSyBERr2dZFk2bNqV+/fqMGzcux9//119/Ua1aNRYtWkSzZs1y/P2//PILM2fO5MYbbyQ4OJg///yTd999lyNHjvDzzz9TsWLFPK35WusVcRcaeCwiXs/hcDBhwgSio6NztQv5nj17GDduXK4Dw6lTp4iNjaVnz560bduWUaNG0bx5c3799ddMA8611nyt9Yq4C7XkiIiIiEdyyZac0aNH07BhQ4KDg4mIiKBz584ZVmG1LIvhw4cTHR1NkSJFaNmyJVu2bLGpYhEREXE1Lhlyli5dyqBBg1i9ejWxsbGkpKTQtm1bTp8+nX7N66+/zpgxYxg3bhzr1q0jKiqKNm3acPLkSRsrFxEREVfhFt1Vhw8fJiIigqVLl3LTTTdhWRbR0dEMHTqUJ598EoDk5GQiIyN57bXX+M9//mNzxSIiImI3t1gMMDExEbi0t018fDwJCQm0bds2/ZqAgABatGjBypUrsww5ycnJJCcnpz9PTU3l2LFjhIWFZboAmoiIiLgey7I4efIk0dHRFCqUdaeUy4ccy7IYNmwYzZo1o2bNmoBZvhzIsNtyZGTkFXcnHj16NC+99FL+FSsiIiIFZu/evZmuQJ7G5UPO4MGD+e2331ixYkWG1/7Z+mJZ1hVbZJ5++mmGDRuW/jwxMZFy5cqxd+9eQkJC8q5oERERyTdJSUmULVuW4ODgK17n0iHn4YcfZs6cOSxbtswpqUVFRQGmRefyXXsPHTqUoXXncgEBARk22wMICQlRyBEREXEzVxtq4pKzqyzLYvDgwcyaNYtFixYRExPj9HpMTAxRUVHExsamnzt//jxLly6ladOmBV2uiIiIuCCXbMkZNGgQU6dO5dtvvyU4ODh9DE5oaChFihTB4XAwdOhQRo0aReXKlalcuTKjRo0iMDCQHj162Fy9iIiIuAKXDDnjx48HoGXLlk7nP/vsM/r06QPAE088wdmzZ3nooYc4fvw4jRs35scff7xq/5yIiIh4B7dYJye/JCUlERoaSmJiosbkiIhIli5evMiFCxfsLsNr+Pn54ePjk+Xr2f357ZItOSIiIq7AsiwSEhILine truncated
|
||||
"text/plain": [
|
||||
"Figure(PyObject <Figure size 640x480 with 1 Axes>)"
|
||||
]
|
||||
},
|
||||
"metadata": {},
|
||||
"output_type": "display_data"
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"(-5.0, 20.0)"
|
||||
]
|
||||
},
|
||||
"execution_count": 5,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"using PyPlot\n",
|
||||
"λ = range(0,100, length=1000)\n",
|
||||
"plot(λ, g.(λ, r=0.5), \"r-\")\n",
|
||||
"plot(λ, g.(λ, r=1), \"b--\")\n",
|
||||
"title(L\"Problem 2(a): $g(\\lambda)$\")\n",
|
||||
"xlabel(L\"\\lambda\")\n",
|
||||
"ylabel(L\"g(\\lambda)\")\n",
|
||||
"legend([L\"r=0.5\", L\"r=1\"])\n",
|
||||
"ylim(-5, 20)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "864bda6b",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"**(b)** If the unconstrained least-square solution satisfies the constraint, then the constraint is **inactive** and we expect that the corresponding Lagrange multiplier λ should be **zero** (corresponding to the \"complementary slackness\" KKT condition). In the dual problem, this corresponds to $$g(\\lambda)$$ for $$\\lambda \\ge 0$$ having its maximum at $\\lambda = 0$, corresponding to $\\boxed{g'(0) < 0}$ for $\\Vert \\hat{x} \\Vert_2 < r$.\n",
|
||||
"\n",
|
||||
"Conversely, if $\\Vert \\hat{x} \\Vert_2 > r$, then the constraint is **active** and we expect a dual maximum for $\\lambda > 0$, which requires $\\boxed{g'(0) > 0}$ in that case.\n",
|
||||
"\n",
|
||||
"To get an example of this, let's check the norm of the unconstrained least-square solution:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 6,
|
||||
"id": "b7facd32",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"0.6508599263498679"
|
||||
]
|
||||
},
|
||||
"execution_count": 6,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"norm(A \\ b)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "32751a1e",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"i.e. $\\Vert \\hat{x} \\Vert_2 = \\Vert (A^T A)^{-1} A^T b \\Vert_2 \\approx 0.651$, so the constraint should be inactive for any radius larger than that.\n",
|
||||
"\n",
|
||||
"For example, the constraint should be *inactive* for $r=1$ and *active* for $r = 0.5$, and we can see by inspection of the plot from part (a) that $g'(0)$ is $< 0$ and $> 0$ in these two cases, respectively. Quantitatively, we can use AD to get the derivative as suggested:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 7,
|
||||
"id": "e3a11067",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"dgdλ(0; r = 1.0) = -0.5763813562718445\n",
|
||||
"dgdλ(0; r = 0.5) = 0.17361864372815555\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"0.17361864372815555"
|
||||
]
|
||||
},
|
||||
"execution_count": 7,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"using ForwardDiff\n",
|
||||
"dgdλ(λ; r=1.0) = ForwardDiff.derivative(λ -> g(λ; r=r), λ)\n",
|
||||
"\n",
|
||||
"@show dgdλ(0; r=1.0)\n",
|
||||
"@show dgdλ(0; r=0.5)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "4bf7db46",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Hurray, the signs match our expectations!"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "41241ac1",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"**(c)** As suggested, we can get a second derivative by:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 8,
|
||||
"id": "ca12a000",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"d²gdλ² (generic function with 1 method)"
|
||||
]
|
||||
},
|
||||
"execution_count": 8,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"d²gdλ²(λ; r=0.5) = ForwardDiff.derivative(λ -> dgdλ(λ; r=0.5), λ)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "484a2d51",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"We can use this to implement a Newton iteration for the dual optimum. A Newton iteration to find a maximum $g'(\\lambda) = 0$ is just\n",
|
||||
"$$\n",
|
||||
"\\lambda \\longleftarrow \\lambda - \\frac{g'(\\lambda)}{g''(\\lambda)}\n",
|
||||
"$$\n",
|
||||
"We can start at $\\lambda = 0$. Of course, we first check whether $g'(0) > 0$ to be sure the constraint is active; otherwise, the dual optimum is simply $\\lambda = 0$.\n",
|
||||
"\n",
|
||||
"Given the dual optimum $\\lambda_*$, by strong duality we then have the corresponding primal optimum $x_* = \\hat{x}_{\\lambda_*}$ from part (a)."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 9,
|
||||
"id": "0d7df6ac",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"dualopt (generic function with 1 method)"
|
||||
]
|
||||
},
|
||||
"execution_count": 9,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"# find the dual optimum (maximizing g(λ))\n",
|
||||
"# and return the dual and primal optima λ,x\n",
|
||||
"function dualopt(r; rtol=10*eps(Float64), maxiter=100)\n",
|
||||
" λ = zero(rtol)\n",
|
||||
" g′ = dgdλ(λ; r=r)\n",
|
||||
" if g′ > 0 # active constraint\n",
|
||||
" for iter=1:maxiter # Newton iterations\n",
|
||||
" g″ = d²gdλ²(λ; r=r)\n",
|
||||
" @show Δλ = g′ / g″\n",
|
||||
" λ -= Δλ\n",
|
||||
" abs(Δλ) ≤ rtol * abs(λ) && break # converged\n",
|
||||
" g′ = dgdλ(λ; r=r)\n",
|
||||
" end\n",
|
||||
" end\n",
|
||||
" x̂_λ = (A'A + λ*I) \\ (A'b) # primal optimum x\n",
|
||||
" return λ, x̂_λ\n",
|
||||
"end"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "fff24d2b",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"For $r = 0.5$, this gives (to nearly machine precision in only 7 Newton steps!)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 10,
|
||||
"id": "4c37fc10",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"Δλ = g′ / g″ = -2.976939515696773\n",
|
||||
"Δλ = g′ / g″ = -2.197446051423321\n",
|
||||
"Δλ = g′ / g″ = -0.7226090155288597\n",
|
||||
"Δλ = g′ / g″ = -0.053830681767815296\n",
|
||||
"Δλ = g′ / g″ = -0.000263613901856008\n",
|
||||
"Δλ = g′ / g″ = -6.262750976944395e-9\n",
|
||||
"Δλ = g′ / g″ = -2.9522501729702895e-15\n",
|
||||
"(λ_opt, x_opt) = dualopt(0.5) = (5.951088884581378, [-0.2084478531413324, 0.23289743085848733, 0.29901616790443714, 0.25079396035796236])\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"(5.951088884581378, [-0.2084478531413324, 0.23289743085848733, 0.29901616790443714, 0.25079396035796236])"
|
||||
]
|
||||
},
|
||||
"execution_count": 10,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"@show λ_opt, x_opt = dualopt(0.5)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "2cd47afd",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"i.e. $\\boxed{\\lambda_* \\approx 5.951088884581378}$ and\n",
|
||||
"$\\boxed{x_* \\approx [-0.2084478531413324, 0.23289743085848733, 0.29901616790443714, 0.25079396035796236]}$."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "b4f2f1ef",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"We should also check that $x_*$ is feasible, and in fact that the constraint is active ($\\Vert x_* \\Vert_2 = r$) since $\\lambda_* > 0$:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 11,
|
||||
"id": "2c3d0137",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"0.5000000000000001"
|
||||
]
|
||||
},
|
||||
"execution_count": 11,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"norm(x_opt)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "530e3856",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Yes, it is feasible and the constraint is active (up to roundoff errors)."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "7878762b",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Problem 3 (5+5+6 points)\n",
|
||||
"\n",
|
||||
"In this problem, you will use ADMM to solve the (primal) optimization problem from problem 2 above, for the parameters from problem 2c, using the equivalent formulation:\n",
|
||||
"$$\n",
|
||||
"\\min_{x \\in \\mathbb{R}^n} \\left( \\Vert b - Ax \\Vert_2^2 + \\begin{cases} 0 & \\Vert x \\Vert_2 \\le r \\\\ \\infty & \\mbox{otherwise} \\end{cases} \\right)\n",
|
||||
"$$\n",
|
||||
"where the second term is the \"indicator\" function of the feasible set (the radius-$r$ ball) as in lecture 24 and section III.4 of the text.\n",
|
||||
"\n",
|
||||
"A basic iteration of ADMM consists of 3 steps, as described in class and in the textbook:\n",
|
||||
"\n",
|
||||
"1. $x^{(k+1)} = \\mbox{arg }\\min_x \\Vert b - Ax \\Vert_2^2 + \\frac{\\rho}{2} \\Vert x - z^{(k)} + s^{(k)} \\Vert_2^2$ (for some penalty parameter $\\rho > 0$)\n",
|
||||
"2. $z^{(k+1)}$ is the projection of $x^{(k+1)} + s^{(k)}$ onto the *closest* point in the feasible set.\n",
|
||||
"3. $s^{(k+1)} = s^{(k)} + x^{(k+1)} - z^{(k+1)}$\n",
|
||||
"\n",
|
||||
"**(a)** Give a closed-form solution for step 1. (Hint: a problem from pset 3 should be helpful.)\n",
|
||||
"\n",
|
||||
"**(b)** Give a closed-form solution for step 2.\n",
|
||||
"\n",
|
||||
"**(c)** Implement this iteration in Julia to solve this problem with the parameters from 2c above, starting from $x = z = s = \\vec{0}$. Make a (semi-log) plot of the error $\\Vert x^{(k)} - x_* \\Vert_2$ versus $k$, where $x_*$ is your solution from 2c, for $\\rho = 1$ and $\\rho = 10$. (The error should converge to zero!)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "ee954031",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"### Solution:\n",
|
||||
"\n",
|
||||
"**(a)** This is essentiallly the same as problem 3 of pset 3, with the \"prior\" $x_p = z^{(k)} - s^{(k)}$ and a penalty strength $\\delta^2 = \\rho/2$. Quoting the solutions from pset 3, we then have:\n",
|
||||
"$$\n",
|
||||
"\\boxed{x^{(k+1)} = \\left(A^T A + \\frac{\\rho}{2} I\\right)^{-1} \\left[A^T b + \\frac{\\rho}{2} (z^{(k)} - s^{(k)}) \\right]}\n",
|
||||
"$$\n",
|
||||
"\n",
|
||||
"**(b)** To find the closest point to a vector $v$ on a sphere of radius $r$, if $\\Vert v \\Vert_2 > r$ we simply shrink $v$ onto the sphere: $v \\to v r / \\Vert v \\Vert_2$. This gives:\n",
|
||||
"\n",
|
||||
"$$\n",
|
||||
"\\boxed{z^{(k+1)} = \\underbrace{\\left( x^{(k+1)} + s^{(k)} \\right)}_{v^{(k)}} \\frac{\\min \\{r, \\Vert v^{(k)} \\Vert_2\\}}{\\Vert v^{(k)} \\Vert_2} }\n",
|
||||
"$$\n",
|
||||
"where the second term shrinks the length if it is $> r$. (You can write this in other ways, of course, e.g. as an \"if–then\" statement.)\n",
|
||||
"\n",
|
||||
"**(c)** We'll implement a function `x,z,s = ADMM(x,z,s)` to take an ADMM step using the formulas above, and then iterate it for 1000 iterations and plot the convergence towards the optimum `x_opt` from problem 2(c) above:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 12,
|
||||
"id": "e0f4080f",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"ADMM (generic function with 1 method)"
|
||||
]
|
||||
},
|
||||
"execution_count": 12,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"# perform 1 ADMM step, returning the new x,z,s vectors\n",
|
||||
"function ADMM(x,z,s; r=0.5, ρ=1)\n",
|
||||
" x = (A'A + (ρ/2)*I) \\ (A'b + (ρ/2)*(z - s))\n",
|
||||
" v = x + s\n",
|
||||
" vlen = norm(v)\n",
|
||||
" z = vlen ≤ r ? v : v * r/vlen\n",
|
||||
" s = s + x - z\n",
|
||||
" return x,z,s\n",
|
||||
"end"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 13,
|
||||
"id": "483e87b6",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"x = z = s = zeros(n)\n",
|
||||
"err1 = [norm(x - x_opt)]\n",
|
||||
"for i = 2:1000\n",
|
||||
" x,z,s = ADMM(x,z,s; ρ=1)\n",
|
||||
" push!(err1, norm(x - x_opt))\n",
|
||||
"end\n",
|
||||
"\n",
|
||||
"x = z = s = zeros(n)\n",
|
||||
"err10 = [norm(x - x_opt)]\n",
|
||||
"for i = 2:1000\n",
|
||||
" x,z,s = ADMM(x,z,s; ρ=10)\n",
|
||||
" push!(err10, norm(x - x_opt))\n",
|
||||
"end"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 14,
|
||||
"id": "f2ce39f1",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAkwAAAHFCAYAAAAAM6ZOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAA9hAAAPYQGoP6dpAABfoElEQVR4nO3deVyU1f4H8M8IsoiAIoniIKJiYSjKQOYuaBCGhV3XXHAryczQunWN0quVVL+uVtfEXbTFNFMqr2WUuHsVEdRcMgtFFFRcQFxQ4fz+mDsjw8ywzDzM+nm/XvOqOc/Deb7zQPHlfM9zjkwIIUBEREREejUwdwBERERElo4JExEREVENmDARERER1YAJExEREVENmDARERER1YAJExEREVENmDARERER1YAJExEREVENmDARERER1YAJE1EVqampkMlk6pejoyPkcjnGjx+P8+fPS3otmUyGqVOn1nje9u3bIZPJsH37dkmvb4xPP/0Ujz/+OLy9veHs7IzWrVtjxIgROHbsmMliKCoqgrOzM2QyGQ4ePKjznHHjxml8P93c3NCmTRs8/fTTWLVqFcrKyrS+pl+/fpDJZGjbti10bYawc+dOdX+pqanq9so/O7q+V0IItG/fHjKZDP369TP0YxORGTBhItJj1apV2LdvH9LT0/H8889j7dq16N27N27evGnu0CzClStXEBMTg+XLl+Pnn3/GnDlzkJ2djW7duuH33383SQyff/457t69CwBYsWKF3vNcXV2xb98+7Nu3D5s3b8bcuXPh5uaG559/HgqFAvn5+Vpf4+7ujtzcXGzbtk3r2MqVK+Hh4aH3eu7u7jrj2bFjB/7880+4u7vX5uMRkQVhwkSkR3BwMB5//HFERERg9uzZeP3115Gbm4u0tDS9X3Pr1i3TBWhmc+bMwezZsxEXF4e+ffti/PjxSEtLw82bN/Hll1+aJIaVK1eiefPmCA8Px9q1a3H79m2d5zVo0ACPP/64+vs5duxYrF27Flu2bMGpU6cwZMgQra9p3bo1Hn/8caxcuVKj/caNG/jmm28wfPhwvXENHz4c3377LUpKSjTaV6xYge7du6N169YGfFrbIYTQ+70islRMmIhq6fHHHwcAnD17FoCy1NO4cWMcPXoUUVFRcHd3R//+/QEAV69exZQpU9CqVSs4OTmhbdu2SEpK0ln+AYAlS5agQ4cOcHZ2RseOHfH111/XKqaDBw/i6aefhpeXF1xcXNC1a1esX79e4xxVmWjbtm14/vnn0axZM3h4eGDs2LG4efMmCgsLMWzYMDRp0gQtW7bEa6+9hnv37hl0jx566CEAgKOjo0b79evX8eqrr6Jt27ZwdnZG8+bNMXDgQJw8edKg6wDA/v378dtvv2HMmDF4/vnnUVxcjG+//bZOfURFReH555/H/v37sXPnTq3jEyZMwMaNG3H9+nV1m+p7M2LECL39jhw5EgCwdu1adZsqvgkTJtQpxq+++grdu3dH48aN0bhxY3Tp0kVr9GrlypUICQmBi4sLvLy8MHjwYJw4cULjHNXP6+nTpzFw4EA0btwYfn5+ePXVV9U/l/fu3UPz5s0xZswYrTiuX78OV1dXzJgxQ91WUlKC1157DQEBAXByckKrVq2QmJioNQqrKj0vXrwYQUFBcHZ2xurVqwEAu3fvRvfu3eHi4oJWrVrh7bffxvLlyyGTyXDmzBmNftatW4fu3bvDzc0NjRs3RnR0NLKzs+v8OVXKysowd+5cBAUFwcXFBc2aNUNERAT27t2rPkcIgUWLFqFLly5wdXVF06ZNMWTIEPz111/VfdvIBjFhIqql06dPA3iQFADA3bt38fTTTyMyMhLfffcd5syZgzt37iAiIgJr1qzBjBkz8J///AejR4/Ghx9+iGeffVar3++//x6ffvop5s6diw0bNsDf3x8jR47Ehg0bqo0nIyMDPXv2xPXr17F48WJ899136NKlC4YPH64xr0Zl0qRJ8PT0xNdff4233noLX331FZ5//nk89dRTCAkJwYYNGxAfH49//etf+Pe//13r+1JeXo6ysjKcPHkSkyZNQvPmzTF+/Hj18Rs3bqBXr15YsmQJxo8fjx9++AGLFy9Ghw4dUFBQoD5PNdeo6i9JfVRJw4QJEzBixAg0atSo2rKcPk8//TQA6EyYRowYAQcHB43EZ8WKFRgyZEi1JTkPDw8MGTJEY3Rq7dq1aNCgQbUjU1XNmjULo0aNgq+vL1JTU7Fp0ybEx8erk3YASE5OxsSJE/Hoo49i48aN+OSTT3DkyBF0794df/zxh0Z/9+7dw9NPP43+/fvju+++w4QJE7BgwQJ88MEHAICGDRti9OjROkfH1q5dizt37qi/t7du3ULfvn2xevVqTJs2DT/++CPeeOMNpKam4umnn9aa+5WWloaUlBTMmjULW7duRe/evXHkyBE88cQTuHXrFlavXo3Fixfj0KFDeO+997Tuxbx58zBy5Eh07NgR69evx+eff44bN26gd+/eOH78eJ0+JwDcv38fMTExeOeddxAbG4tNmzYhNTUVPXr0QF5envq8yZMnIzExEQMGDEBaWhoWLVqEY8eOoUePHrh48WKtv5dkAwQRaVi1apUAIP773/+Ke/fuiRs3bojNmzeLhx56SLi7u4vCwkIhhBDx8fECgFi5cqXG1y9evFgAEOvXr9do/+CDDwQA8fPPP6vbAAhXV1d1n0IIcf/+ffHII4+I9u3bq9syMjIEAJGRkaFue+SRR0TXrl3FvXv3NK4TGxsrWrZsKcrLyzU+z8svv6xxXlxcnAAg5s+fr9HepUsXERoaWtvbJZydnQUAAUB06NBBHD9+XOP43LlzBQCRnp5ebT8TJkwQDg4O4syZMzVe8+bNm8LDw0M8/vjj6rb4+Hghk8nE6dOnNc6Nj48Xbm5uevs6ceKEACBefPFFdVvfvn3Fo48+qv76sLAwIYQQx44dEwDE9u3bRWZmpgAgVq1apf461b3OzMxUf89+++03IYQQ4eHhYty4cUIIIR599FHRt2/faj/jX3/9JRwcHMSoUaP0nnPt2jXh6uoqBg4cqNGel5cnnJ2dxXPPPadxH3T9XA4cOFA8/PDD6vdHjhwRAMTSpUs1znvssceEQqFQv09OThYNGjQQmZmZGudt2LBBABBbtmxRtwEQnp6e4urVqxrnDh06VLi5uYnLly+r28rLy0XHjh0FAJGbm6v+PI6Ojlo/wzdu3BAtWrQQw4YNq/PnXLNmjQAgli1bJvTZt2+fACD+9a9/abSfO3dOuLq6itdff13v15Lt4QgTkR6PP/44GjZsCHd3d8TGxqJFixb48ccf4ePjo3He3/72N43327Ztg5ubm9a8mHHjxgEAfv31V432/v37a/Tp4OCA4cOH4/Tp0zonIwPK0a6TJ09i1KhRAJR/LateAwcOREFBgdbE69jYWI33QUFBAICnnnpKq73yCEZN9u7di3379uGLL76Au7s7IiIiNJ6U+/HHH9GhQwcMGDCg2n5WrFiB+/fvw9/fv8Zrrl+/HiUlJRrlrQkTJkAIgVWrVtU6dgA6n4KrbMKECTh48CCOHj2KFStWoF27dujTp0+N/fbt2xft2rXDypUrcfToUWRmZtapHJeeno7y8nK89NJLes/Zt28fbt++rf7ZUvHz80NkZKTWz5pMJsOgQYM02jp37qzx/e7UqRMUCoXGfTxx4gQOHDigEf/mzZsRHByMLl26aPz8RUdH63xKMDIyEk2bNtVo27FjByIjI+Ht7a1ua9CgAYYNG6Zx3tatW3H//n2MHTtW41ouLi7o27ev1rVq8zl//PFHuLi4VPs92bx5M2QyGUaPHq1x3RYtWiAkJMSinlql+udY8ylE9mnNmjUICgqCo6MjfHx80LJlS61zGjVqpFWauXLlClq0aAGZTKbR3rx5czg6OuLKlSsa7S1atNDqV9V25coVyOVyreOqUsBrr72G1157TWf8RUVFGu+9vLw03js5Oeltv3Pnjs4+dQkNDQWgTDCffvpptG/fHm+++Sa+++47AMDly5cln+S8YsUKuLi44Mknn1TPL+rcuTPatGmD1NRUzJkzBw4ODrXqS/VL1NfXV+fxPn36IDAwEEuWLMH69euRmJio9b3VRSaTYfz48fj0009x584ddOjQAb17967dB4TyvgHQ+f1XUf0s6frZ9PX1RXp6ukZbo0aN4OLiotHm7Oys9f2eMGECXnrpJZw8eRKPPPIIVq1aBWdnZ/XcLED5M3j69Gk0bNhQZ2xVf/50xXjlyhWtP0AAaLWpft7Dw8N1XqtBA82//WvzOS9fvgxfX1+tr616XSGEzhgBoG3btnq/lmwPEyYiPYKCghAWFlbtObp+cTZr1gz79++HEELj+KVLl3D//n2Nv6YBoLCwUKsPVVuzZs10XlfVx8yZM3XOiwKAhx9+uNrY64O7uzseeeQRnDp1St320EMP6R0pM8SpU6ewe/duANCbiG3duhUDBw6sVX/ff/89AFS7LtL48ePx1ltvQSaTIT4+vtaxjhs3DrNmzcLixYt1zsupjmquXH5+Pvz8/HSeo/r5qDwXTOXChQtaP2u1NXLkSMyYMQOpqal477338PnnnyMuLk5jhMjb2xuurq5aTxFWPl6Zvv9WdM0DqvrfhKov1Rw/KTz00EPYvXs3Kioq9CZN3t7ekMlk2LVrF5ydnbWO62oj28WEiUhi/fv3x/r165GWlobBgwer29esWaM+Xtmvv/6Kixcvqv+KLS8vx7p169CuXTu9owsPP/wwAgMDcfjwYcybN6+ePkndFRUV4ejRo+jZs6e6LSYmBrNmzcK2bdsQGRlp9DVUE7uXLVuG9u3baxy7ffs2nnnmGaxcubJWCVN6ejqWL1+OHj16oFevXnrPi4+Px/79+xEUFIRWrVrVOtZWrVrh73//O06ePFmnRAtQPsHn4OCAlJQUdO/eXec53bt3h6urK7744gsMHTpU3Z6Line truncated
|
||||
"text/plain": [
|
||||
"Figure(PyObject <Figure size 640x480 with 1 Axes>)"
|
||||
]
|
||||
},
|
||||
"metadata": {},
|
||||
"output_type": "display_data"
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"PyObject <matplotlib.legend.Legend object at 0x290b7ddf0>"
|
||||
]
|
||||
},
|
||||
"execution_count": 14,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"semilogy(err1, \"b.-\")\n",
|
||||
"semilogy(err10, \"r.-\")\n",
|
||||
"xlabel(\"iteration\")\n",
|
||||
"ylabel(L\"error $\\Vert x - x_*\\Vert_2$\")\n",
|
||||
"title(\"Problem 3c: ADMM convergence\")\n",
|
||||
"legend([L\"\\rho = 1\", L\"\\rho = 10\"])"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "5132c6e8",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"As desired, the error is clearly converging towards zero — ADMM is converging towards the same solution as 2(c) above.\n",
|
||||
"\n",
|
||||
"Of course, it's not converging nearly as fast as a Newton iteration on the dual problem. In this problem, $\\rho = 10$ is converging significantly faster than $\\rho = 1$. \n",
|
||||
"\n",
|
||||
"**Optional commentary:**\n",
|
||||
"\n",
|
||||
"However, it's not simply a case of \"bigger is better\". For any given problem, there is an optimum $\\rho$ for ADMM convergence. Here, if we increase to $\\rho = 100$ the convergence is better, but if we go to $\\rho = 1000$ it gets *worse*. (You were not required to investigate this, however.) This kind of hyper-parameter tuning is painful to do by hand; a more sophisticated approach might be to try to update $\\rho$ dynamically depending on the convergence of the objective function and the constraint."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 15,
|
||||
"id": "040588ae",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAkwAAAHFCAYAAAAAM6ZOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAA9hAAAPYQGoP6dpAAB/XklEQVR4nO3deViUVfsH8O+wCwguKIqD4IK5oAi4JC4JlkZqqa9mlkqpr5oLrplLtphCpeWCua/tWpk/K5dQSER4FRHMrdQCYUTcUhBcUOb8/qCZGGaGmYFhFvh+rmuu4syZ5znzsMztc59zH4kQQoCIiIiItLIx9wCIiIiILB0DJiIiIiIdGDARERER6cCAiYiIiEgHBkxEREREOjBgIiIiItKBARMRERGRDgyYiIiIiHRgwERERESkAwMmK/e///0Pw4YNQ+PGjeHg4IBGjRph6NChSE5OrtRx16xZg23btqm1Z2ZmQiKRaHzOHN59911IJBK9+ykeDg4OaNasGaZNm4Y7d+4YbTyK67Ns2TKdfbdt2waJRILMzEyjnb+yFixYgMDAQNSrVw9OTk5o3rw5xo8fj8uXL5tsDKdPn4ZEIoG9vT2uXr2qsU/v3r2V30sbGxvUrl0bLVu2xLBhw/Ddd99BLpervcbX1xcSiQS9e/fWeMzPPvtMecxff/1V2a742bGxscFff/2l9rrCwkK4ublBIpHg1VdfrchbJisTExODli1bwsHBARKJxKh/Q8qyxL8TNRUDJisWExOD7t27QyaT4aOPPsLBgwexbNkyXLlyBT169MDq1asrfGxtAVPjxo2RnJyM/v37V2Lk5rN//34kJyfj559/xqBBgxATE4Pw8HBwh6ASd+7cwYgRI7B9+3bs378fs2fPxk8//YSuXbvi1q1bJhnDpk2bAACPHz/GZ599prVf8+bNkZycjKSkJOzevRtz587F/fv3MWzYMPTu3Rt5eXlqr6lduzYSEhLw559/qj23ZcsWuLm5aT2fq6srtm7dqtb+7bff4tGjR7C3t9fn7ZGVS09PR2RkJEJDQxEXF4fk5GTUrl3b3MMiUxBklRITE4WNjY0YMGCAePTokcpzjx49EgMGDBA2NjYiMTGxQsdv166deOqpp4ww0qr1zjvvCH1+jBX9bty4odI+atQoAaDc61RYWKj3eDIyMgQAsXTpUp19t27dKgCIjIwMvY9vDnv37hUAxObNm6v8XA8ePBD169cXAQEBokmTJqJVq1Ya+z311FOiXbt2Gp/bsmWLACBefPFFlXYfHx8RHh4upFKpmD9/vspzly5dEhKJRPz3v/8VAER8fLzyOcXPzrhx44S3t7coLi5WeW2PHj3EiBEjhIuLi4iIiDD8TVcTjx8/Fg8ePDD3MKrcF198IQCIY8eOGe2Y5f2NMfbfCUP+nlVWdfuZ4B0mKxUdHQ2JRIK1a9fCzs5O5Tk7OzusWbMGEokEH3zwgbJdkVpIS0vDkCFD4ObmBnd3d4wcORI3btxQ9vP19cXZs2dx+PBhZYrC19cXgPaUXGJiIvr06YPatWvD2dkZISEh+Pnnn1X6KG4tx8fH4/XXX4eHhwfq16+PIUOGICcnR6Xvjh070LdvXzRu3Bi1atVCmzZtMHfuXBQWFhrh6v3rySefBABlyql3797w9/dHQkICQkJC4OzsjDFjxgAAsrKyMHLkSDRs2BCOjo5o06YNPv74Y43pH7lcjiVLlqBp06ZwcnJCp06dcOjQIb3GdPDgQfTp0wdubm5wdnZG9+7d1V6r+F7+9ttvGDZsGNzd3VGvXj3MnDkTjx8/xh9//IFnn30WtWvXhq+vLz766KMKX6MGDRoAgNrP2ZUrVzB+/Hh4e3vDwcEBXl5eGDp0KK5du1bhc+3evRu3bt3CuHHjEBERgQsXLiAxMdGgY7z22mt47rnn8O2336qlEm1sbDB69Ghs375d5fu2ZcsWeHt74+mnn9Z63DFjxiA7OxuxsbHKNsX4FD8j+pDL5YiJiUHHjh1Rq1Yt1KlTB08++ST27Nmj0uejjz5C69at4ejoiIYNG2L06NGQyWQqx1L8vKakpKBnz55wdnZG8+bN8cEHHyjf340bN+Dg4ICFCxeqjeX333+HRCLBqlWrlG25ubmYMGECpFKpMnX93nvv4fHjx8o+ir8DH330ERYvXoxmzZrB0dER8fHxAID/+7//Q4cOHeDo6IjmzZtj5cqVGtPnQgisWbNGeS3q1q2LoUOHqqU+9XmfCnfu3MGsWbPQvHlz5bV77rnn8Pvvvyv7FBUVYfHixcrr26BBA7z22msqfwc16d27N0aOHAkA6Nq1q1oadsuWLQgICICTkxPq1auHwYMH4/z58yrHePXVV+Hq6orTp0+jb9++qF27Nvr06VPuecuKjY3FCy+8AKlUCicnJ7Rs2RITJkzAzZs3VfoprvnJkycxdOhQ1K1bFy1atAAAPHz4ELNmzUKjRo3g7OyMXr16ITU1Fb6+vmqpZWP8TGiTmZmJ/v37w9XVFfb29vD398d3331n0PUwGXNHbGS4x48fC2dnZ9G1a9dy+3Xp0kU4OzuLx48fCyH+/Zeyj4+PeOONN8SBAwfEJ598IlxcXERgYKAoKioSQghx8uRJ0bx5cxEYGCiSk5NFcnKyOHnypBDi3zsoW7duVZ7n119/Ffb29iI4OFjs2LFD7N69W/Tt21dIJBLxzTffKPsp/qXUvHlzMXXqVHHgwAGxadMmUbduXREaGqoy9vfff18sX75c/Pzzz+LXX38V69atE82aNVPrV9k7TDNmzBAAxC+//CKEKLlzUa9ePeHt7S1iYmJEfHy8OHz4sLh+/bpo0qSJaNCggVi3bp3Yv3+/mDJligAgXn/9deXxFNfH29tb9OjRQ3z//ffi22+/FZ07dxb29vYiKSlJ7XqU/pfj559/LiQSiRg0aJDYtWuX+PHHH8WAAQOEra2tOHjwoNr7eeKJJ8T7778vYmNjxZw5cwQAMWXKFNG6dWuxatUqERsbK1577TUBQHz//fc6r5PCo0ePxL1798TJkydF9+7dRatWrcTdu3eVz8tkMtG4cWPh4eEhPvnkE3Hw4EGxY8cOMWbMGHH+/Hllv6eeekqv74/CM888IxwdHcXff/+tvOvz6quvqvUr7w6TEEKsW7dOABCff/65ss3Hx0f0799fedy9e/cKIUp+n5o0aSLefvtt8e2332q9w3Tjxg3Rs2dPlTtXb775pvD19RVyuVzvO0yjRo0SEolEjBs3Tvzf//2f2Ldvn1iyZIlYuXKlss/48eOV38v9+/eLdevWiQYNGghvb2+Vn+GnnnpK1K9fX/j5+Yl169aJ2NhYMWnSJAFAbN++Xdlv8ODBGu+OzZkzRzg4OIibN28KIYS4evWq8Pb2Fj4+PmL9+vXi4MGD4v333xeOjo4q3wfFz3mTJk1EaGio+O6778Qvv/wiMjIyxL59+4SNjY3o3bu3+OGHH8S3334runbtKnx9fdV+Fv773/8Ke3t7MWvWLLF//37x1VdfidatWwtPT0+Rm5tr8PvMz88X7dq1Ey4uLmLRokXiwIED4vvvvxfTpk0TcXFxQgghiouLxbPPPitcXFzEe++9J2JjY8WmTZtEkyZNRNu2bcW9e/e0fu/Onj0r3nrrLeXfwOTkZHHp0iUhhBBRUVECgBgxYoT4+eefxWeffSaaN28u3N3dxYULF5THiIiIEPb29sLX11dER0eLQ4cOiQMHDmg9p6a/E2vXrhXR0dFiz5494vDhw2L79u0iICBAPPHEE8q/40Ko/s1/8803RWxsrNi9e7cQQogRI0YIGxsbMXfuXPHLL7+IFStWCG9vb+Hu7q7yc2yMnwltMjIyRMOGDUXbtm3F119/Lfbs2SOeffZZIZFIDPp7ZSoMmKxQbm6uACBeeumlcvsNHz5cABDXrl0TQvz7yzNjxgyVfl9++aUAIL744gtlm7aUnKaA6cknnxQNGzZU+UB9/Pix8Pf3F1KpVMjlciHEv7/4kyZNUjnmRx99JACIq1evanwfcrlcPHr0SBw+fFgAEKdOnVI+Z2jAlJubKx49eiRu374tvvjiC1GrVi3h7e0t7t+/L4T49wP+0KFDKq+fO3euxtvwr7/+upBIJOKPP/5QuT5eXl7KYwpR8oe8Xr164umnn1a2lf1DWFhYKOrVqycGDhyoco7i4mIREBAgunTpovZ+Pv74Y5W+HTt2FADErl27lG2PHj0SDRo0EEOGDNF5nYQo+QMJQPno2rWruHLlikqfMWPGCHt7e3Hu3LlyjxUWFiZsbW31Om9mZqawsbFR+bl+6qmnhIuLi8jPz1fpqytg2rdvnwAgPvzwQ2WbImBSvH7o0KFCCCF+/vlnIZFIREZGhs6AaevWrcLR0VHcunVLPH78WDRu3Fi8++67QgihV8CUkJAgAIgFCxZo7XP+/HmNvyfHjh0TAFTSiYqf17I/l23bthX9+vVTfr1nzx6VfxgIUfI76uXlJf7zn/8o2yZMmCBcXV3F5cuXVY63bNkyAUCcPXtWCPHvz3mLFi1UPqCFEKJz587C29tbPHz4UNl29+5dUb9+fZXf1eTkZI0/w9nZ2aJWrVpizpw5Br/PRYsWCQAiNjZWaPP1119r/AdESkqKACDWrFmj9bVC/Pt7m5KSomy7ffu2qFWrlnjuuedU+mZlZQlHR0fx8ssvK9siIiIEALFly5Zyz1P2fNoCD8Xfx8uXLwsA4v/+7/+Uzyl+dt9++22V15w9e1YAEG+++aZKu+LalP45NsbPhDbjxo0TtWrVUgmOHz9+LFq3bi38/Pz0OoYpMSVXjYl/JjKXvQ3+yiuvqHz94osvws7OTuetU00KCwtx7NgxDB06FK6ursp2W1tbjBo1CjKZDH/88YfKa55//nmVrzt06AAAKumTv/76Cy+//DIaNWoEW1tb2Nvb46mnngIAtVvchmjUqBHs7e1Rt25djBw5EkFBQdi/fz+cnJyUferLine truncated
|
||||
"text/plain": [
|
||||
"Figure(PyObject <Figure size 640x480 with 1 Axes>)"
|
||||
]
|
||||
},
|
||||
"metadata": {},
|
||||
"output_type": "display_data"
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"PyObject <matplotlib.legend.Legend object at 0x2a053fd00>"
|
||||
]
|
||||
},
|
||||
"execution_count": 15,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"x = z = s = zeros(n)\n",
|
||||
"err100 = [norm(x - x_opt)]\n",
|
||||
"for i = 2:1000\n",
|
||||
" x,z,s = ADMM(x,z,s; ρ=100)\n",
|
||||
" push!(err100, norm(x - x_opt))\n",
|
||||
"end\n",
|
||||
"\n",
|
||||
"x = z = s = zeros(n)\n",
|
||||
"err1000 = [norm(x - x_opt)]\n",
|
||||
"for i = 2:1000\n",
|
||||
" x,z,s = ADMM(x,z,s; ρ=1000)\n",
|
||||
" push!(err1000, norm(x - x_opt))\n",
|
||||
"end\n",
|
||||
"\n",
|
||||
"semilogy(err10, \"r.-\")\n",
|
||||
"semilogy(err100, \"k.-\")\n",
|
||||
"semilogy(err1000, \"g.-\")\n",
|
||||
"xlabel(\"iteration\")\n",
|
||||
"ylabel(L\"error $\\Vert x - x_*\\Vert_2$\")\n",
|
||||
"title(\"Optional Problem 3c: ADMM convergence for larger ρ\")\n",
|
||||
"legend([L\"\\rho = 10\", L\"\\rho = 100\", L\"\\rho = 1000\"])"
|
||||
]
|
||||
}
|
||||
],
|
||||
"metadata": {
|
||||
"@webio": {
|
||||
"lastCommId": null,
|
||||
"lastKernelId": null
|
||||
},
|
||||
"kernelspec": {
|
||||
"display_name": "Julia 1.8.3",
|
||||
"language": "julia",
|
||||
"name": "julia-1.8"
|
||||
},
|
||||
"language_info": {
|
||||
"file_extension": ".jl",
|
||||
"mimetype": "application/julia",
|
||||
"name": "julia",
|
||||
"version": "1.8.5"
|
||||
}
|
||||
},
|
||||
"nbformat": 4,
|
||||
"nbformat_minor": 5
|
||||
}
|
||||
Reference in new issue
Block a user