PG code discussion
This commit is contained in:
NT committed
1 parent
9a3f1cc46a
commit
fc6547f78f
4 files changed
+327
-168
No files matched your search
+321
-164
@@ -6,15 +6,72 @@
|
||||
"source": [
|
||||
"# Optimizations with Physical Gradients\n",
|
||||
"\n",
|
||||
"Test equations: properly define z in R2 to R2 etc. <br>\n",
|
||||
" <br>\n",
|
||||
" \n",
|
||||
"$\\mathbf{z}(\\mathbf{x}) = \\mathbf{z}(x_0,x_1) = [x_0 \\ x_1^2]^T$ <br>\n",
|
||||
"$L(\\mathbf{z}) = |\\mathbf{z}|^2 = z_0^2 + z_1^2$ <br>\n",
|
||||
"The previous section has made many comments about the advantages and disadvantages of different optimization methods. Below we'll show with a practical example how much differences these properties actually make.\n",
|
||||
"\n",
|
||||
"(Note, to be in sync with python arrays, indices start at 0 here.)\n",
|
||||
"\n",
|
||||
"..."
|
||||
"## Problem Formulation\n",
|
||||
"\n",
|
||||
"We'll consider a very simple setup to clearly illustrate what's happening: we have a two-dimensional input space $\\mathbf{x}$, a mock \"physical model\" likewise with two dimensions $\\mathbf{z}$, and a scalar loss $L$, i.e. \n",
|
||||
"$\\mathbf{x} \\in \\mathbb{R}^2$, \n",
|
||||
"$\\mathbf{z}: \\mathbb{R}^2 \\rightarrow \\mathbb{R}^2 $, and \n",
|
||||
"$L: \\mathbb{R}^2 \\rightarrow \\mathbb{R} $.\n",
|
||||
"The components of a vector like $\\mathbf{x}$ are denoted with $x_i$, and to be in sync with python arrays the indices start at 0.\n",
|
||||
"\n",
|
||||
"Specifically, we'll use the following $\\mathbf{z}$ and $L$:\n",
|
||||
"\n",
|
||||
"$\\quad \\mathbf{z}(\\mathbf{x}) = \\mathbf{z}(x_0,x_1) = \\begin{bmatrix} x_0 \\\\ x_1^2 \\end{bmatrix}$, \n",
|
||||
"i.e. $\\mathbf{z}$ only squares the second component of its input, and\n",
|
||||
"\n",
|
||||
"$\\quad L(\\mathbf{z}) = |\\mathbf{z}|^2 = z_0^2 + z_1^2 \\ $ \n",
|
||||
"represents a simple squared $L^2$ loss.\n",
|
||||
"\n",
|
||||
"As starting point for some example optimizations we'll use \n",
|
||||
"$\\mathbf{x} = \\begin{bmatrix} \n",
|
||||
" 3 \\\\ 3\n",
|
||||
"\\end{bmatrix}$ as initial guess for solving the following simple minimization problem:\n",
|
||||
"\n",
|
||||
"$\\quad \\text{arg min}_{\\mathbf{x}} \\ L(\\mathbf{x}).$\n",
|
||||
"\n",
|
||||
"For us as humans it's quite obvious that $[0 \\ 0]^T$ is the right answer, but let's see how quickly the different optimization algorithms discussed in the previous section can find that solution. And while $\\mathbf{z}$ is a very simple function, it is nonlinear due to its $x_1^2$.\n",
|
||||
"\n",
|
||||
"\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"\n",
|
||||
"## 3 Spaces\n",
|
||||
"\n",
|
||||
"In order to understand the following examples, it's important to keep in mind that we'll deal with mappings between the three _spaces_ we've introduced here:\n",
|
||||
"$\\mathbf{x}$, $\\mathbf{z}$ and $L$. A regular forward pass maps an\n",
|
||||
"$\\mathbf{x}$ to $L$, while for the optimization we'll need to associate values\n",
|
||||
"and changes in $L$ with positions in $\\mathbf{x}$. While doing this, it will \n",
|
||||
"be interesting how this influences the positions in $\\mathbf{z}$ that develop while searching for\n",
|
||||
"the right position in $\\mathbf{x}$.\n",
|
||||
"\n",
|
||||
"```{figure} resources/placeholder.png\n",
|
||||
"---\n",
|
||||
"height: 220px\n",
|
||||
"name: pg-three-spaces\n",
|
||||
"---\n",
|
||||
"TODO, visual overview of 3 spaces\n",
|
||||
"```\n",
|
||||
"\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Implementation\n",
|
||||
"\n",
|
||||
"For this example we'll use the [JAX framework](https://github.com/google/jax), which represents a nice alternative to pytorch and tensorflow for efficiently working with differentiable functions.\n",
|
||||
"\n",
|
||||
"JAX also has a nice numpy wrapper that implements most of numpy's functions. Below we'll use this wrapper as `np`, and the _original_ numpy as `onp`.\n",
|
||||
"\n",
|
||||
"We'll start by defining the $\\mathbf{z}$ and $L$ functions, together with a single composite function `fun` which calls L and z. Having a single native python function is necessary for many of the JAX operations."
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -33,7 +90,9 @@
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"Some test calls:\n"
|
||||
"Starting point x = [3. 3.]\n",
|
||||
"\n",
|
||||
"Some test calls of the functions we defined so far, from top to bottom, z, Lz (manual), Lz:\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -41,8 +100,7 @@
|
||||
"text/plain": [
|
||||
"(DeviceArray([3., 9.], dtype=float32),\n",
|
||||
" DeviceArray(90., dtype=float32),\n",
|
||||
" DeviceArray(90., dtype=float32),\n",
|
||||
" DeviceArray(0, dtype=int32))"
|
||||
" DeviceArray(90., dtype=float32))"
|
||||
]
|
||||
},
|
||||
"execution_count": 1,
|
||||
@@ -71,9 +129,17 @@
|
||||
"\n",
|
||||
"\n",
|
||||
"x = np.asarray([3,3], dtype=np.float32)\n",
|
||||
"print(\"Starting point x = \"+format(x) +\"\\n\")\n",
|
||||
"\n",
|
||||
"print(\"Some test calls:\") \n",
|
||||
"fun_z(x) , fun_L( fun_z(x) ), fun(x), fun([0,0])"
|
||||
"print(\"Some test calls of the functions we defined so far, from top to bottom, z, Lz (manual), Lz:\") \n",
|
||||
"fun_z(x) , fun_L( fun_z(x) ), fun(x) "
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Now we can evaluate the derivatives of our function via `jax.grad`. E.g., `jax.grad(fun_L)(fun_z(x))` evaluates the Jacobian $\\partial L / \\partial z$. The cell below evaluates this and a few variants, together with a sanity check for the inverse of the Jacobian of $\\mathbf{z}$:"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -86,19 +152,21 @@
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"Jacobian L(z): [ 6. 18.]\n",
|
||||
"\n",
|
||||
"Jacobian z(x): \n",
|
||||
"[[1. 0.]\n",
|
||||
" [0. 6.]]\n",
|
||||
"Gradient for full L(x): [ 6. 108.]\n",
|
||||
"\n",
|
||||
"Sanity check with inverse Jacobian of z, this should give x again: [3. 3.]\n"
|
||||
"Sanity check with inverse Jacobian of z, this should give x again: [3. 3.]\n",
|
||||
"\n",
|
||||
"Gradient for full L(x): [ 6. 108.]\n",
|
||||
"\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"\n",
|
||||
"# this works:\n",
|
||||
"print(\"Jacobian L(z): \" + format(jax.grad(fun_L)(fun_z(x))) )\n",
|
||||
"print(\"Jacobian L(z): \" + format(jax.grad(fun_L)(fun_z(x))) +\"\\n\")\n",
|
||||
"\n",
|
||||
"# the following would give an error as z (and hence fun_z) is not scalar\n",
|
||||
"#jax.grad(fun_z)(x) \n",
|
||||
@@ -110,44 +178,43 @@
|
||||
"# the following also gives error, JAX grad needs a single function object\n",
|
||||
"#jax.grad( fun_L(fun_z) )(x) \n",
|
||||
"\n",
|
||||
"print( \"\\nSanity check with inverse Jacobian of z, this should give x again: \" + format(np.linalg.solve(J, np.matmul(J,x) )) +\"\\n\")\n",
|
||||
"\n",
|
||||
"# instead use composite 'fun' from above\n",
|
||||
"print(\"Gradient for full L(x): \" + format( jax.grad(fun)(x) ) )\n",
|
||||
"\n",
|
||||
"print( \"\\nSanity check with inverse Jacobian of z, this should give x again: \" + format(np.linalg.solve(J, np.matmul(J,x) )) )\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": []
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"This already shows a problem, the gradient mostly points along $x_2$ with 108! \n",
|
||||
"We know both dimensions should go towards zero, this points to problems..."
|
||||
"print(\"Gradient for full L(x): \" + format( jax.grad(fun)(x) ) +\"\\n\")\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Gradient descent\n",
|
||||
"The last line is worth a closer look: here we print the gradient $\\partial L / \\partial \\mathbf{x}$ at our initial position. And while we know that we should just move diagonally towards the origin (with the zero vector being the minimizer), this gradient is not very diagonal - it has a strongly dominant component along $x_1$ with an entry of 108.\n",
|
||||
"\n",
|
||||
"Update step:<br>\n",
|
||||
"Let's see how the different methods cope with this situation. We'll compare the first order method _gradient descent_ (i.e., regular, non-stochastic, \"steepest gradient descent\"), _Newton's method_ as a representative of the second order methods, and _physical gradients_.\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Gradient Descent\n",
|
||||
"\n",
|
||||
"$\\Delta \\mathbf{x} = \n",
|
||||
"- \\eta ( J_{L} J_{\\mathbf{z}} )^T\n",
|
||||
"=\n",
|
||||
"- \\eta ( \\frac{\\partial L }{ \\partial \\mathbf{z} } \\frac{\\partial \\mathbf{z} }{ \\partial \\mathbf{x} } )^T$\n",
|
||||
"For gradient descent, the simple gradient based update from equation {eq}`GD-update`\n",
|
||||
"in our setting gives the following update step in $\\mathbf{x}$:\n",
|
||||
"\n",
|
||||
"with step size $\\eta$, \n",
|
||||
"start at $x=[3,3]$\n",
|
||||
"$\\begin{aligned}\n",
|
||||
"\\quad\n",
|
||||
"\\Delta \\mathbf{x} \n",
|
||||
"&= \n",
|
||||
"- \\eta ( J_{L} J_{\\mathbf{z}} )^T \\\\\n",
|
||||
"&=\n",
|
||||
"- \\eta ( \\frac{\\partial L }{ \\partial \\mathbf{z} } \\frac{\\partial \\mathbf{z} }{ \\partial \\mathbf{x} } )^T\n",
|
||||
"\\end{aligned}$\n",
|
||||
"\n",
|
||||
"..."
|
||||
"where $\\eta$ denotes the step size parameter .\n",
|
||||
"\n",
|
||||
"Let's start the optimization via gradient descent at $x=[3,3]$, and update our solution ten times with\n",
|
||||
"$\\eta = 0.01$:\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -175,12 +242,12 @@
|
||||
"source": [
|
||||
"x = np.asarray([3.,3.])\n",
|
||||
"eta = 0.01\n",
|
||||
"history = [x]; updates = []\n",
|
||||
"historyGD = [x]; updatesGD = []\n",
|
||||
"\n",
|
||||
"for i in range(10):\n",
|
||||
" G = jax.grad(fun)(x)\n",
|
||||
" x += -eta * G\n",
|
||||
" history.append(x); updates.append(G)\n",
|
||||
" historyGD.append(x); updatesGD.append(G)\n",
|
||||
" print( \"GD iter %d: \"%i + format(x) )\n"
|
||||
]
|
||||
},
|
||||
@@ -188,8 +255,9 @@
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Let's look at the progression of this...\n",
|
||||
"Target shown in black"
|
||||
"Here we've already printed the resulting positions in $\\mathbf{x}$, and they seem to be going down, i.e. moving in the right direction. The last point, $[2.451 \\ 0.962]$ still has a fair distance of 2.63 to the origin.\n",
|
||||
"\n",
|
||||
"Let's take a look at the progression over the course of the iterations (the evolution was stored in the `history` list above). The blue points denote the positions in $\\mathbf{x}$ from the GD iterations, with the target at the origin shown with a thin black cross."
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -209,7 +277,7 @@
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAFtCAYAAADrr7rKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAcRElEQVR4nO3df5RcZZ3n8fc3TchxId16kEkQgsnIWUThDCBHDesAyrAqkmh2ZfWIcyawnhHHnKjr6poFXWcWJjN6FnEWRlZ0gqxzHP9ZEBRFGSd4ViIqyq81srKCA0iCona3uxIz5Lt/3Gotiqrqru6qvlX9vF/n3NPUc59b95tr+enbz33q3shMJEllWFZ3AZKkxWPoS1JBDH1JKoihL0kFMfQlqSCGviQVxNCXpIIY+pJUkIPqLmCxRUQAzwGm665FkvpoJfDjnOUbt8WFPlXgP1x3EZI0AEcBj3TrUGLoTwM89NBDjI+P112LJC3Y1NQUa9asgTmMYJQY+gCMj48b+pKK44VcSSqIoS9JBTH0Jakghr4kFcTQl6SCGPqSNEQG/TDDWkM/It4WEXdHxFRj2RURr55lm3Mj4vsR8URE3BMRZy9WvZI0CNPTsHUrrFsHa9ZUP7durdr7Lep8Rm5EbACeBH4ABPBHwHuAkzLzf7XpfyrwNWAb8HngTcB/AE7OzHvnuM9xYHJyctJ5+pJqNz0N69fD7t1w4MBv25ctg+OOg127YOXK7u8xNTXFxMQEwERmTnXrW+uZfmbemJk3ZeYPMvN/Z+ZFwC+Bl3bY5B3AlzLzw5m5OzPfD3wH2LJYNUtSP1100dMDH6rXu3fDxRf3d39DM6YfEWMR8UbgEGBXh27rgVta2m5utHd63xURMT6zUN2USJKGwo03Pj3wZxw4ADfc0N/91R76EXFCRPwS2AdcBWzKzO916L4a2NvStrfR3sk2YLJp8WZrkoZCJuzf373P/v39vbhbe+gD9wEnAi8BPgZ8KiJe0Mf33w5MNC1H9fG9JWneImD58u59li+v+vVL7aGfmb/OzPsz847M3AbcRTV2384eYFVL26pGe6f335eZUzML3kdf0hDZsKG6aNvOsmWwcWN/91d76LexDFjRYd0u4MyWtrPofA1AkobapZdWs3Rag39m9s4ll/R3f3XP098eEadFxNrG2P524Azgbxvrr220zfgo8KqIeHdEPD8iPgicAlyx2LVLUj+sXFlNy9yyBdauhSOPrH5u2TK36Zq9qnue/iepztyPoLrIejfwl5n5lcb6ncCDmbm5aZtzgUuAtVTz+9+bmTf1sE/n6UsaWpm9j+H3Mk+/1tCvg6EvaakZmS9nSZIWl6EvSQUx9CWpIIa+JBXE0Jekghj6klQQQ1+SCmLoS1JBDH1JKoihL0kFMfQlqSCGviQVxNCXpIIY+pJUEENfkgpi6EtSQQx9SSqIoS9JBTH0Jakghr4kFcTQl6SCGPqSVBBDX5IKYuhLUkEMfUkqiKEvSQUx9CWpIIa+JBXE0Jekghj6klQQQ1+SCmLoS1JBDH1JKoihL0kFMfQlqSC1hn5EbIuIb0XEdEQ8FhHXR8Sxs2yzOSKyZXlisWqWpFFW95n+6cCVwEuBs4DlwJcj4pBZtpsCjmhanjvIIiVpqTiozp1n5quaX0fEZuAx4EXA17pvmnsGWJokLUl1n+m3mmj8/Nks/Q6NiB9FxEMR8bmIeGGnjhGxIiLGZxZgZd+qlaQRMzShHxHLgMuBr2fmvV263gdcALwWeDPVv+G2iDiqQ/9twGTT8nC/apakUROZWXcNAETEx4BXAy/LzDkHc0QsB3YDn8nM97dZvwJY0dS0Enh4cnKS8fHxBVYtSfWbmppiYmICYCIzp7r1rXVMf0ZEXAGcA5zWS+ADZOb+iPgucEyH9fuAfU37WkipkjTS6p6yGY3A3wS8IjMfmMd7jAEnAI/2uz5JWmrqPtO/EngT1fj8dESsbrRPZuavACLiWuCRzNzWeP0B4BvA/cAzgfdQTdn8xOKWLkmjp+7Qf1vj586W9vOBaxr/fTRwoGnds4CrgdXAz4E7gFMz83sDq1KSloihuZC7WBrTNie9kCtpqejlQu7QTNmUJA2eoS9JBTH0Jakghr4kFcTQl6SCGPqSVBBDX5IKYuhLUkEMfUkqiKEvSQUx9CWpIIa+JBXE0Jekghj6klQQQ1+SCmLoS1JBDH1JKoihL0kFMfQlqSCGviQVxNCXpIIY+pJUEENfkgpi6EtSQQx9SSqIoS9JBTH0Jakghr4kFcTQl6SCGPqSVBBDX5IKYuhLUkEMfUkqiKEvSQUx9CWpILWGfkRsi4hvRcR0RDwWEddHxLFz2O7ciPh+RDwREfdExNmLUa8kjbq6z/RPB64EXgqcBSwHvhwRh3TaICJOBT4DfBI4CbgeuD4ijh94tZI04iIz667hNyLicOAx4PTM/FqHPp8FDsnMc5ravgHcmZkXzmEf48Dk5OQk4+PjfapckuozNTXFxMQEwERmTnXrW/eZfquJxs+fdemzHrilpe3mRvvTRMSKiBifWYCVCy9TkkbT0IR+RCwDLge+npn3dum6Gtjb0ra30d7ONmCyaXl4YZVK0ugamtCnGts/Hnhjn993O9VfEDPLUX1+f0kaGQfVXQBARFwBnAOclpmznYnvAVa1tK1qtD9NZu4D9jXtawGVSlL/ZMJiR1LdUzajEfibgFdk5gNz2GwXcGZL21mNdkkaatPTsHUrrFsHa9ZUP7durdoXQ62zdyLir4E3Aa8F7mtaNZmZv2r0uRZ4JDO3NV6fCtwKvA/4AtVw0H8ETp7lWsDMPp29I6kW09Owfj3s3g0HDvy2fdkyOO442LULVs5jqskozd55G9U4+07g0ablDU19jgaOmHmRmbdR/aL4Y+Au4PXA6+YS+JJUp4suenrgQ/V69264+OLB1zBU8/QXg2f6kuqybh08+GDn9WvXwgNzGeRuMUpn+pJUhEzYv797n/37q36DZOhL0iKIgOXLu/dZvnzws3kMfUlaJBs2VBdt21m2DDZuHHwNhr4kLZJLL61m6bQG/8zsnUsuGXwNhr4kLZKVK6tpmVu2VBdtjzyy+rlly/yna/bK2TuSVJN+fSPX2TuSNALquCuMoS9JBTH0Jakghr4kFcTQl6SCGPqSVBBDX5IKYuhLUkEMfUkqiKEvSQUx9CWpIIa+JBXE0Jekghj6klQQQ1+SCmLoS9KADdNjSwx9SRqA6WnYuhXWrYM1a6qfW7dW7XXyyVmS1GfT07B+PezeDQcO/LZ95lm4/X40ok/OkqQaXXTR0wMfqte7d8PFF9dTFxj6ktR3N9749MCfceAA3HDD4tbTzNCXpD7KhP37u/fZv7++i7uGviT1UQQsX969z/Ll9TwUHQx9Seq7DRuqi7btLFsGGzcubj1P2X99u5akpenSS6tZOq3BPzN755JL6qkLDH1J6ruVK6tpmVu2wNq1cOSR1c8tW/o/XbNXztOXpAHLHOwYvvP0JWmI1HXRth1DX5IK0rfQj4jjIuKHPW5zWkTcGBE/joiMiNfN0v+MRr/WZfWCipekQvTzTP9g4Lk9bnMIcBfw9h63OxY4oml5rMftJalIB821Y0RcNkuXw3vdeWZ+Efhi4/172fSxzPxFr/uTpNLNOfSBdwB3Ap2uDB+64Grm7s6IWAHcC3wwM7++iPuWpI4GPVNnoXoJ/fuBj2Tmp9utjIgTgTv6UVQXjwIXAt8GVgBvAXZGxEsy8zsd6lrR6Dujxhmykpai6enqzpo33ljdV2f58upbuZdeWu+c/HZ6Cf1vAy8C2oY+kMBAf79l5n3AfU1Nt0XE84B3AX/YYbNtwH8aZF2SytXp3vlXXglf/Wr9X8Zq1cuF3HcDl3damZl3ZWYdU0C/CRzTZf12YKJpOWoxipJUhmG+d347cw7pzNyTmT+KiJd36hMRb+1PWT05kWrYp63M3JeZUzMLUPPDyiQtJcN87/x25nNm/qWI+HBE/ObmoRHx7Ii4EfiLXt4oIg6NiBMb1wMA1jVeH91Yvz0irm3q/86IeG1EHBMRx0fE5cArgCvn8e+QpAUZ9nvntzOf0H85sAn4VkS8ICJeQzWLZpzqrLsXpwDfbSwAlzX++88ar48Ajm7qfzDwX4B7gFuB3wP+IDP/vvd/hiQtzLDfO7+dXi7kApCZtzXOzK8CvkP1i+P9wIeyx7u3ZeZOulz8zczNLa8/BHyot4olaXA2bKgu2rYb4qn73vntzPfC6z+nOkt/GPgnqm/I/rN+FSVJo2KY753fTs+hHxHvA3YBXwGOB14MnATcHRHr+1ueJA23Yb53fjs9308/Ih4FLmjcQmGmbTnw58DWzFzRceMh4P30JQ1SHd/I7eV++j2P6QMnZOZPmxsycz/wnoj4/DzeT5KWjGG6aNtOz8M7rYHfsu7WhZUjSaNhmKZh9sKHqEjSHE1Pw9atsG4drFlT/dy6tWofFT4jV5LmoNM9dmZm6dR50dZn5EpSn43aPXY6MfQlaQ5G7R47nRj6kjSLUbzHTieGviTNYhTvsdOJoS9Jc7Bhw9NvtTBjGO+x04mhL0mzyBy9e+x0Mp9v5ErSktfuubevfCX8/u/Dl77027aNG6vAH7Z77HRi6EtSi05z8q++ujqrv/tuOPTQ0RjDb+XwjiS1mMuc/FEMfDD0Jelplsqc/HYMfUlqspTm5Ldj6EtSk6U0J78dQ1+SmmQunTn57Th7R1LxWqdnjo3BM58JP//5U4dxRm1OfjuGvqSidbtl8rOeVU3NfPLJ0ZyT346hL6lo3aZn/uIX8OY3w+WXj+4YfivH9CUVbS7TM5dK4IOhL6lgS316ZjuGvqSiLeXpme0Y+pKK0vpw88cf79x31KdntuOFXEnF6DRTp52lMD2zHUNfUjE6zdSBagjn0ENhfHzpTM9sx9CXVIxuM3Uy4bDD4Ic/XFpj+K0c05dUhAMH5jZTZ6kz9CUtWc0XbY8+GvbLine truncated
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAFtCAYAAADrr7rKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAcCElEQVR4nO3df5BdZZ3n8fc3IQQ3STfWqAlCQjJSqyjU8KvUsA6ibFZFgrIrq+U4RdadXXBMRV1X1yzo6CxMZrQWcZaMUOgo7Ew5/rOFxF8I46C1ElGj/FojKys4oCSMP+hud03MkO/+cW7r5ebe23277+1zbz/vV9Wp5jz3Ofd8c7x++vRznntOZCaSpDIsqbsASdLCMfQlqSCGviQVxNCXpIIY+pJUEENfkgpi6EtSQQx9SSrIUXUXsNAiIoBnA1N11yJJfbQK+HHO8I3b4kKfKvAfrbsISRqAE4AfdetQYuhPATzyyCOMjY3VXYskzdvk5CRr166FWYxglBj6AIyNjRn6korjhVxJKoihL0kFMfQlqSCGviQVxNCXpIIY+pI0RAb9MMNaQz8i3hIR90bEZGPZHRGvmmGbiyPiexFxICLui4jzF6peSRqEqSnYtg02bIC1a6uf27ZV7f0WdT4jNyI2A08C3wcCuAR4F3B6Zv6vNv3PBr4KbAc+C7wR+E/AGZl5/yz3OQZMTExMOE9fUu2mpmDjRti7Fw4f/k37kiVw8smwezesWtX9PSYnJxkfHwcYz8zJbn1rPdPPzF2Z+fnM/H5m/u/MvBz4BfDiDpu8DfhiZn4oM/dm5nuBbwNbF6pmSeqnyy8/MvChWt+7F664or/7G5ox/YhYGhFvAFYAuzt02wjc3tJ2a6O90/suj4ix6YXqpkSSNBR27Toy8KcdPgy33NLf/dUe+hFxakT8AjgIXAdclJnf7dB9DbC/pW1/o72T7cBE0+LN1iQNhUw4dKh7n0OH+ntxt/bQBx4ATgNeBHwUuDEint/H998BjDctJ/TxvSVpziJg2bLufZYtq/r1S+2hn5m/yswHM3NPZm4H7qEau29nH7C6pW11o73T+x/MzMnpBe+jL2mIbN5cXbRtZ8kSuPDC/u6v9tBvYwmwvMNru4HzWto20fkagCQNtauuqmbptAb/9OydK6/s7/7qnqe/IyLOiYj1jbH9HcC5wF83Xr+p0TbtI8ArI+KdEfG8iHg/cBZw7ULXLkn9sGpVNS1z61ZYvx6OP776uXXr7KZr9qruefofpzpzP47qIuu9wJ9l5m2N1+8AHs7MLU3bXAxcCaynmt//7sz8fA/7dJ6+pKGV2fsYfi/z9GsN/ToY+pIWm5H5cpYkaWEZ+pJUEENfkgpi6EtSQQx9SSqIoS9JBTH0Jakghr4kFcTQl6SCGPqSVBBDX5IKYuhLUkEMfUkqiKEvSQUx9CWpIIa+JBXE0Jekghj6klQQQ1+SCmLoS1JBDH1JKoihL0kFMfQlqSCGviQVxNCXpIIY+pJUEENfkgpi6EtSQQx9SSqIoS9JBTH0Jakghr4kFcTQl6SCGPqSVBBDX5IKUmvoR8T2iPhmRExFxOMRcXNEPHeGbbZERLYsBxaqZkkaZXWf6b8U2Am8GNgELAO+FBErZthuEjiuaTlxkEVK0mJxVJ07z8xXNq9HxBbgceBM4KvdN819AyxNkhalus/0W403fv5shn4rI+KHEfFIRHwmIl7QqWNELI+IsekFWNW3aiVpxAxN6EfEEuAa4GuZeX+Xrg8AbwZeA7yJ6t9wZ0Sc0KH/dmCiaXm0XzVL0qiJzKy7BgAi4qPAq4CXZOasgzkilgF7gU9l5nvbvL4cWN7UtAp4dGJigrGxsXlWLUn1m5ycZHx8HGA8Mye79a11TH9aRFwLXACc00vgA2TmoYj4DnBSh9cPAgeb9jWfUiVppNU9ZTMagX8R8PLMfGgO77EUOBV4rN/1SdJiU/eZ/k7gjVTj81MRsabRPpGZvwSIiJuAH2Xm9sb6+4CvAw8CxwLvopqy+bGFLV2SRk/dof+Wxs87Wtr/DfDJxn+vAw43vfZ04AZgDfBzYA9wdmZ+d2BVStIiMTQXchdKY9rmhBdyJS0WvVzIHZopm5KkwTP0Jakghr4kFcTQl6SCGPqSVBBDX5IKYuhLUkEMfUkqiKEvSQUx9CWpIIa+JBXE0Jekghj6klQQQ1+SCmLoS1JBDH1JKoihL0kFMfQlqSCGviQVxNCXpIIY+pJUEENfkgpi6EtSQQx9SSqIoS9JBTH0Jakghr4kFcTQl6SCGPqSVBBDX5IKYuhLUkEMfUkqiKEvSQUx9CWpIIa+JBWk1tCPiO0R8c2ImIqIxyPi5oh47iy2uzgivhcRByLivog4fyHqlaRRV/eZ/kuBncCLgU3AMuBLEbGi0wYRcTbwKeDjwOnAzcDNEXHKwKuVpBEXmVl3Db8WEc8EHgdemplf7dDn08CKzLygqe3rwN2Zedks9jEGTExMTDA2NtanyiWpPpOTk4yPjwOMZ+Zkt751n+m3Gm/8/FmXPhuB21vabm20HyEilkfE2PQCrJp/mZI0moYm9CNiCXAN8LXMvL9L1zXA/pa2/Y32drYDE03Lo/OrVJJG19CEPtXY/inAG/r8vjuo/oKYXk7o8/tL0sg4qu4CACLiWuAC4JzMnOlMfB+wuqVtdaP9CJl5EDjYtK95VCpJ/ZMJCx1JdU/ZjEbgXwS8PDMfmsVmu4HzWto2NdolaahNTcG2bbBhA6xdW/3ctq1qXwi1zt6JiL8A3gi8Bnig6aWJzPxlo89NwI8yc3tj/WzgK8B7gM9RDQf9Z+CMGa4FTO/T2TuSajE1BRs3wt69cPjwb9qXLIGTT4bdu2HVHKaajNLsnbdQjbPfATzWtLy+qc864Ljplcy8k+oXxb8H7gFeB7x2NoEvSXW6/PIjAx+q9b174YorBl/DUM3TXwie6Uuqy4YN8PDDnV9fvx4ems0gd4tROtOXpCJkwqFD3fscOlT1GyRDX5IWQAQsW9a9z7Jlg5/NY+hL0gLZvLm6aNvOkiVw4YWDr8HQl6QFctVV1Syd1uCfnr1z5ZWDr8HQl6QFsmpVNS1z69bqou3xx1c/t26d+3TNXjl7R5Jq0q9v5Dp7R5JGQB13hTH0Jakghr4kFcTQl6SCGPqSVBBDX5IKYuhLUkEMfUkqiKEvSQUx9CWpIIa+JBXE0Jekghj6klQQQ1+SCmLoS1JBDH1JGrBhemyJoS9JAzA1Bdu2wYYNsHZt9XPbtqq9Tj45S5L6bGoKNm6EvXvh8OHftE8/C7ffj0b0yVmSVKPLLz8y8KFa37sXrriinrrA0Jekvtu168jAn3b4MNxyy8LW08zQl6Q+yoRDh7r3OXSovou7hr4k9VEELFvWvc+yZfU8FB0MfUnqu82bq4u27SxZAhdeuLD1PGX/9e1akhanq66qZum0Bv/07J0rr6ynLjD0JanvVq2qpmVu3Qrr18Pxx1c/t27t/3TNXjlPX5IGLHOwY/jO05ekIVLXRdt2DH1JKkjfQj8iTo6IH/S4zTkRsSsifhwRGRGvnaH/uY1+rcuaeRUvSYXo55n+0cCJPW6zArgHeGuP2z0XOK5pebzH7SWpSEfNtmNEXD1Dl2f2uvPM/ALwhcb797Lp45n5RK/7k6TSzTr0gbcBdwOdrgyvnHc1s3d3RCwH7gfen5lfW8B9S1JHg56pM1+9hP6DwIcz86/avRgRpwF7+lFUF48BlwHfApYDfwDcEREvysxvd6hreaPvtBpnyEpajKamqjtr7tpV3Vdn2bLqW7lXXVXvnPx2egn9bwFnAm1DH0hgoL/fMvMB4IGmpjsj4jnAO4Df77DZduCPBlmXpHJ1unf+zp3w5S/X/2WsVr1cyH0ncE2nFzPznsysYwroN4CTury+AxhvWk5YiKIklWGY753fzqxDOjP3ZeYPI+JlnfpExKX9Kasnp1EN+7SVmQczc3J6AWp+WJmkxWSY753fzlzOzL8YER+KiF/fPDQinhERu4A/7eWNImJlRJzWuB4AsKGxvq7x+o6IuKmp/9sj4jURcVJEnBIR1wAvB3bO4d8hSfMy7PfOb2cuof8y4CLgmxHx/Ih4NdUsmjGqs+5enAV8p7EAXN347z9urB8HrGvqfzTwX4H7gK8AvwP888z8297/GZI0P8N+7/x2ermQC0Bm3tk4M78O+DbVL473Ah/MHu/elpl30OXib2ZuaVn/IPDB3iqWpMHZvLm6aNtuiKfue+e3M9cLr/+U6iz9UeAfqb4h+0/6VZQkjYphvnd+Oz2HfkS8B9gN3AacArwQOB24NyI29rc8SRpuw3zv/HZ6vp9+RDwGvLlxC4XptmXAnwDbMnN5x42HgPfTlzRIdXwjt5f76fc8pg+cmpk/aW7IzEPAuyLis3N4P0laNIbpom07PQ/vtAZ+y2tfmV85kjQahmkaZi98iIokzdLUFGzbBhs2wNq11c9t26r2UeEzciVpFjrdY2d6lk6dF219Rq4k9dmo3WOnE0NfkmZh1O6x04mhL0kzGMV77HRi6EvSDEbxHjudGPqSNAubNx95q4Vpw3iPnU4MfUmaQebo3WOnk7l8I1eSFr12z719xSvgd38XvvjF37RdeGEV+MN2j51ODH1JatFpTv4NN1Rn9ffeCytXjsYYfiuHdySpxWzm5I9i4IOhL0lHWCxz8tsx9CWpyWKak9+OoS9JTRbTnPx2DH1JapK5eObkt+PsHUnFa52euXQpHHss/PznTx3GGbU5+e0Y+pKK1u2WyU9/ejU188knR3NOfjuGvqSidZue+cQT8KY3wTXXjO4YfivH9CUVbTbTMxdL4IOhL6lgi316ZjuGvqSiLebpme0Y+pKK0vpw85/+tHPfUZ+e2Y4XciUVo9NMnXYWw/TMdgx9ScXoNFMHqiGclSthbGzxTM9sx9CXVIxuM3Uy4bd+C37wg8U1ht/KMX1JRTh8eHYzdRY7Q1/SotV80XbdOtiLine truncated
|
||||
"text/plain": [
|
||||
"<Figure size 400x400 with 1 Axes>"
|
||||
]
|
||||
@@ -223,34 +291,46 @@
|
||||
"source": [
|
||||
"import matplotlib.pyplot as plt\n",
|
||||
"axes = plt.figure(figsize=(4, 4), dpi=100).gca()\n",
|
||||
"historyGD = onp.asarray(history)\n",
|
||||
"updatesGD = onp.asarray(updates) # for later\n",
|
||||
"historyGD = onp.asarray(historyGD)\n",
|
||||
"updatesGD = onp.asarray(updatesGD) # for later\n",
|
||||
"axes.scatter(historyGD[:,0], historyGD[:,1], lw=0.5, color='blue')\n",
|
||||
"axes.scatter([0], [0], lw=0.25, color='black') # target at 0,0\n",
|
||||
"axes.scatter([0], [0], lw=0.25, color='black', marker='x') # target at 0,0\n",
|
||||
"axes.set_xlabel('x0')\n",
|
||||
"axes.set_ylabel('x1')"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"No surprise here: the initial step mostly moves downwards along $x_1$ (in top right corner), and the updates afterwards curve towards the origin. But they don't get very far. It's still quite a distance to the solution in the bottom left corner."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Newton\n",
|
||||
"\n",
|
||||
"For Newton's method, let's consider the combined function, i.e., \n",
|
||||
"$L(\\mathbf{x}) = x_0^2 + x_1^4$\n",
|
||||
"...\n",
|
||||
"For Newton's method, the update step is given by\n",
|
||||
"\n",
|
||||
"Newton's method update step is given by\n",
|
||||
"$\n",
|
||||
"\\Delta \\mathbf{x} = \n",
|
||||
"$\\begin{aligned}\n",
|
||||
"\\quad\n",
|
||||
"\\Delta \\mathbf{x} &= \n",
|
||||
"- \\eta \\left( \\frac{\\partial^2 L }{ \\partial \\mathbf{x}^2 } \\right)^{-1}\n",
|
||||
" \\frac{\\partial L }{ \\partial \\mathbf{x} }\n",
|
||||
"=\n",
|
||||
"\\\\\n",
|
||||
"&=\n",
|
||||
"- \\eta \\ H_L^{-1} \\ ( J_{L} J_{\\mathbf{z}} )^T\n",
|
||||
"$\n",
|
||||
"\\end{aligned}$\n",
|
||||
"\n",
|
||||
"Hence we need to compute the Hessian and its inverse, easy in JAX, and very cheap here, but not so simple in practice...\n"
|
||||
"Hence, in addition to the same gradient as for GD, we now need to evaluate and invert the Hessian of $\\frac{\\partial^2 L }{ \\partial \\mathbf{x}^2 }$.\n",
|
||||
"\n",
|
||||
"This is quite straightforward in JAX: we can call `jax.jacobian` two times, and then use the JAX version of `linalg.inv` to invert the matrix.\n",
|
||||
"\n",
|
||||
"For the optimization with Newton's method we'll use a larger step size of $\\eta =1/3$. For this example and the following one, we've chosen the stepsize such that the magnitude of the first update step is roughly the same as the one of GD. In this way, we can compare the trajectories of all three methods relative to each other. Note that this is by no means meant to illustrate or compare the stability of the methods in this example. Stability and upper limits for $\\eta$ are separate topics. Here we're focusing on convergence properties.\n",
|
||||
"\n",
|
||||
"In the next cell, we apply the Newton updates ten times starting from the same initial guess:"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -278,16 +358,16 @@
|
||||
"source": [
|
||||
"x = np.asarray([3.,3.])\n",
|
||||
"eta = 1./3.\n",
|
||||
"history = [x]; updates = []\n",
|
||||
"historyNt = [x]; updatesNt = []\n",
|
||||
"\n",
|
||||
"for i in range(10):\n",
|
||||
" G = jax.grad(fun)(x)\n",
|
||||
" H = jax.jacobian(jax.jacobian(fun))(x)\n",
|
||||
" #H = jax.jacfwd(jax.jacrev(fun_Lz))(x) # more efficient? \n",
|
||||
" #H = jax.jacfwd(jax.jacrev(fun_Lz))(x) # alternative\n",
|
||||
" Hinv = np.linalg.inv(H)\n",
|
||||
" \n",
|
||||
" x += -eta * np.matmul( Hinv , G)\n",
|
||||
" history.append(x); updates.append( np.matmul( Hinv , G) )\n",
|
||||
" historyNt.append(x); updatesNt.append( np.matmul( Hinv , G) )\n",
|
||||
" print( \"Newton iter %d: \"%i + format(x) )\n",
|
||||
"\n"
|
||||
]
|
||||
@@ -296,7 +376,9 @@
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Let's compare, Newton actually does better, note - uses different step size..."
|
||||
"The last line already indicates: Newton's method does quite a bit better. The last point $[0.052 \\ 0.924]$ only has a distance of 0.925 to the origin (compared to 2.63 for GD)\n",
|
||||
"\n",
|
||||
"Below, we plot the Newton trajectory in orange next to the GD version in blue."
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -316,7 +398,7 @@
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAFtCAYAAADrr7rKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAg5ElEQVR4nO3df7TcdX3n8ef7XkKWDfdeLdKAAUwKBxcLDQqnNawHUJa2oglkt6wesW1gPRbXnKjr2jUL7bZdUvrjVKVrWqu2Qeoe6z8LJoqFWhp6KhEFJZU1QrGgEElQwHuvXYkX7nv/+M7V4Wbm3js3M/P9zv0+H+d8zzCf+Xxn3hkmr3zn+/3M5xOZiSSpHobKLkCS1D+GviTViKEvSTVi6EtSjRj6klQjhr4k1YihL0k1YuhLUo0cVXYB/RYRAbwYmCy7FknqohHg2znPL25rF/oUgf9Y2UVIUg+cBOyfq0MdQ38S4NFHH2V0dLTsWiTpiE1MTHDyySfDAs5g1DH0ARgdHTX0JdWOF3IlqUYMfUmqEUNfkmrE0JekGjH0JalGDH1JqpCc7u1qhqWGfkS8LSL+MSImGtueiHjtPPtcHhFfj4hnIuKrEXFJv+qVpF6YfHqS3X+0hUe3r+HAn53Mo9vXsPuPtjD5dPcnDogy18iNiPXAc8A/AQH8KvAe4OWZ+X9b9D8P+HtgK/Bp4E3AfwNekZn3L/A1R4Hx8fFxx+lLKt3k05McuGkdP/UT+xgenv5R+7PPDfHwU2dwwq/sYeSFI3M+x8TEBGNjYwBjmTkxV99Sj/Qzc1dm3pqZ/5SZD2bmNcD3gVe22eUdwF9n5h9m5r7M/A3gy8DmftUsSd10719cc1jgAxw1PM2an9jHvTuu7errVeacfkQMR8QbgRXAnjbd1gGfm9V2W6O93fMuj4jRmY1iUiJJqoRT/9WuwwJ/xlHD05x69M6uvl7poR8RZ0XE94FDwIeAjZn5tTbdTwAOzmo72GhvZysw3rQ52ZqkSsjp5KihqTn7HDU01dWLu6WHPvAAcDbwc8CfAh+LiJd18fmvB8aatpO6+NyStGgxFDw7vWzOPs/mMmIouvaapYd+Zv4wMx/KzHszcyuwl+LcfSsHgJWz2lY22ts9/6HMnJjZcB59SRXyjWfW8+xzraP42eeG+MahDV19vdJDv4UhYHmbx/YAF81qu5j21wAkqdLOuWobDz91xmHB/+xzQzz85Bmcc+V1XX29UqdWjojrgc8C36K4wPom4ELgFxqP3wTsb3wDALgBuDMi3g18BngjcC7w1v5WLkndMfLCEfiVPfzDjms59eidHDU0xbPTy/jGDzdwzpXXzTtcs1Nlj9P/c4oj9xMpLrL+I/D7mfk3jcd3A49k5qamfS4HrgNWU4zv//XMvLWD13ScvqTKyuns+Bx+J+P0Sw39Mhj6kpaagflxliSpvwx9SaoRQ1+SasTQl6QaMfQlqUYMfUmqEUNfkmrE0JekGjH0JalGDH1JqhFDX5JqxNCXpBox9CWpRgx9SaoRQ1+SasTQl6QaMfQlqUYMfUmqEUNfkmrE0JfUWs3Wz66Lo8ouQFKFTE3C3mtg/y6YnoKhZbBqPazdBstGyq5OXWDoSypMTcLt62B8HzD94/YHt8PBO+Dn9xj8S4CndyQV9l5zeOBDcX98H+y9toyq1GWGvqTC/l0cHvgzpmH/zn5Wox4x9CUVF22np+buMz3lxd0lwNCXBBHFRdu5DC0r+mmgGfqSCqvW0z4ShmDVhn5Wox4x9CUV1m6DsTM4PBaGiva115VRlbrM0JdUWDZSDMs8fTOsWA3HrCpuT9/scM0lJLJmF2YiYhQYHx8fZ3R0tOxypOrK9Bz+gJiYmGBsbAxgLDMn5urrkb6k1gz8JcnQl6QaMfQlqUYMfUmqEUNfkmqk1NCPiK0R8aWImIyIJyLiloh46Tz7bIqInLU906+aJWmQlX2kfwGwHXglcDGwDLg9IlbMs98EcGLT9pJeFilJS0Wp8+ln5i8234+ITcATwDnA38+9ax7oYWmStCSVfaQ/21jj9ql5+h0bEd+MiEcj4lMR8dPtOkbE8ogYndkAf1YoqbYqE/oRMQR8APh8Zt4/R9cHgKuAS4E3U/wZ7oqIk9r03wqMN22PdatmSRo0lZmGISL+FHgt8KrMXHAwR8QyYB/wicz8jRaPLweWNzWNAI85DYOkpaKTaRgqsUZuRHwQeD1wfieBD5CZUxHxFeC0No8fAg41vdaRlCpJA63sIZvRCPyNwGsy8+FFPMcwcBbweLfrk6Slpuwj/e3AmyjOz09GxAmN9vHM/AFARNwE7M/MrY37vwl8AXgIeAHwHoohmx/tb+mSNHjKDv23NW53z2q/Erix8d+n8PzVml8IfAQ4AXgauBc4LzO/1rMqJWmJqMyF3H5xPn1JS43z6UvdVLMDIy1tZZ/ekappahL2XgP7d8H0FAwtKxYOX7vNZQM10Ax9abapSbh9HYzv43mXkx7cDgfvcL1YDTRP70iz7b3m8MCH4v74Pth7bRlVSV1h6Euz7d/F4YE/Yxr27+xnNVJXGfpSs8ziHP5cpqe8uKuBZehLzSKKi7ZzGVpW9JMGkKEvzbZqPe3/agzBqg39rEbqKkNfmm3tNhg7g8P/egwV7WuvK6MqqSsMfWm2ZSPFsMzTN8OK1XDMquL29M0O19TAcxoGaT6ZnsNXpTkNg9RNBr6WEENfkmrE0JekGjH0JalGDH1JqhFDX5JqxNCXpBox9CWpRgx9SaoRQ1+SasTQl6QaMfQlqUYMfUmqEUNf5anZDK9SFRxVdgGqmalJ2HtNsfj49FSx9OCq9cXCJc5TL/Wcoa/+mZqE29fB+D5g+sftD26Hg3e4QInUB57eUf/svebwwIfi/vg+2HttGVVJtWLoq3/27+LwwJ8xDft39rMaqZYMffVHZnEOfy7TU17clXrM0Fd/RBQXbecytMylCaUeM/TVP6vW0/4jNwSrNvSzGqmWDH31z9ptMHYGh3/shor2tdeVUZVUK4a++mfZSDEs8/TNsGI1HLOquD19s8M1pT6JrNmFs4gYBcbHx8cZHR0tu5x6y/QcvtQFExMTjI2NAYxl5sRcfUs90o+IrRHxpYiYjIgnIuKWiHjpAva7PCK+HhHPRMRXI+KSftSrLjPwpb4r+/TOBcB24JXAxcAy4PaIWNFuh4g4D/gE8OfAy4FbgFsi4syeVytJA65Sp3ci4njgCeCCzPz7Nn0+CazIzNc3tX0BuC8zr17Aa3h6R9KSMjCnd1oYa9w+NUefdcDnZrXd1mg/TEQsj4jRmQ3waqGk2qpM6EfEEPAB4POZef8cXU8ADs5qO9hob2UrMN60PXZklUrS4KpM6FOc2z8TeGOXn/d6im8QM9tJXX5+SRoYlZhaOSI+CLweOD8z5zsSPwCsnNW2stF+mMw8BBxqeq0jqFSSuqeMUctlD9mMRuBvBF6TmQ8vYLc9wEWz2i5utEtSpU1OwpYtsGYNnHxycbtlS9HeD6WO3omIPwHeBFwKPND00Hhm/qDR5yZgf2Zubdw/D7gTeC/wGYrTQf8deMU81wJmXtPRO5JKMTkJ69bBvn0w3TTL+NAQnHEG7NkDI4sYajJIo3feRnGefTfweNP2hqY+pwAnztzJzLso/qF4K7AX+CXgsoUEviSV6ZprDg98KO7v2wfX9mEdoUqN0+8Hj/QllWXNGnjkkfaPr14NDy/kJPcsg3Skryqp2QGA1E+ZMDXPOkJTfVhHqBKjd1Siqcli7dr9u4qVq4aWFfPer93mrJdSF0XAsnnWEVrWh3WEPNKvs6lJuH0dPLgd/uUR+MH+4vbB7UX7VJ+GE0g1sX59cdG2laEh2NCHdYQM/Trbew2M7+Pwxcqni/a9fbiqJNXItm3FKJ3ZwT8zeue6PqwjZOjX2f5dHB74M6Zh/85+ViMteSMjxbDMzZuLi7arVhW3mzcvfrhmpzynX1eZxTn8uUxPudCJ1GUjI3DDDcVWu1/kqkQRxUXbuQz14aqSVGNl/PUy9Ots1XrafwSGYFUfripJ6itDv87WboOxMzj8YzBUtK/tw1UlSX1l6NfZshH4+T1w+mZYsRqOWVXcnr65aHecvrTkOA2DfsyLttJAchoGLY6BLy15hr4k1YihL0k1YuhLUo0Y+pJUI4a+JNWIoS9JNWLoS1KNGPqDrmY/rpN0ZJxaeRC5xKGkRTL0B83MEoezV7x6cDscvMM5c6QKqtIMJ57eGTQucSgNhMlJ2LIF1qyBk08ubrdsKdrL5IRrg+ZTa4rFy9tZsRoufbhf1UhqYXIS1q2Dfftguun4bGYt3G4vjeiEa0tVJ0scSirNNdccHvhQ3N+3D64t8Qu5oT9IXOJQGgi7dh0e+DOmp2Hnzv7W08zQHzQucShVWiZMzfOFfKrEL+SG/qBxiUOp0iJg2TxfyJeV+IXc0B80LnEoVd769cVF21aGhmBDiV/IHb0z6Ko0AFgS4Ogd9ZKBL1XOyEgR7Js3w+rVsGpVcbt5c/cDv1Me6UtSj/X6C7lH+pJUIVX6Qm7oS1KNdC30I+KMiPjnDvc5PyJ2RcS3IyIj4rJ5+l/Y6Dd7O+GIipekmujmkf7RwEs63GcFsBd4e4f7vRQ4sWl7osP9JamWFjy1ckS8b54ux3f64pn5WeCzjefvZNcnMvN7nb7eQHAIpqQe6mQ+/XcA9wHtrgwfe8TVLNx9EbEcuB/4rcz8fB9fu/tcFEVaMqp+3NZJ6D8EvD8zP97qwYg4G7i3G0XN4XHgauAeYDnwFmB3RPxcZn65TV3LG31nVCtFXRRFGniTk8XMmrt2FfPqLFtW/Cp327Zyx+S30kno3wOcA7QMfSCBnv77lpkPAA80Nd0VEacC7wJ+uc1uW4H/0cu6jshCFkU594YyKpO0AO1+fbt9O9xxR/k/xpqtkwu57wY+0O7BzNybmWUMAf0icNocj18PjDVtJ/WjqAXbv4vDA3/GNOwLine truncated
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAFtCAYAAADrr7rKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAgpklEQVR4nO3df5RkZX3n8fe3h2Fgh+7GrGbAAZwJHAwGMiicxGE9iJJJFB3A3bB6jIkTNyu69hnjumadhSSaZUJ+nKhknSirRiHZY/xnF6ejBjAEcyIj0VEmso4QDOgwMkMU6G6zztjQ3/3jVmvRXdXd1VNV91bf9+uce4q69dyq7xQ1n7l1n6eeJzITSVI9DJVdgCSpfwx9SaoRQ1+SasTQl6QaMfQlqUYMfUmqEUNfkmrE0JekGjmu7AL6LSICeDYwVXYtktRFw8C3c5Ff3NYu9CkC/+Gyi5CkHjgNOLhQgzqG/hTAgQMHGBkZKbsWSTpmk5OTnH766bCEKxh1DH0ARkZGDH1JtWNHriTViKEvSTVi6EtSjRj6klQjhr4k1YihL0kVkjO9Xc2w1NCPiDdHxD9ExGRj2xMRL1/kmKsi4usRcSQivhoRl/WrXknqhanHp7jzj7ZzYNdGDt14Ogd2beTOP9rO1OPdnzggylwjNyK2Ak8B/wgE8HrgHcDzM/P/tmh/EfC3wA7gL4HXAv8VeEFm3rvE1xwBJiYmJhynL6l0U49PcejmzfzEj+1n1aqZH+5/8qkhHnzsHE75lT0MP2N4weeYnJxkdHQUYDQzJxdqW+qZfmaOZ+anM/MfM/P+zLwG+B7wwjaHvBX4q8z8w8zcn5m/CXwZGOtXzZLUTXv/9Jp5gQ9w3KoZNv7YfvZ+9Nquvl5lrulHxKqIeA2wFtjTptlm4LNz9t3a2N/ueddExMjsRjEpkSRVwpknjM8L/FnHrZrhzON3d/X1Sg/9iDgvIr4HHAU+CLwqM7/WpvkpwOE5+w439rezA5ho2pxsTVIl5Exy3ND0gm2OG5ruaudu6aEP3AecD/ws8AHgpoh4Xhef/3pgtGk7rYvPLUnLFkPBkzOrF2zzZK4mhqJrr1l66GfmDzLzgczcm5k7gH0U1+5bOQSsm7NvXWN/u+c/mpmTsxvOoy+pQr5xZCtPPtU6ip98aohvHL28q69Xeui3MASsafPYHuDSOfu20L4PQJIq7YI37OTBx86ZF/xPPjXEg989hwt+9bquvl6pUytHxPXAZ4BvUXSwvha4BPiFxuM3Awcb3wAAbgA+FxFvBz4FvAa4EHhjfyuXpO4YfsYw/Moe/u6j13Lm8bs5bmiaJ2dW840fXM4Fv3rdosM1O1X2OP2PUJy5n0rRyfoPwO9n5u2Nx+8EHsrMbU3HXAVcB2ygGN//G5n56Q5e03H6kiorZ7Lja/idjNMvNfTLYOhLWmkG5sdZkqT+MvQlqUYMfUmqEUNfkmrE0JekGjH0JalGDH1JqhFDX5JqxNCXpBox9CWpRgx9SaoRQ1+SasTQl6QaMfQlqUYMfUmqEUNfkmrE0JekGjH0JalGDH1JqhFDX1JrNVs/uy6OK7sASRUyPQX7roGD4zAzDUOrYf1W2LQTVg+XXZ26wNCXVJiegts2w8R+YOZH++/fBYfvgJ/fY/CvAF7ekVTYd838wIfi/sR+2HdtGVWpywx9SYWD48wP/FkzcHB3P6tRjxj6kopO25nphdvMTNu5uwIY+pIgoui0XcjQ6qKdBpqhL6mwfivtI2EI1l/ez2rUI4a+pMKmnTB6DvNjYajYv+m6MqpSlxn6kgqrh4thmWePwdoNcOL64vbsMYdrriCRNeuYiYgRYGJiYoKRkZGyy5GqK9Nr+ANicnKS0dFRgNHMnFyorWf6kloz8FckQ1+SasTQl6QaMfQlqUYMfUmqkVJDPyJ2RMQXI2IqIh6NiFsi4rmLHLMtInLOdqRfNUvSICv7TP/FwC7ghcAWYDVwW0SsXeS4SeDUpu05vSxSklaKUufTz8yXNd+PiG3Ao8AFwN8ufGge6mFpkrQilX2mP9do4/axRdqdFBHfjIgDEfHJiPipdg0jYk1EjMxugD8rlFRblQn9iBgC3gd8PjPvXaDpfcAbgCuA11H8Ge6KiNPatN8BTDRtD3erZkkaNJWZhiEiPgC8HHhRZi45mCNiNbAf+Hhm/maLx9cAa5p2DQMPOw2DpJWik2kYKrFGbkS8H3glcHEngQ+QmdMR8RXgrDaPHwWONr3WsZQqSQOt7CGb0Qj8VwEvzcwHl/Ecq4DzgEe6XZ8krTRln+nvAl5LcX1+KiJOaeyfyMzvA0TEzcDBzNzRuP9bwBeAB4CTgXdQDNn8cH9Ll6TBU3bov7lxe+ec/b8KfKzx32fw9NWanwF8CDgFeBzYC1yUmV/rWZWStEJUpiO3X5xPX9JK43z6UjfV7MRIK1vZl3ekapqegn3XwMFxmJmGodXFwuGbdrpsoAaaoS/NNT0Ft22Gif08rTvp/l1w+A7Xi9VA8/KONNe+a+YHPhT3J/bDvmvLqErqCkNfmuvgOPMDf9YMHNzdz2qkrjL0pWaZxTX8hcxM27mrgWXoS80iik7bhQytLtpJA8jQl+Zav5X2fzWGYP3l/axG6ipDX5pr004YPYf5fz2Giv2briujKqkrDH1prtXDxbDMs8dg7QY4cX1xe/aYwzU18JyGQVpMptfwVWlOwyB1k4GvFcTQl6QaMfQlqUYMfUmqEUNfkmrE0JekGjH0JalGDH1JqhFDX5JqxNCXpBox9CWpRgx9SaoRQ1+SasTQV3lqNsOrVAXHlV2AamZ6CvZdUyw+PjNdLD24fmuxcInz1Es9Z+irf6an4LbNMLEfmPnR/vt3weE7XKBE6gMv76h/9l0zP/ChuD+xH/ZdW0ZVUq0Y+uqfg+PMD/xZM3Bwdz+rkWrJ0Fd/ZBbX8BcyM23nrtRjhr76I6LotF3I0GqXJpR6zNBX/6zfSvuP3BCsv7yf1Ui1ZOirfzbthNFzmP+xGyr2b7qujKqkWjH01T+rh4thmWePwdoNcOL64vbsMYdrSn0SWbOOs4gYASYmJiYYGRkpu5x6y/QavtQFk5OTjI6OAoxm5uRCbUs904+IHRHxxYiYiohHI+KWiHjuEo67KiK+HhFHIuKrEXFZP+pVlxn4Ut+VfXnnxcAu4IXAFmA1cFtErG13QERcBHwc+AjwfOAW4JaIOLfn1UrSgKvU5Z2IeBbwKPDizPzbNm0+AazNzFc27fsCcE9mvmkJr+HlHUkrysBc3mlhtHH72AJtNgOfnbPv1sb+eSJiTUSMzG6AvYWSaqsyoR8RQ8D7gM9n5r0LND0FODxn3+HG/lZ2ABNN28PHVqkkDa7KhD7Ftf1zgdd0+Xmvp/gGMbud1uXnl6SBUYmplSPi/cArgYszc7Ez8UPAujn71jX2z5OZR4GjTa91DJVKUveUMWq57CGb0Qj8VwEvzcwHl3DYHuDSOfu2NPZLUqVNTcH27bBxI5x+enG7fXuxvx9KHb0TEX8CvBa4Ariv6aGJzPx+o83NwMHM3NG4fxHwOeCdwKcoLgf9N+AFi/QFzL6mo3cklWJqCjZvhv37YaZplvGhITjnHNizB4aXMdRkkEbvvJniOvudwCNN26ub2pwBnDp7JzPvoviH4o3APuAXgSuXEviSVKZrrpkf+FDc378fru3DOkKVGqffD57pSyrLxo3w0EPtH9+wAR5cykXuOQbpTF9VUrMTAKmfMmF6kXWEpvuwjlAlRu+oRNNTxdq1B8eLlauGVhfz3m/a6ayXUhdFwOpF1hFa3Yd1hDzTr7PpKbhtM9y/C/7lIfj+weL2/l3F/uk+DSeQamLr1qLTtpWhIbi8D+sIGfp1tu8amNjP/MXKZ4r9+/rQqyTVyM6dxSiducE/O3rnuj6sI2To19nBceYH/qwZOLi7n9VIK97wcDEsc2ys6LRdv764HRtb/nDNTnlNv64yi2v4C5mZdqETqcuGh+GGG4qtdr/IVYkiik7bhQz1oVdJqrEy/noZ+nW2fivtPwJDsL4PvUqS+srQr7NNO2H0HOZ/DIaK/Zv60Kskqa8M/TpbPQw/vwfOHoO1G+DE9cXt2WPFfsfpSyuO0zDoR+y0lQaS0zBoeQx8acUz9CWpRgx9SaoRQ1+SasTQl6QaMfQlqUYMfUmqEUNfkmrE0B90NftxnaRj49TKg8glDiUtk6E/aGaXOJy74tX9u+DwHc6ZI1VQlWY48fLOoHGJQ2kgTE3B9u2wcSOcfnpxu317sb9MTrg2aD65sVi8vJ21G+CKB/tVjaQWpqZg82bYvx9mms7PZtfC7fbSiE64tlJ1ssShpNJcc838wIfi/v79cG2JX8gN/UHiEofSQBgfnx/4s2ZmYPfu/tbTzNAfNC5xKFVaJkwv8oV8usQv5Ib+oHGJQ6nSImD1Il/IV5f4hdzQHzQucShV3tatRadtK0NDcHmJX8gdvTPoqjQAWBLg6B31koEvVc7wcBHsY2OwYQOsX1/cjo11P/A75Zm+JPVYr7+Qe6YvSRVSpS/khr4k1UjXQj8izomIf+rwmIsjYjwivh0RGRFXLtL+kka7udspx1S8JNVEN8/0jwee0+Exa4F9wFs6PO65wKlN26MdHi9JtbTkqZUj4j2LNHlWpy+emZ8BPtN4/k4OfTQzn+j09QaCQzAl9VAn8+m/FbgHaNczfNIxV7N090TEGuBe4F2Z+fk+vnb3uSiKtGJU/bytk9B/AHhvZv55qwcj4nxgbzeKWsAjwJuALwFrgF8D7oyIn83ML7epa02j7axqpaiLokgDb2qqmFlzfLyYV2f16uJXuTt3ljsmv5VOQv9LwAVAy9AHEujpv2+ZeR9wX9OuuyLiTOBtwC+3OWwH8Nu9rOuYLGVRlAtvKKMySUvQ7te3u3bBHXeU/2OsuTrpyH078L52D2bmvswsYwjo3wNnLfD49cBo03ZaP4pasoPjzA/8WTNwsMQ5WCULine truncated
|
||||
"text/plain": [
|
||||
"<Figure size 400x400 with 1 Axes>"
|
||||
]
|
||||
@@ -329,42 +411,46 @@
|
||||
],
|
||||
"source": [
|
||||
"axes = plt.figure(figsize=(4, 4), dpi=100).gca()\n",
|
||||
"historyNt = onp.asarray(history)\n",
|
||||
"updatesNt = onp.asarray(updates) \n",
|
||||
"historyNt = onp.asarray(historyNt)\n",
|
||||
"updatesNt = onp.asarray(updatesNt) \n",
|
||||
"axes.scatter(historyGD[:,0], historyGD[:,1], lw=0.5, color='blue')\n",
|
||||
"axes.scatter(historyNt[:,0], historyNt[:,1], lw=0.5, color='orange')\n",
|
||||
"axes.scatter([0], [0], lw=0.25, color='black') # target at 0,0\n",
|
||||
"axes.scatter([0], [0], lw=0.25, color='black', marker='x') # target at 0,0\n",
|
||||
"axes.set_xlabel('x0')\n",
|
||||
"axes.set_ylabel('x1')"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Not completely surprising: for this simple example we can reliably evaluate the Hessian, and Newtons's method profits from the second order information. It's trajectory is much more diagonal (that would be the ideal, shortest path to the solution), and does not slow down as much as GD."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Physical Gradients\n",
|
||||
"\n",
|
||||
"Now we also use the _inverse physics_, i.e. the inverse of z:\n",
|
||||
"$\\mathbf{z}^{-1}(\\mathbf{x}) = [x_0 \\ x_1^{1/2}]^T$\n",
|
||||
"Now we also use inverse physics, i.e. the inverse of z:\n",
|
||||
"$\\mathbf{z}^{-1}(\\mathbf{x}) = [x_0 \\ x_1^{1/2}]^T$, to compute the _physical gradient_. As a slight look-ahead to the next section, we'll use a Newton's step for $L$, and combine it with the inverse physics function to get an overall update.\n",
|
||||
"\n",
|
||||
"Let’s call intermediate physical result $\\mathbf{y}$, i.e. \n",
|
||||
"$\\mathbf{z}(\\mathbf{x}) = \\mathbf{y}$\n",
|
||||
"\n",
|
||||
"With an update step:\n",
|
||||
"$\n",
|
||||
"\\Delta \\mathbf{x} = \n",
|
||||
"\\mathbf{z}^{-1} \\left( \\mathbf{y} - \\eta_L \n",
|
||||
"This gives an update step:\n",
|
||||
"$\\begin{aligned}\n",
|
||||
"\\Delta \\mathbf{x} &= \n",
|
||||
"\\mathbf{z}^{-1} \\left( \\mathbf{z}(\\mathbf{x}) - \\eta\n",
|
||||
" \\left( \\frac{\\partial^2 L }{ \\partial \\mathbf{z}^2 } \\right)^{-1}\n",
|
||||
" \\frac{\\partial L }{ \\partial \\mathbf{z} }\n",
|
||||
"\\right) - \\mathbf{x}\n",
|
||||
"$\n",
|
||||
"\\end{aligned}$\n",
|
||||
"\n",
|
||||
"\n"
|
||||
"Below, we define our inverse function `fun_z_inv_analytic` (we'll come to a variant below), and then evaluate an optimization with the physical gradient for ten steps:\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 7,
|
||||
"execution_count": 8,
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
@@ -387,24 +473,27 @@
|
||||
"source": [
|
||||
"x = np.asarray([3.,3.])\n",
|
||||
"eta = 0.3\n",
|
||||
"history = [x]; updates = []\n",
|
||||
"historyPG = [x]; historyPGz = []; updatesPG = []\n",
|
||||
"\n",
|
||||
"def fun_z_inv_analytic(y):\n",
|
||||
" return np.array( [y[0], np.power(y[1],0.5)] )\n",
|
||||
"\n",
|
||||
"for i in range(10):\n",
|
||||
" y = fun_z(x)\n",
|
||||
" \n",
|
||||
" # Newton step for L(y)\n",
|
||||
" GL = jax.grad(fun_L)(y)\n",
|
||||
" HL = jax.jacobian(jax.jacobian(fun_L))(y)\n",
|
||||
" # Newton step for L(z)\n",
|
||||
" zForw = fun_z(x)\n",
|
||||
" GL = jax.grad(fun_L)(zForw)\n",
|
||||
" HL = jax.jacobian(jax.jacobian(fun_L))(zForw)\n",
|
||||
" HLinv = np.linalg.inv(HL)\n",
|
||||
" y += -eta * np.matmul( HLinv , GL)\n",
|
||||
" \n",
|
||||
" # step in z space\n",
|
||||
" zBack = zForw -eta * np.matmul( HLinv , GL)\n",
|
||||
" historyPGz.append(zBack)\n",
|
||||
"\n",
|
||||
" # \"inverse physics\" step via z-inverse\n",
|
||||
" x = fun_z_inv_analytic(y)\n",
|
||||
" history.append(x)\n",
|
||||
" updates.append( history[-2] - history[-1] )\n",
|
||||
" x = fun_z_inv_analytic(zBack)\n",
|
||||
" historyPG.append(x)\n",
|
||||
" updatesPG.append( historyPG[-2] - historyPG[-1] )\n",
|
||||
" print( \"PG iter %d: \"%i + format(x) )\n",
|
||||
"\n"
|
||||
]
|
||||
@@ -413,12 +502,14 @@
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"All together now"
|
||||
"Now we obtain $[0.084 \\ 0.504]$ as the final position, with a distance of only 0.51! This is clearly better than both Newton and GD.\n",
|
||||
"\n",
|
||||
"Let's directly visualize how the PGs (in red) fare in comparison to Newton's method (orange) and GD (blue)."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 8,
|
||||
"execution_count": 9,
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
@@ -427,13 +518,13 @@
|
||||
"Text(0, 0.5, 'x1')"
|
||||
]
|
||||
},
|
||||
"execution_count": 8,
|
||||
"execution_count": 9,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAFtCAYAAADrr7rKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAjEElEQVR4nO3df5RcZZ3n8fe3Ok0Om3Q3I2ICncQEOcEg2CgcNVkOIAz+QBKIK6tHnJ3AuopLn6jrOks2jMcZiZlxziKZpdXxxyRmcBj+IqQVB2QY2F0To6AJohFkTSRpSECQ7taF0KS++8dThZXqquqq7qq6z637eZ1zT1G3nlv1pKj+1K17n/t9zN0REZFsyCXdARERaR+FvohIhij0RUQyRKEvIpIhCn0RkQxR6IuIZIhCX0QkQxT6IiIZMivpDrSbmRlwCjCedF9ERJqoB3jSp7jiNnOhTwj8g0l3QkSkBRYAI7UaZDH0xwEOHDhAb29v0n0REZmxsbExFi5cCHUcwchi6APQ29ur0BeRzNGJXBGRDFHoi4hkiEJfRCRDFPoiIhmi0BcRyRCFvohIRPxovqXPn2jom9nHzOxhMxsrLDvN7N1TbHOlmf3CzF40s5+a2aXt6q+ISCuMP/EkB88b4OUTZ5F/dTcvnziLg+cNMP7Ek01/raT39A8C1wPnAOcC9wF3mtkbKjU2sxXAbcA3gDcB24BtZnZmW3orItJk4088yXFvPZX+HQ8z67mjdD2fZ9ZzR+nf+TDHvfXUpge/xTYxupk9B3za3b9R4bHbgTnuflnJuh8Au9392jqfvxcYHR0d1cVZIpK4g+cN0L/jYaxCFHsORlYMsOB/7675HGNjY/T19QH0uftYrbZJ7+m/wsy6zOwDwBxgZ5Vmy4F7y9bdXVhf7Xlnm1lvcSEUJRIRicL8vT+rGPgAlof5P/9ZU18v8dA3s7PM7HfAEeArwGp3/3mV5vOBw2XrDhfWV7MOGC1ZVGxNRKLgR/NYvvbRFsvnm3pyN/HQBx4FzgbeCnwZ+KaZndHE598I9JUsC5r43CIi02ZdOTxnNdt4zrCu5kV14qHv7i+5++Pu/pC7rwP2AB+v0vwQMK9s3bzC+mrPf8Tdx4oLqqMvIhE5tOwNeJUk9hwcOqO541QSD/0KcsDsKo/tBC4uW3cJ1c8BiIhEre8fv8tLr5k9Kfg9By+9ZjZ937qrqa+X9Dj9jWZ2vpktLhzb3whcCHyr8PjWwrqiTcC7zOxTZvZ6M/ssYajnLe3uu4hIM/QsOoWXdv2KkRUDvPyqWRw9IcfLr5rFyIoBXtr1K3oWndLU10u6nv5rgK3AyYSTrA8D73T37xUeXwS8cgbD3XeY2QeBG4HPA78ErnD3R9raaxGRJupZdAo9hWGZfjSPdeVadvIxunH6raZx+iLSaVI5Tl9ERFpPoS8ikiEKfRGRDFHoi4hkiEJfRCRDFPoiIhmi0BcRyRCFvohIhij0RUQyRKEvIpIhCn0RkQxR6IuIZIhCX0QkQxT6IiIZotAXEckQhb6ISIYo9EVEMkShLyKSIQp9EZEMUeiLSGUZmz87K2Yl3QERicjEOOxZDyPDkJ+AXDf0r4SBDdDdk3TvpAkU+iISTIzDPcthdC+Q/8P6x4bg8H3wjp0K/g6gwzsiEuxZPznwIdwf3Qt7bkiiV9JkCn0RCUaGmRz4RXkY2d7O3kiLKPRFJJy0zU/UbpOf0MndDqDQFxEwCydta8l1h3aSagp9EQn6V1I9EnLQv6qdvZEWUeiLSDCwAfqWMTkWcmH9wI1J9EqaTKEvIkF3TxiWuXQQ5iyG4/vD7dJBDdfsIOYZOzFjZr3A6OjoKL29vUl3RyRe7jqGnxJjY2P09fUB9Ln7WK222tMXkcoU+B1JoS8ikiEKfRGRDFHoi8jMZOy8YNop9EWkcePjsHYtLFkCCxeG27Vrw3qJWqKhb2brzOxHZjZuZk+b2TYzO32KbdaYmZctL7arzyKZNz4Oy5fD0BDs3w8jI+F2aCisV/BHLek9/QuAIeBtwCVAN3CPmc2ZYrsx4OSS5bWt7KSIlFi/HvbuhXxZcbZ8Pqy/QdU4Y5Zo6Lv7u9x9i7v/zN33AGuARcA5U2/qh0qWwy3vrIgEw8OTA78on4ftqsYZs6T39Mv1FW6fm6LdXDP7tZkdMLM7zewN1Rqa2Wwz6y0ugC4rFJkud5iYohrnhKpxxiya0DezHHAz8H13f6RG00eBa4DLgQ8R/g07zGxBlfbrgNGS5WCz+iySOWbQPUU1zm5V44xZNKFPOLZ/JvCBWo3cfae7b3X33e7+APBe4Bngo1U22Uj4BVFcqn05iEg9Vq6EXJXoyOVglapxxiyK0DezW4DLgLe7e0N74u4+AfwEOK3K40fcfay4ABpaIDITGzbAsmWTgz+XC+tvVDXOmCU9ZNMKgb8auMjd903jObqAs4Cnmt0/Eamgpwd27oTBQVi8GPr7w+3gYFjfo9NmMUu0yqaZfQn4IOH4/KMlD426+wuFNluBEXdfV7j/GeAHwOPACcCngSuAc9z953W8pqpsijSTqnEmrpEqm7Pa06WqPla4vb9s/dXAlsJ/L+LY2Zr/CPgaMB/4LfAQsKKewBeRFlDgp4rq6YuIpJzq6Ys0U8Z2jKSzJX14RyROE+OwZz2MDEN+AnLdYeLwgQ2aNlBSTaEvUm5iHO5ZDqN7OeZ00mNDcPg+zRcrqabDOyLl9qyfHPgQ7o/uhT0qKCbppdAXKTcyzOTAL8rDiAqKSXop9EVKuYdj+LXkVVBM0kuhL1LKLJy0rSWngmKSXgp9kXL9K6n+p5GDfhUUk/RS6IuUG9gAfcuY/OeRC+sHVFBM0kuhL1KuuycMy1w6CHMWw/H94XbpoIZrSuqpDIPIVFRQTCKnMgwizaTAlw6i0BcRyRCFvohIhij0RUQyRKEvIpIhCn0RkQxR6IuIZIhCX0QkQxT6IiIZotAXEckQhb6ISIYo9EVEMkShLyKSIQp9SU7GKryKxGBW0h2QjJkYhz3rw+Tj+Ykw9WD/yjBxSafUqVcpZomY9vSlfSbG4Z7l8NgQ/H4/vDASbh8bCusnxpPu4fSNj8PatbBkCSxcGG7Xrg3rRSKiSVSkfR5cGwKefIUHc2FmqnM3tbtXMzc+DsuXw969kC/5t+VysGwZ7NwJPR3yK0aipElUJE4jw1QOfML6ke3t7E3zrF8/OfAh3N+7F264IZl+iVSg0Jf2cA/H8GvJT6Tz5O7w8OTAL8rnYXtKv8ykIyn0pT3MwknbWnLd6TsB6g4TU3yZTaT0y0w6kkJf2qd/JdU/cjnoX9XO3jSHGXRP8WXWncIvM+lYCn1pn4EN0LeMyR+7XFg/cGMSvZq5lSvDSdtKcjlYlcIvM+lYGr0j7TUxDntuCCdtXxmnvyoEflrH6Wv0jiSskdE7Cn1JTiddxDQ+HkbpbN8ejuF3d4c9/BtvVOBLy6Um9M1sHfBe4PXAC8AO4L+5+6NTbHcl8DlgMfDLwjZ31fmaCn1prU76MpNUSNM4/QuAIeBtwCVAN3CPmc2ptoGZrQBuA74BvAnYBmwzszNb3luReijwJWJRHd4xs5OAp4EL3P1/VWlzOzDH3S8rWfcDYLe7X1vHa2hPX0Q6Spr29Mv1FW6fq9FmOXBv2bq7C+snMbPZZtZbXAAdYBWRzIom9M0sB9wMfN/dH6nRdD5wuGzd4cL6StYBoyXLwZn1VEQkvaIJfcKx/TOBDzT5eTcSfkEUlwVNfn4RkdSIop6+md0CXAac7+5T7YkfAuaVrZtXWD+Jux8BjpS81gx6KiLSPEkM9Ep0T9+CW4DVwEXuvq+OzXYCF5etu6SwXkQkaklPvZD0OP0vAR8ELgdKx+aPuvsLhTZbgRF3X1e4vwJ4ALge+A7hcNB/B948xbmA4mtq9I6IJKJVF2+nafTOxwjH2e8HnipZ3l/SZhFwcvGOu+8gfFF8BNgDvA+4op7AFxFJUgxTL0Q1Tr8dtKcvIklZsgT276/++OLFsK+eg9xl0rSnLzHJ2A6ASDvFMvVCFKN3JEET47BnfZjK8JWqlytDGeS0Vr0UiVAsUy9oTz/LJsbhnuVhsvLf74cXRsLtY0Nh/USbhhOIZEQMUy8o9LNsz3oY3cvkycrzYf0eTegt0kwbNoRROuXBXxy9c2Mb5hFS6GfZyDCTA78oHyY6EZGm6ekJwzIHB8NJ2/7+cDs42L65dnRMP6vcwzH8WvITqg0v0mQ9PbBpU1gyd0WuJMgsnLStJacJvUVaKYk/L4V+lvWvpPpHIBfmrhWRjqLQz7KBDdC3jMkfg1xYP9CGs0oi0lYK/Szr7oF37ISlgzBnMRzfH26XDob1Gqcv0nFUhkH+IA0nbdPQR5E2UxkGmZ5YwzTpWrQiHUR7+hK3VtWiFekg2tOXzhFDLVqRDqLQl7gND08O/KJ8HrbrqmGRRij0JV6x1KIV6SAKfYlXLLVoRTqIQl/iFkMtWpEOotCXuMVQi1akgyj0067Tj2fHUItWpINonH4aZXmKQ12RKzJJI+P0VU8/bYpTHJbPePXYEBy+r/Nr5ijwJYVi2lfR4Z200RSHIqkQa/UQHd5JmzuXhMnLq5mzGC7f167eiEgF7a4eojIMnaqRKQ5FJDExVw9R6KeJpjgUSYWYq4co9NNGUxyKRC326iEK/bTRFIciUYu9eohCP200xaFI9GKuHqLRO2kX0wBgEQE0ekdaSYEvEp2Yq4doT19EpMVa/YNce/oiIhGJ6Qe5Ql9EJEOaFvpmtszMftXgNueb2bCZPWlmbmZXTNH+wkK78mX+jDovIpIRzdzTPw54bYPbzAH2ANc1uN3pwMkly9MNbi8ikkl1l1Y2s5umaHJSoy/u7t8Fvlt4/kY2fdrdn2/09VIhTUMw09RXEQEaq6f/cWA3UO3M8NwZ96Z+u81sNvAI8Fl3/34bX7v50jQpyvh4qCY1PByuJe/uDleibNigWaxEiH9fqO4hm2bLine truncated
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAFtCAYAAADrr7rKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAi1UlEQVR4nO3df5RcZZ3n8fe3kiawSXczjppAk5ggJxgFG4SjJstBhI0iEn44MnpcZ8m6M4JLn6jrOks2jKsjMTPOWSSztJpVR8joOPwX0qMOP4aB3TUtI9E0ohFkTSRpSFCR7tYloUl/94+nCivVdauruqvqPrfu53XOPUXdem7Vk6L6U7fufe73MXdHRETyoZB2B0REpH0U+iIiOaLQFxHJEYW+iEiOKPRFRHJEoS8ikiMKfRGRHFHoi4jkyPy0O9BuZmbAqcBE2n0REWmibuApn+GK29yFPiHwD6bdCRGRFjgNGK3VII+hPwFw4MABenp60u6LiMicjY+Ps3TpUqjjCEYeQx+Anp4ehb6I5I5O5IqI5IhCX0QkRxT6IiI5otAXEckRhb6ISI4o9EVEIuLHplr6/KmGvpl9yMweMbPx4jJsZu+YYZtrzOwnZnbEzH5oZpe1q78iIq0w8eRTHLygnxd/fz5TL+/ixd+fz8EL+pl48qmmv1bae/oHgRuB84DzgfuBu8zsddUam9ka4BvAV4BzgR3ADjM7qy29FRFpsoknn+KEN51O365HmP/sMeY9N8X8Z4/RN/wIJ7zp9KYHv8U2MbqZPQt83N2/UuWxO4GF7n552brvAnvc/fo6n78HGBsbG9PFWSKSuoMX9NO36xGsShR7AUbX9HPa/95T8znGx8fp7e0F6HX38Vpt097Tf4mZzTOz9wILgeGEZquB+yrW3V1cn/S8C8ysp7QQihKJiERhyd4fVQ18AJuCJT/+UVNfL/XQN7Ozzew3wFHgi8DV7v7jhOZLgMMV6w4X1yfZCIyVLSq2JiJR8GNT2FTtoy02NdXUk7uphz7wGHAO8CbgC8AdZvbaJj7/FqC3bDmtic8tIjJrNq+AF6xmGy8YNq95UZ166Lv7C+7+hLvvdveNwAjw4YTmh4DFFesWF9cnPf9Rdx8vLaiOvohE5NCq1+EJSewFOPTa5o5TST30qygACxIeGwYuqVi3luRzACIiUev9u2/zwisXTAt+L8ALr1xA79e/1dTXS3uc/hYzu9DMlheP7W8BLgK+Xnx8e3FdyVbgUjP7mJm9xsw+SRjqeVu7+y4i0gzdy07lhYd+xuiafl582XyOnVzgxZfNZ3RNPy889DO6l53a1NdLu57+K4HtwCmEk6yPAG9393uLjy8DXjqD4e67zOx9wM3AZ4CfAle5+6Nt7bWISBN1LzuV7uKwTD82hc0rtOzkY3Tj9FtN4/RFpNNkcpy+iIi0nkJfRCRHFPoiIjmi0BcRyRGFvohIjij0RURyRKEvIpIjCn0RkRxR6IuI5IhCX0QkRxT6IiI5otAXEckRhb6ISI4o9EVEckShLyKSIwp9EZEcUeiLiOSIQl9EJEcU+iIiOaLQF5HqcjZ/dl7MT7sDIhKRyQkY2QSjQzA1CYUu6FsH/Zuhqzvt3kkTKPRFJJicgHtWw9heYOp36x8fhMP3w9uGFfwdQId3RCQY2TQ98CHcH9sLIzel0StpMoW+iASjQ0wP/JIpGN3Zzt5Iiyj0RSSctJ2arN1malIndzuAQl9EwCyctK2l0BXaSaYp9EUk6FtHciQUoO+KdvZGWkShLyJB/2boXcX0WCiE9f03p9EraTKFvogEXd1hWObKAVi4HE7qC7crBzRcs4OY5+zEjJn1AGNjY2P09PSk3R2ReLnrGH5GjI+P09vbC9Dr7uO12mpPX0SqU+B3JIW+iEiOKPRFRHJEoS8ic5Oz84JZp9AXkcZNTMCGDbBiBSxdGm43bAjrJWqphr6ZbTSz75nZhJk9Y2Y7zOzMGbZZb2ZesRxpV59Fcm9iAlavhsFB2L8fRkfD7eBgWK/gj1rae/pvAQaBNwNrgS7gHjNbOMN248ApZcurWtlJESmzaRPs3QtTFcXZpqbC+ptUjTNmqYa+u1/q7re7+4/cfQRYDywDzpt5Uz9UthxueWdFJBgamh74JVNTsFPVOGOW9p5+pd7i7bMztFtkZj83swNmdpeZvS6poZktMLOe0gLoskKR2XKHyRmqcU6qGmfMogl9MysAtwLfcfdHazR9DPgAcCXwfsK/YZeZnZbQfiMwVrYcbFafRXLHDLpmqMbZpWqcMYsm9AnH9s8C3lurkbsPu/t2d9/j7g8C7wJ+AVyXsMkWwi+I0pL05SAi9Vi3DgoJ0VEowBWqxhmzKELfzG4DLgfe6u4N7Ym7+yTwA+CMhMePuvt4aQE0tEBkLjZvhlWrpgd/oRDW36xqnDFLe8imFQP/auBid983i+eYB5wNPN3s/olIFd3dMDwMAwOwfDn09YXbgYGwvlunzWKWapVNM/s88D7C8fnHyh4ac/fni222A6PuvrF4/xPAd4EngJOBjwNXAee5+4/reE1V2RRpJlXjTF0jVTbnt6dLiT5UvH2gYv2/B24v/vcyjp+t+feALwFLgF8Du4E19QS+iLSAAj9TVE9fRCTjVE9fpJlytmMknS3twzsicZqcgJFNMDoEU5NQ6AoTh/dv1rSBkmkKfZFKkxNwz2oY28txp5MeH4TD92u+WMk0Hd4RqTSyaXrgQ7g/thdGVFBMskuhL1JpdIjpgV8yBaMqKCbZpdAXKecejuHXMqWCYpJdCn2RcmbhpG0tBRUUk+xS6ItU6ltH8p9GAfpUUEyyS6EvUql/M/SuYvqfRyGs71dBMckuhb5Ipa7uMCxz5QAsXA4n9YXblQMarimZpzIMIjNRQTGJnMowiDSTAl86iEJfRCRHFPoiIjmi0BcRyRGFvohIjij0RURyRKEvIpIjCn0RkRxR6IuI5IhCX0QkRxT6IiI5otAXEckRhb6ISI4o9CU9OavwKhKD+Wl3QHJmcgJGNoXJx6cmw9SDfevCxCWdUqdepZglYtrTl/aZnIB7VsPjg/Db/fD8aLh9fDCsn5xIu4ezNzEBGzbAihWwdGm43bAhrBeJiCZRkfZ5eEMIeKaqPFgIM1Odv7XdvZq7iQlYvRr27oWpsn9boQCrVsHwMHR3yK8YiZImUZE4jQ5RPfAJ60d3trM3zbNp0/TAh3B/71646aZ0+iVShUJf2sM9HMOvZWoymyd3h4amB37J1BTszOiXmXQkhb60h1k4aVtLoSt7J0DdYXKGL7PJjH6ZSUdS6Ev79K0j+SNXgL4r2tmb5jCDrhm+zLoy+GUmHUuhL+3Tvxl6VzH9Y1cI6/tvTqNXc7duXThpW02hAFdk8MtMOpZG70h7TU7AyE3hpO1L4/SvCIGf1XH6Gr0jKWtk9I5CX9LTSRcxTUyEUTo7d4Zj+F1dYQ//5psV+NJymQl9M9sIvAt4DfA8sAv4L+7+2AzbXQN8GlgO/LS4zbfqfE2FvrRWJ32ZSSZkaZz+W4BB4M3AWqALuMfMFiZtYGZrgG8AXwHOBXYAO8zsrJb3VqQeCnyJWFSHd8zsFcAzwFvc/X8ltLkTWOjul5et+y6wx92vr+M1tKcvIh0lS3v6lXqLt8/WaLMauK9i3d3F9dOY2QIz6yktgA6wikhuRRP6ZlYAbgW+4+6P1mi6BDhcse5wcX01G4GxsuXg3HoqIpJd0YQ+4dj+WcB7m/y8Wwi/IErLaU1+fhGRzIiinr6Z3QZcDlzo7jPtiR8CFlesW1xcP427HwWOlr3WHHoqItI8aQz0SnVP34LbgKuBi919Xx2bDQOXVKxbW1wvIhK1tKdeSHuc/ueB9wFXAuVj88fc/flim+3AqLtvLN5fAzwI3Ah8k3A46L8Cb5jhXEDpNTV6R0RS0aqLt7M0eudDhOPsDwBPly3vKWuzDDildMfddxG+KD4IjADvBq6qJ/BFRNIUw9QLUY3Tbwft6YtIWlasgP37kx9fvhz21XOQu0KW9vQlJjnbARBpp1imXohi9I6kaHICRjaFqQxfqnq5LpRBzmrVS5EIxTL1gvb082xyAu5ZHSYr/+1+eH403D4+GNZPtmk4gUhOxDD1gkI/z0Y2wdhepk9WPhXWj2hCb5Fm2rw5jNKpDP7S6J2b2zCPkEI/z0aHmB74JVNhohMRaZru7jAsc2AgnLTt6wu3AwPtm2tHx/Tzyj0cw69lalK14UWarLsbtm4NS+6uyJUUmYWTtrUUNKG3SCul8eel0M+zvnUkfwQKYe5aEekoCv08698MvauY/jEohPX9bTirJCJtpdDPs65ueNswrByAhcvhpL5wu3IgrNc4fZGOozIM8jtZOGmbhT6KtJnKMMjsxBqmadeiFekg2tOXuLWqFq1IB9GevnSOGGrRinQQhb7EbWhoeuCXTE3BTl01LNIIhb7EK5ZatCIdRKEv8YqlFq1IB1HoS9xiqEUr0kEU+hK3GGrRinQQhX7Wdfrx7Bhq0Yp0EI3Tz6I8T3GoK3JFpmlknL7q6WdNaYrDyhmvHh+Ew/d3fs0cBb5kUEz7Kjq8kzWa4lAkE2KtHqLDO1lz14oweXmShcvhyn3t6o2IVNHu6iEqw9CpGpniUERSE3P1EIV+lmiKQ5FMiLl6iEI/azTFoUjUYq8eotDPGk1xKBK12KuHKPSzRlMcikQv5uohGr2TdTENABYRQKN3pJUU+CLRibl6iPb0RURarNU/yLWnLyISkZh+kCv0RURypGmhb2arzOxnDW5zoZkNmdlTZuZmdtUM7S8qtqtclsyp8yIiOdHMPf0TgFc1uM1CYAS4ocHtzgROKVueaXB7EZFcqru0spndMkOTVzT64u7+beDbxedvZNNn3P25Rl8vE7I0BDNLfRURoLF6+h8G9gBJZ4YXzbk39dtjZguAR4FPuvt32vjazZelSVEmJkI1qaGhcC15V1e4EmXzZs1iJUL8+0J1D9kLine truncated
|
||||
"text/plain": [
|
||||
"<Figure size 400x400 with 1 Axes>"
|
||||
]
|
||||
@@ -445,13 +536,14 @@
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"historyPG = onp.asarray(historyPG)\n",
|
||||
"updatesPG = onp.asarray(updatesPG) \n",
|
||||
"\n",
|
||||
"axes = plt.figure(figsize=(4, 4), dpi=100).gca()\n",
|
||||
"historyPG = onp.asarray(history)\n",
|
||||
"updatesPG = onp.asarray(updates) \n",
|
||||
"axes.scatter(historyGD[:,0], historyGD[:,1], lw=0.5, color='blue')\n",
|
||||
"axes.scatter(historyNt[:,0], historyNt[:,1], lw=0.5, color='orange')\n",
|
||||
"axes.scatter(historyPG[:,0], historyPG[:,1], lw=0.5, color='red')\n",
|
||||
"axes.scatter([0], [0], lw=0.25, color='black') # target at 0,0\n",
|
||||
"axes.scatter([0], [0], lw=0.25, color='black', marker='x') # target at 0,0\n",
|
||||
"axes.set_xlabel('x0')\n",
|
||||
"axes.set_ylabel('x1')"
|
||||
]
|
||||
@@ -460,16 +552,14 @@
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"So, PG does best...\n",
|
||||
"Also most _diagonal_, i.e., direct path to target\n",
|
||||
"This illustrates that the PG variant does significantly better than Newton's method. It yields a trajectory that is better aligned with the _diagonal_, and its final state is closer to the origin. A key ingredient here is the inverse function for $\\mathbf{z}$, which yielded higher order terms than the second-order approximation for Newton's method. Despite the simplicity of the problem, Newton's method has problems finding the right search direction. For PG, on the other hand, the higher order information yields an improved direction for the optimization.\n",
|
||||
"\n",
|
||||
"Difference also shows in first update steps: compare...\n",
|
||||
"GD mostly \"down\" along x1, Newton mostly along x0, PG is nicely diagonal."
|
||||
"This difference also shows in first update step for each method: below we measure how well it is aligned with the diagonal."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 9,
|
||||
"execution_count": 10,
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
@@ -495,40 +585,133 @@
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Approximate Inversions\n",
|
||||
"The largest value of 1.356 for PG confirms what we've seen above: the PG gradient was the closest one to the diagonal direction from our starting point to the origin."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"\n",
|
||||
"We can use BFGS if inverse of z is not readily available\n",
|
||||
"---\n",
|
||||
"\n",
|
||||
"Just a bit ugly for calling `fmin_l_bfgs_b()` from scipy:"
|
||||
"<br>\n",
|
||||
"\n",
|
||||
"## Z Space\n",
|
||||
"\n",
|
||||
"To understand the behavior and differences of the methods here, it's important to keep in mind that we're not dealing with a black box that maps between $\\mathbf{x}$ and $L$, but rather there are spaces inbetween that matter. In our case, we only have a single $\\mathbf{z}$ space, but for DL settings, we might have a large number of latent spaces, over which we have a certain amount of control.\n",
|
||||
"\n",
|
||||
"We will return to NNs soon, but for now let's focus on $\\mathbf{z}$. One first thing to note is that for PG, we explicitly map from $L$ to $\\mathbf{z}$, and then continue with a mapping to $\\mathbf{x}$. Thus we already obtained the trajectory in $\\mathbf{z}$ space, and not conincidentally, we already stored it in the `historyPGz` list above.\n",
|
||||
"\n",
|
||||
"Let's directly take a look what PG did in $\\mathbf{z}$ space:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 10,
|
||||
"execution_count": 14,
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"BFGS optimization test run, find x such that y=[2,2]:\n"
|
||||
]
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"Text(0, 0.5, 'z1')"
|
||||
]
|
||||
},
|
||||
"execution_count": 14,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAGDCAYAAAAyH0uIAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAcxklEQVR4nO3de5zddX3n8deHBCY+IDOIA1IbLTha0ha7XIJotgWMCtq1raigZd0lG1svFbdiH7VS7+6yrFsfDMVgTSPetongrpU+xFWLlVxQ4gVvg3dToFCuI8wZIjCJ8bN//M6Qk2FmMjM5Z37nm3k9H4/f4+R8z/d3fp/z4/DON9/zu0RmIkkqz0F1FyBJmhsDXJIKZYBLUqEMcEkqlAEuSYUywCWpUAa4JBXKAJekQhngklQoA1ySCmWAS1KhDHBJKpQBLkmFMsBVlIg4JiJyqmUf6y6NiMsi4taIGIuIeyPiuog4qaXPpoi4OSJOjoivRMTDEXFLRLx2wnsdEhHviYibIqIRET+PiK0R8ZxJtntQRPxZRAxFxCMRcV9EfD4iVkzo98rm+z0cEfdHxFUR8eT93Wc6cIWXk1VJIuJQ4OwJzQcDg8DOzDxqmnU3AC8D1gLfB54A/A5wdWZuaPbZBDwdWAx8EvgxcG6z36sy88PNfv3Ad4FPAD8BlgKvAp4KPDMzv92y3Y8Aq4HPAV9ovvfvAl/MzLXNPm8F/ltzm5uBI4E3ADuAEzNzZBa7SQtFZrq4FL0AVwC/AJ6zj34jwNp99NkEJPCmlrZDgG8B9wAHN9sWAYdMWPdw4G7gypa25zTf728m2db4AOrXmvX/1YTXjwd2TWx3cRlfnEJR0SLiPwN/Crw5M6/fR/cR4NSIeNI++v0CWDf+JDN3Np8fBZzcbNvdbB+fIjmCamT9DeCklvd6KVWAv3viRjJz/J+/L6GazvxkRPSPL1R/GfyE6i8B6TEMcBUrIk4APgh8IjMvncEqb6Ya1d4eEV+LiHdFxFMn6XdnZv58QtuPm4/HtGz//Ij4LvAI8DPgPuA/AH0t6w003+/+aep6OhBUYX3fhOU3qP7ikB5jcd0FSHMREY8HPkUVrH88k3Uy85MRsZVqDv1M4C+Av4yIl2Tm52a5/VcCHwWuAf4auBfYDVxEFdqzcRDVKP2FzfeYaMcs308LhAGu4kTEQcAGqjnn52XmQzNdNzPvAj4AfCAijgK+CbyV6gfGcU+KiEMnjMJ/vfl4a/PxZcC/AC9pmQohIiZOlWwHzoqII6YZhW+nGoHfkpk/nqKP9BhOoahE7wTOAv4oM2+ZyQoRsSgiWqc2yMx7gTuBngndFwOvaVn3kObz+4Cbms3jI+Vo6Xcq8OwJ7/WpZp93TlLT+Lr/0Hy/d7a0PdonIp6wr8+nhckRuIoSEc8A3g5sAY5qTmU8KjP/fopVlwJ3RMT/Bb5DNS3xPOAU4M8n9L2TamrlGKopmpcDJwCvzsxdzT7XUv34+OmI+CxwLPBaqsMTD2up5/qI+N/Af42IpwOfpxo4/S5wPdVRMdsj4m3AJcAxEXEN8GDzPc8G/g543wx3kRYQjwNXUSLiDKrgm1RmxmTtzVH0f6ea+34qVYj+FFiXmX/b0m8T0A+cD7wfOJHq8MG/zswrWvoF8BaqkfnRVMH9duAc4IzMPKal7yLgQvYcJ96gOlrlbZn5zZZ+L2n2O7HZdDvwz8DlTq1oMga41GI8wDPz+LprkfbFOXBJKpQBLkmFMsAlqVDOgUtSoRyBS1KhDHBJKlTRJ/I0j8V9EtVJD5J0oFhKdRG0aee4iw5wqvC+o+4iJKkDlgH/Nl2H0gP8QYDbb7+d3t7eumuRpP02OjrKk5/8ZJjBzELpAQ5Ab2+vAS5pwfFHTEkqlAEuSYUywCWpUAa4JBXKAJekQhngktQpHb7WlAEuSZ2wezesWFE9dogBLkmdsHUrDA3BDTd0bBMHxIk8ktQVhoeh0aj+PDgIu3ZVj8uWVW19fdDf37bNOQKXpHa5+GJYvhxWroQtW6q2zZur58uXV6+3kQEuSe0yOAjr1sGiRTAyUrWNjMDixbB+ffV6GxngktROa9bAqlV7t61aBatXt31TzoFLUjtlVtMm/f0wMADbt8OmTVV7RFs35QhcktppaAh27ICNG2HbNtiwoXo+NNT2TRV9U+OI6AUajUbDy8lK6g5jY9XSmkmNBixZAj09+1x9dHSUvr4+gL7MHJ2ur1MoktROPT2PDeoqkNvOKRRJKpQBLkmFMsAlqVAGuCQVqvYAj4hfjYi/j4ifRcTDETEUESvqrkuSul2tR6FExOOBLwPXAy8E7gOeDjxQZ12SVIK6DyP8S+D2zPwvLW231FWMJJWk7imUPwC+ERH/JyLujYhvRcSfTNU5Inoiond8AZbOX6mS1F3qDvCnAq8DfgKcBfwtcHlEnD9F/4uARstyx3wUKUndqNZT6SNiJ/CNzFzZ0nY5cEpmPnuS/j1A6ylOS4E7PJVe0oGipFPp7wK+P6HtB8BLJ+ucmWPA2PjzaPOVvSSpJHVPoXwZOG5C268Dt9VQiyQVpe4AHwSeFRF/FRFPi4jzgFcDV9RclyR1vVoDPDO/DpwN/BFwM/B24I2ZuaHOuiSpBHXPgZOZ1wLX1l2HJJWm7ikUSdIcGeCSVCgDXJIKZYBLUqEMcEkqlAEuSYUywCWpUAa4JBXKAJekQhngklQoA1ySCmWAS1KhDHBJKpQBLkmFMsAlqVAGuCQVygCXpEIZ4JJUKANckgplgEtSoQxwSSqUAS5JhTLAJalQBrgkFcoAl6RCGeCSVCgDXJIKZYBLUqEMcEkqlAEuSYUywCWpUAa4JBWq1gCPiHdFRE5YflhnTZJUisV1FwB8D3hey/Nf1FWIJJWkGwL8F5l5d91FSFJpumEO/OkRcWdE/EtEbIiIp0zVMSJ6IqJ3fAGWzmOdktRV6g7wrwKrgRcArwOOBbZGxFTBfBHQaFnumIcaJakrRWbWXcOjIuJw4DbgTZl55SSv9wA9LU1LgTsajQa9vb3zU6QkddDo6Ch9fX0AfZk5Ol3fbpgDf1RmjkTEj4GnTfH6GDA2/jwi5qs0Seo6dU+h7CUiDgMGgLvqrkWSul3dx4G/LyJOj4hjImIl8GlgN/CJOuuSpBLUPYWyjCqsnwDcB9wAPCsz76u1KkkqQK0BnpmvqHP7klSyrpoDlyTNnAEuSYUywCWpUAa4JBXKAJekQhngklQoA1ySCmWAS1KhDHBJKpQBLkmFMsAlqVAGuCQVygCXpEIZ4JJUKANc0h5ddI9c7ZsBLqmyezesWFE9qggGuKTK1q0wNAQ33FB3JZqhum+pJqlOw8PQaFR/HhyEXbuqx2XLqra+Pujvr68+TcsRuLSQXXwxLF8OK1fCli1V2+bN1fPly6vX1bUMcGkhGxyEdetg0SIYGanaRkZg8WJYv756XV3LAJcWujVrYNWqvdtWrYLVq2spRzPnHLi00GVW0yb9/TAwANu3w6ZNVXtE3dVpGo7ApYVuaAh27ICNG2HbNtiwoXo+NFR3ZdqHyIIP3I+IXqDRaDTo7e2tuxypTGNj1dL6/1CjAUuWQE9PfXUtUKOjo/T19QH0ZebodH2dQpEWup6exwZ1FSDqck6hSFKhDHBJKpQBLkmFMsAlqVAGuCQVygCXpEIZ4JJUKANckgrVNQEeEW+JiIyIy+quRZJK0BUBHhGnAK8Bvlt3LZJUitoDPCIOAzYAfwI8UHM5klSM2gMcuAL4bGZ+cV8dI6InInrHF2Bp58uTpO5U68WsIuIVwEnAKTNc5SLgnZ2rSJLKUdsIPCKeDPwN8B8z85EZrnYJ0NeyLOtQeZLU9eocgZ8MHAV8M/bc9WMRcFpEXAD0ZObu1hUycwwYG38e3i1E0gJWZ4D/M/CMCW0fAX4IvHdieEuS9lZbgGfmg8DNrW0R8XPgZ5l58+RrSZLGdcNRKJKkOeiqW6pl5hl11yBJpXAELkmFMsAlqVAGuCQVygCXpEIZ4JJUKANckgplgEtSoQxwSSqUAS5JhTLAJalQBrgkFcoAl6RCGeCSVCgDXJIKZYBLUqEMcEkqlAEuSYUywCWpUAa4JBXKAJekQhngklQoA1ySCmWAS1KhDHBJKpQBLkmFamuAR8S/i4jd7XxPSdLkOjECjw68pyRpgsWz6RwR/7CPLn1Azr0cqctkQjgmUXea7Qj894ElQGOKZUdbq5PqtHs3rFhRPUpdaFYjcOAHwKcy88rJXoyIE4AX7W9RUlfYuhWGhuCGG+D00+uuRnqM2Qb4TcBJwKQBDowB/7pfFUl1Gh6GRqP68+Ag7NpVPS5bVrX19UF/f331SS0ic+ZT1hHRAyzKzIc6V9LMRUQv0Gg0GvT29tZdjg4EF14Ia9fCEUfAzp0wMgKHHw6HHAL33w8XXFAFutQho6Oj9PX1AfRl5uh0fWc1B56ZY5n5UER8KSLeOfH1iHh8RHxppu8XEa+LiO9GxGhzuTEiXjibmqS2GhyEdetg0aIqvKF6XLwY1q83vNVV5noY4RnABRFxTUQc2tJ+CDCbycI7gLcAJwMrgC8B/xgRvzXHuqT9t2YNrFq1d9uqVbB6dS3lSFPZn+PAnwccDWyLiGPm8gaZ+ZnM/H+Z+ZPM/HFmvpXqSJZn7Udd0v7JhM2bq7nuU0+tHjdtqtqlLrI/AX4X1Wh7CPh6RJyxP4VExKKIeAVwKHDjFH16IqJ3fAGW7s82pUkNDcGOHbBxI2zbBhs2VM+HhuquTNrLbI9CGZdQzYkD50XE24DPA++d7RtFxDOoAnsJ1ej77Mz8/hTdLwIeM/cutdVxx8Ftt8H4D+Nnngm33gpLltRaljTRrI5CeXSliF8CR2fmvS1tLwU+BjwuMxfN4r0OAZ5CdRbny4A/Bk6fLMSbR8H0tDQtBe7wKBRJB4rZHIUy1xH4scB9rQ2Z+amI+CHVj5Ezlpk7gZ82n94UEacAfwa8ZpK+Y1THmgMQnuILine truncated
|
||||
"text/plain": [
|
||||
"array([2.00000003, 1.41421353])"
|
||||
"<Figure size 400x400 with 1 Axes>"
|
||||
]
|
||||
},
|
||||
"execution_count": 10,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
"metadata": {
|
||||
"needs_background": "light"
|
||||
},
|
||||
"output_type": "display_data"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"#import numpy as np\n",
|
||||
"import jax.numpy as np\n",
|
||||
"historyPGz = onp.asarray(historyPGz)\n",
|
||||
"\n",
|
||||
"axes = plt.figure(figsize=(4, 4), dpi=100).gca()\n",
|
||||
"axes.set_title('z space')\n",
|
||||
"axes.scatter(historyPGz[:,0], historyPGz[:,1], lw=0.5, color='red', marker='*')\n",
|
||||
"axes.scatter([0], [0], lw=0.25, color='black', marker='*') \n",
|
||||
"axes.set_xlabel('z0')\n",
|
||||
"axes.set_ylabel('z1')"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"For PG we're making explicit steps in $\\mathbf{z}$ space, which progress in a straight diagonal line to the origin. (Note that in $\\mathbf{z}$ space the origin is likewise the solution.)\n",
|
||||
"\n",
|
||||
"Interestingly, neither GD nor Newton's method give us information about progress in intermediate spaces like $\\mathbf{z}$: for GD we're concatenating the Jacobians, so we're moving in directions that locally should decrease the loss. However, the $\\mathbf{z}$ position is influenced by $\\mathbf{x}$, and hence we don't know where we end up in $\\mathbf{z}$ space until we have the definite point in $\\mathbf{x}$ space.\n",
|
||||
"\n",
|
||||
"With PGs we do not have this problem, as we directly map points in $\\mathbf{z}$ to $\\mathbf{x}$ via an inverse function. Hence we know eactly where we started in $\\mathbf{z}$ space, as this position is crucial for the inversion."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"historyGDz = onp.asarray(historyGDz)\n",
|
||||
"historyNtz = onp.asarray(historyNtz)\n",
|
||||
"\n",
|
||||
"axes = plt.figure(figsize=(4, 4), dpi=100).gca()\n",
|
||||
"axes.set_title('z space')\n",
|
||||
"axes.scatter(historyGDz[:,0], historyGDz[:,1], lw=0.5, marker='*', color='blue')\n",
|
||||
"axes.scatter(historyNtz[:,0], historyNtz[:,1], lw=0.5, marker='*', color='orange')\n",
|
||||
"axes.scatter(historyPGz[:,0], historyPGz[:,1], lw=0.5, marker='*', color='red')\n",
|
||||
"axes.scatter([0], [0], lw=0.25, color='black', marker='*') \n",
|
||||
"axes.set_xlabel('z0')\n",
|
||||
"axes.set_ylabel('z1')"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Conclusions \n",
|
||||
"\n",
|
||||
"That concludes our simple example. And despite its simplicity, it already surprisingly large differences between gradient descent, Newton's method, and the physical gradients emerged.\n",
|
||||
"\n",
|
||||
"The main takeaways were:\n",
|
||||
"* GD easily yields \"unbalanced\" updates\n",
|
||||
"* Newtons method does better, but \n",
|
||||
"* PGs outperform both if an inverse function is available\n",
|
||||
"* Be aware of how an optimizer progresses in latent spaces\n",
|
||||
" \n",
|
||||
"In the next sections we can build on these observations to use PGs for training NNs via invertible physical models."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"\n",
|
||||
"---\n",
|
||||
"\n",
|
||||
"## Approximate Inversions\n",
|
||||
"\n",
|
||||
"If an analytic inverse like the `fun_z_inv_analytic` above is not readily available, we can actually resort to optimization schemes like Newton's method or BFGS to approximate it numerically. This is a topic that is orthogonal to the comparison of different optimization methods, but it can be easily illustrated based on the PG example above.\n",
|
||||
"\n",
|
||||
"Below, we'll use the BFGS variant `fmin_l_bfgs_b` from `scipy` to compute the inverse. It's not very complicated, but we'll use numpy and scipy directly here, which makes the code a bit messier than it should be."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"def fun_z_inv_opt(target_y, x_ini):\n",
|
||||
" # a bit ugly, we switch to pure scipy here inside each iteration for BFGS\n",
|
||||
" import numpy as np\n",
|
||||
@@ -542,7 +725,7 @@
|
||||
" return ret\n",
|
||||
" \n",
|
||||
" ret = fmin_l_bfgs_b(lambda x: fun_z_opt(x,target_y), x_ini, approx_grad=True )\n",
|
||||
" #print( ret ) # full BFGS details\n",
|
||||
" #print( ret ) # return full BFGS details\n",
|
||||
" return ret[0]\n",
|
||||
"\n",
|
||||
"print(\"BFGS optimization test run, find x such that y=[2,2]:\")\n",
|
||||
@@ -550,27 +733,17 @@
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 11,
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"PG iter 0: [2.09999967 2.50998022]\n",
|
||||
"PG iter 1: [1.46999859 2.10000011]\n",
|
||||
"PG iter 2: [1.02899871 1.75698602]\n",
|
||||
"PG iter 3: [0.72029824 1.4699998 ]\n",
|
||||
"PG iter 4: [0.50420733 1.22988982]\n",
|
||||
"PG iter 5: [0.35294448 1.02899957]\n",
|
||||
"PG iter 6: [0.24705997 0.86092355]\n",
|
||||
"PG iter 7: [0.17294205 0.72030026]\n",
|
||||
"PG iter 8: [0.12106103 0.60264817]\n",
|
||||
"PG iter 9: [0.08474171 0.50421247]\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"Nonetheless, we can now use this numerically inverted $\\mathbf{z}$ function to perform the PG optimization. Apart from calling `fun_z_inv_opt`, the rest of the code is unchanged."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"x = np.asarray([3.,3.])\n",
|
||||
"eta = 0.3\n",
|
||||
@@ -598,47 +771,31 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 12,
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"Diff analytic PG versus approximated: 0.000022\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"historyPGa = onp.asarray(history)\n",
|
||||
"updatesPGa = onp.asarray(updates) \n",
|
||||
"\n",
|
||||
"print(\"Diff analytic PG versus approximated: %f\" % (np.sum(np.abs(historyPGa-historyPG))) )\n"
|
||||
"print(\"MAE difference between analytic PG and approximate inversion: %f\" % (np.average(np.abs(historyPGa-historyPG))) )\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Conclusions \n",
|
||||
"That concludes simple example...\n",
|
||||
"<br>\n",
|
||||
"\n",
|
||||
"---\n",
|
||||
"\n",
|
||||
"Main takeaways:\n",
|
||||
"* PGs work neatly\n",
|
||||
"* much faster convergence\n",
|
||||
" \n",
|
||||
"Coming up, complex example."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Next Steps\n",
|
||||
"\n",
|
||||
"Try\n",
|
||||
"- X\n",
|
||||
"- Y"
|
||||
"Based on this code example you can try the following modifications:\n",
|
||||
"\n",
|
||||
"- instead of the simple L(z(x)) contrsuction above, try other, more complicated functions\n",
|
||||
"\n",
|
||||
"- instead of the simple \"regular\" gradient descent, compare the versions above to commonly used DL optimizers such as AdaGrad, RmsProp or Adam."
|
||||
]
|
||||
},
|
||||
{
|
||||
|
||||
+1
-1
@@ -1,4 +1,4 @@
|
||||
Physical Gradients
|
||||
Physical Gradients and NNs
|
||||
=======================
|
||||
|
||||
Re-cap?
|
||||
|
||||
+4
-2
@@ -340,10 +340,12 @@ Even when the Jacobian is singular (because the function is not injective, chaot
|
||||
The update obtained with a regular gradient descent method has surprising shortcomings.
|
||||
The physical gradient instead allows us to more accurately backpropagate through nonlinear functions, provided that we have access to good inverse functions.
|
||||
|
||||
Before moving on to including PGs in NN training processes, the next example will illustrate ...
|
||||
Before moving on to including PGs in NN training processes, the next example will illustrate the differences between these approaches with a practical example.
|
||||
|
||||
|
||||
**todo, integrate comments below?**
|
||||
|
||||
|
||||
**TODO, sometime, integrate comments below?**
|
||||
|
||||
Old Note:
|
||||
The inverse function to a simulator is typically the time-reversed physical process.
|
||||
|
||||
@@ -269,7 +269,7 @@
|
||||
" if not transposed:\n",
|
||||
" block.add_module('%s_conv' % name, nn.Conv2d(in_c, out_c, kernel_size=size, stride=2, padding=pad, bias=True))\n",
|
||||
" else:\n",
|
||||
" block.add_module('%s_upsam' % name, nn.Upsample(scale_factor=2))\n",
|
||||
" block.add_module('%s_upsam' % name, nn.Upsample(scale_factor=2, mode='bilinear'))\n",
|
||||
" # reduce kernel size by one for the upsampling (ie decoder part)\n",
|
||||
" block.add_module('%s_tconv' % name, nn.Conv2d(in_c, out_c, kernel_size=(size-1), stride=1, padding=pad, bias=True))\n",
|
||||
"\n",
|
||||
|
||||
Reference in new issue
Block a user