daviddoji/pbdl-book
1
0 Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

renamed z to y, sensitivity section updates

Browse Source
This commit is contained in:
NT committed 2022-02-05 16:53:49 +01:00
1 parent 571ec6cbb8
commit 3c5bcd3dda
3 files changed
+261 -247
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

physgrad-comparison.ipynb
+134 -127
View File
@@ -11,18 +11,18 @@
"\n", "\n",
"## Problem formulation\n", "## Problem formulation\n",
"\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", "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{y}$, and a scalar loss $L$, i.e. \n",
"$\\mathbf{x} \\in \\mathbb{R}^2$, \n", "$\\mathbf{x} \\in \\mathbb{R}^2$, \n",
"$\\mathbf{z}: \\mathbb{R}^2 \\rightarrow \\mathbb{R}^2 $, and \n", "$\\mathbf{y}: \\mathbb{R}^2 \\rightarrow \\mathbb{R}^2 $, and \n",
"$L: \\mathbb{R}^2 \\rightarrow \\mathbb{R} $.\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", "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", "\n",
"Specifically, we'll use the following $\\mathbf{z}$ and $L$:\n", "Specifically, we'll use the following $\\mathbf{y}$ and $L$:\n",
"\n", "\n",
"$\\quad \\mathbf{z}(\\mathbf{x}) = \\mathbf{z}(x_0,x_1) = \\begin{bmatrix} x_0 \\\\ x_1^2 \\end{bmatrix}$, \n", "$\\quad \\mathbf{y}(\\mathbf{x}) = \\mathbf{y}(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", "i.e. $\\mathbf{y}$ only squares the second component of its input, and\n",
"\n", "\n",
"$\\quad L(\\mathbf{z}) = |\\mathbf{z}|^2 = z_0^2 + z_1^2 \\ $ \n", "$\\quad L(\\mathbf{y}) = |\\mathbf{y}|^2 = y_0^2 + y_1^2 \\ $ \n",
"represents a simple squared $L^2$ loss.\n", "represents a simple squared $L^2$ loss.\n",
"\n", "\n",
"As starting point for some example optimizations we'll use \n", "As starting point for some example optimizations we'll use \n",
@@ -32,7 +32,7 @@
"\n", "\n",
"$\\quad \\text{arg min}_{\\mathbf{x}} \\ L(\\mathbf{x}).$\n", "$\\quad \\text{arg min}_{\\mathbf{x}} \\ L(\\mathbf{x}).$\n",
"\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", "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{y}$ is a very simple function, it is nonlinear due to its $x_1^2$.\n",
"\n", "\n",
"\n" "\n"
] ]
@@ -45,10 +45,10 @@
"## 3 Spaces\n", "## 3 Spaces\n",
"\n", "\n",
"In order to understand the following examples, it's important to keep in mind that we're dealing with mappings between the three _spaces_ we've introduced here:\n", "In order to understand the following examples, it's important to keep in mind that we're dealing 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}$, $\\mathbf{y}$ and $L$. A regular forward pass maps an\n",
"$\\mathbf{x}$ to $L$, while for the optimization we'll need to associate values\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", "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", "be interesting how this influences the positions in $\\mathbf{y}$ that develop while searching for\n",
"the right position in $\\mathbf{x}$.\n", "the right position in $\\mathbf{x}$.\n",
"\n", "\n",
"```{figure} resources/placeholder.png\n", "```{figure} resources/placeholder.png\n",
@@ -74,7 +74,7 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 16, "execution_count": 1,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
@@ -87,21 +87,28 @@
"cell_type": "markdown", "cell_type": "markdown",
"metadata": {}, "metadata": {},
"source": [ "source": [
"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." "We'll start by defining the $\\mathbf{y}$ and $L$ functions, together with a single composite function `fun` which calls L and y. Having a single native python function is necessary for many of the JAX operations."
] ]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 17, "execution_count": 2,
"metadata": {}, "metadata": {},
"outputs": [ "outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"WARNING:absl:No GPU/TPU found, falling back to CPU. (Set TF_CPP_MIN_LOG_LEVEL=0 and rerun for more info.)\n"
]
},
{ {
"name": "stdout", "name": "stdout",
"output_type": "stream", "output_type": "stream",
"text": [ "text": [
"Starting point x = [3. 3.]\n", "Starting point x = [3. 3.]\n",
"\n", "\n",
"Some test calls of the functions we defined so far, from top to bottom, z, Lz (manual), Lz:\n" "Some test calls of the functions we defined so far, from top to bottom, y, manual L(y), L(y):\n"
] ]
}, },
{ {
@@ -112,56 +119,56 @@
" DeviceArray(90., dtype=float32))" " DeviceArray(90., dtype=float32))"
] ]
}, },
"execution_count": 17, "execution_count": 2,
"metadata": {}, "metadata": {},
"output_type": "execute_result" "output_type": "execute_result"
} }
], ],
"source": [ "source": [
"# \"physics\" function z\n", "# \"physics\" function y\n",
"def fun_z(x):\n", "def fun_y(x):\n",
" return np.array( [x[0], x[1]*x[1]] )\n", " return np.array( [x[0], x[1]*x[1]] )\n",
"\n", "\n",
"# simple L2 loss\n", "# simple L2 loss\n",
"def fun_L(z):\n", "def fun_L(y):\n",
" #return z[0]*z[0] + z[1]*z[1] # \"manual version\"\n", " #return y[0]*y[0] + y[1]*y[1] # \"manual version\"\n",
" return np.sum( np.square(z) )\n", " return np.sum( np.square(y) )\n",
"\n", "\n",
"# composite function with L & z\n", "# composite function with L & y\n",
"def fun(x):\n", "def fun(x):\n",
" return fun_L(fun_z(x))\n", " return fun_L(fun_y(x))\n",
"\n", "\n",
"\n", "\n",
"x = np.asarray([3,3], dtype=np.float32)\n", "x = np.asarray([3,3], dtype=np.float32)\n",
"print(\"Starting point x = \"+format(x) +\"\\n\")\n", "print(\"Starting point x = \"+format(x) +\"\\n\")\n",
"\n", "\n",
"print(\"Some test calls of the functions we defined so far, from top to bottom, z, Lz (manual), Lz:\") \n", "print(\"Some test calls of the functions we defined so far, from top to bottom, y, manual L(y), L(y):\") \n",
"fun_z(x) , fun_L( fun_z(x) ), fun(x) " "fun_y(x) , fun_L( fun_y(x) ), fun(x) "
] ]
}, },
{ {
"cell_type": "markdown", "cell_type": "markdown",
"metadata": {}, "metadata": {},
"source": [ "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 \\mathbf{z}$. The cell below evaluates this and a few variants, together with a sanity check for the inverse of the Jacobian of $\\mathbf{z}$:" "Now we can evaluate the derivatives of our function via `jax.grad`. E.g., `jax.grad(fun_L)(fun_y(x))` evaluates the Jacobian $\\partial L / \\partial \\mathbf{y}$. The cell below evaluates this and a few variants, together with a sanity check for the inverse of the Jacobian of $\\mathbf{y}$:"
] ]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 2, "execution_count": 3,
"metadata": {}, "metadata": {},
"outputs": [ "outputs": [
{ {
"name": "stdout", "name": "stdout",
"output_type": "stream", "output_type": "stream",
"text": [ "text": [
"Jacobian L(z): [ 6. 18.]\n", "Jacobian L(y): [ 6. 18.]\n",
"\n", "\n",
"Jacobian z(x): \n", "Jacobian y(x): \n",
"[[1. 0.]\n", "[[1. 0.]\n",
" [0. 6.]]\n", " [0. 6.]]\n",
"\n", "\n",
"Sanity check with inverse Jacobian of z, this should give x again: [3. 3.]\n", "Sanity check with inverse Jacobian of y, this should give x again: [3. 3.]\n",
"\n", "\n",
"Gradient for full L(x): [ 6. 108.]\n", "Gradient for full L(x): [ 6. 108.]\n",
"\n" "\n"
@@ -170,19 +177,19 @@
], ],
"source": [ "source": [
"# this works:\n", "# this works:\n",
"print(\"Jacobian L(z): \" + format(jax.grad(fun_L)(fun_z(x))) +\"\\n\")\n", "print(\"Jacobian L(y): \" + format(jax.grad(fun_L)(fun_y(x))) +\"\\n\")\n",
"\n", "\n",
"# the following would give an error as z (and hence fun_z) is not scalar\n", "# the following would give an error as y (and hence fun_y) is not scalar\n",
"#jax.grad(fun_z)(x) \n", "#jax.grad(fun_y)(x) \n",
"\n", "\n",
"# computing the jacobian of z is a valid operation:\n", "# computing the jacobian of y is a valid operation:\n",
"J = jax.jacobian(fun_z)(x)\n", "J = jax.jacobian(fun_y)(x)\n",
"print( \"Jacobian z(x): \\n\" + format(J) ) \n", "print( \"Jacobian y(x): \\n\" + format(J) ) \n",
"\n", "\n",
"# the following also gives error, JAX grad needs a single function object\n", "# the following also gives error, JAX grad needs a single function object\n",
"#jax.grad( fun_L(fun_z) )(x) \n", "#jax.grad( fun_L(fun_y) )(x) \n",
"\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", "print( \"\\nSanity check with inverse Jacobian of y, this should give x again: \" + format(np.linalg.solve(J, np.matmul(J,x) )) +\"\\n\")\n",
"\n", "\n",
"# instead use composite 'fun' from above\n", "# instead use composite 'fun' from above\n",
"print(\"Gradient for full L(x): \" + format( jax.grad(fun)(x) ) +\"\\n\")\n" "print(\"Gradient for full L(x): \" + format( jax.grad(fun)(x) ) +\"\\n\")\n"
@@ -215,9 +222,9 @@
"$$\\begin{aligned}\n", "$$\\begin{aligned}\n",
"\\Delta \\mathbf{x} \n", "\\Delta \\mathbf{x} \n",
"&= \n", "&= \n",
"- \\eta ( J_{L} J_{\\mathbf{z}} )^T \\\\\n", "- \\eta ( J_{L} J_{\\mathbf{y}} )^T \\\\\n",
"&=\n", "&=\n",
"- \\eta ( \\frac{\\partial L }{ \\partial \\mathbf{z} } \\frac{\\partial \\mathbf{z} }{ \\partial \\mathbf{x} } )^T\n", "- \\eta ( \\frac{\\partial L }{ \\partial \\mathbf{y} } \\frac{\\partial \\mathbf{y} }{ \\partial \\mathbf{x} } )^T\n",
"\\end{aligned}$$\n", "\\end{aligned}$$\n",
"\n", "\n",
"where $\\eta$ denotes the step size parameter .\n", "where $\\eta$ denotes the step size parameter .\n",
@@ -228,7 +235,7 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 3, "execution_count": 4,
"metadata": {}, "metadata": {},
"outputs": [ "outputs": [
{ {
@@ -271,7 +278,7 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 4, "execution_count": 5,
"metadata": {}, "metadata": {},
"outputs": [ "outputs": [
{ {
@@ -280,7 +287,7 @@
"Text(0, 0.5, 'x1')" "Text(0, 0.5, 'x1')"
] ]
}, },
"execution_count": 4, "execution_count": 5,
"metadata": {}, "metadata": {},
"output_type": "execute_result" "output_type": "execute_result"
}, },
@@ -330,7 +337,7 @@
" \\frac{\\partial L }{ \\partial \\mathbf{x} }\n", " \\frac{\\partial L }{ \\partial \\mathbf{x} }\n",
"\\\\\n", "\\\\\n",
"&=\n", "&=\n",
"- \\eta \\ H_L^{-1} \\ ( J_{L} J_{\\mathbf{z}} )^T\n", "- \\eta \\ H_L^{-1} \\ ( J_{L} J_{\\mathbf{y}} )^T\n",
"\\end{aligned}\n", "\\end{aligned}\n",
"$$\n", "$$\n",
"\n", "\n",
@@ -345,7 +352,7 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 5, "execution_count": 6,
"metadata": {}, "metadata": {},
"outputs": [ "outputs": [
{ {
@@ -373,7 +380,7 @@
"for i in range(10):\n", "for i in range(10):\n",
" G = jax.grad(fun)(x)\n", " G = jax.grad(fun)(x)\n",
" H = jax.jacobian(jax.jacobian(fun))(x)\n", " H = jax.jacobian(jax.jacobian(fun))(x)\n",
" #H = jax.jacfwd(jax.jacrev(fun_Lz))(x) # alternative\n", " #H = jax.jacfwd(jax.jacrev(fun_Ly))(x) # alternative\n",
" Hinv = np.linalg.inv(H)\n", " Hinv = np.linalg.inv(H)\n",
" \n", " \n",
" x += -eta * np.matmul( Hinv , G)\n", " x += -eta * np.matmul( Hinv , G)\n",
@@ -393,7 +400,7 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 6, "execution_count": 7,
"metadata": {}, "metadata": {},
"outputs": [ "outputs": [
{ {
@@ -402,7 +409,7 @@
"Text(0, 0.5, 'x1')" "Text(0, 0.5, 'x1')"
] ]
}, },
"execution_count": 6, "execution_count": 7,
"metadata": {}, "metadata": {},
"output_type": "execute_result" "output_type": "execute_result"
}, },
@@ -443,23 +450,23 @@
"source": [ "source": [
"## Physical Gradients\n", "## Physical Gradients\n",
"\n", "\n",
"Now we also use inverse physics, i.e. the inverse of z:\n", "Now we also use inverse physics, i.e. the inverse of y:\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. This gives an update step:\n", "$\\mathbf{y}^{-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. This gives an update step:\n",
"\n", "\n",
"$$\\begin{aligned}\n", "$$\\begin{aligned}\n",
"\\Delta \\mathbf{x} &= \n", "\\Delta \\mathbf{x} &= \n",
"\\mathbf{z}^{-1} \\left( \\mathbf{z}(\\mathbf{x}) - \\eta\n", "\\mathbf{y}^{-1} \\left( \\mathbf{y}(\\mathbf{x}) - \\eta\n",
" \\left( \\frac{\\partial^2 L }{ \\partial \\mathbf{z}^2 } \\right)^{-1}\n", " \\left( \\frac{\\partial^2 L }{ \\partial \\mathbf{y}^2 } \\right)^{-1}\n",
" \\frac{\\partial L }{ \\partial \\mathbf{z} }\n", " \\frac{\\partial L }{ \\partial \\mathbf{y} }\n",
"\\right) - \\mathbf{x}\n", "\\right) - \\mathbf{x}\n",
"\\end{aligned}$$\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" "Below, we define our inverse function `fun_y_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", "cell_type": "code",
"execution_count": 7, "execution_count": 8,
"metadata": {}, "metadata": {},
"outputs": [ "outputs": [
{ {
@@ -482,25 +489,25 @@
"source": [ "source": [
"x = np.asarray([3.,3.])\n", "x = np.asarray([3.,3.])\n",
"eta = 0.3\n", "eta = 0.3\n",
"historyPG = [x]; historyPGz = []; updatesPG = []\n", "historyPG = [x]; historyPGy = []; updatesPG = []\n",
"\n", "\n",
"def fun_z_inv_analytic(y):\n", "def fun_y_inv_analytic(y):\n",
" return np.array( [y[0], np.power(y[1],0.5)] )\n", " return np.array( [y[0], np.power(y[1],0.5)] )\n",
"\n", "\n",
"for i in range(10):\n", "for i in range(10):\n",
" \n", " \n",
" # Newton step for L(z)\n", " # Newton step for L(y)\n",
" zForw = fun_z(x)\n", " zForw = fun_y(x)\n",
" GL = jax.grad(fun_L)(zForw)\n", " GL = jax.grad(fun_L)(zForw)\n",
" HL = jax.jacobian(jax.jacobian(fun_L))(zForw)\n", " HL = jax.jacobian(jax.jacobian(fun_L))(zForw)\n",
" HLinv = np.linalg.inv(HL)\n", " HLinv = np.linalg.inv(HL)\n",
" \n", " \n",
" # step in z space\n", " # step in y space\n",
" zBack = zForw -eta * np.matmul( HLinv , GL)\n", " zBack = zForw -eta * np.matmul( HLinv , GL)\n",
" historyPGz.append(zBack)\n", " historyPGy.append(zBack)\n",
"\n", "\n",
" # \"inverse physics\" step via z-inverse\n", " # \"inverse physics\" step via y-inverse\n",
" x = fun_z_inv_analytic(zBack)\n", " x = fun_y_inv_analytic(zBack)\n",
" historyPG.append(x)\n", " historyPG.append(x)\n",
" updatesPG.append( historyPG[-2] - historyPG[-1] )\n", " updatesPG.append( historyPG[-2] - historyPG[-1] )\n",
" print( \"PG iter %d: \"%i + format(x) )\n", " print( \"PG iter %d: \"%i + format(x) )\n",
@@ -518,7 +525,7 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 8, "execution_count": 9,
"metadata": {}, "metadata": {},
"outputs": [ "outputs": [
{ {
@@ -527,7 +534,7 @@
"Text(0, 0.5, 'x1')" "Text(0, 0.5, 'x1')"
] ]
}, },
"execution_count": 8, "execution_count": 9,
"metadata": {}, "metadata": {},
"output_type": "execute_result" "output_type": "execute_result"
}, },
@@ -561,14 +568,14 @@
"cell_type": "markdown", "cell_type": "markdown",
"metadata": {}, "metadata": {},
"source": [ "source": [
"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", "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{y}$, 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", "\n",
"This difference also shows in first update step for each method: below we measure how well it is aligned with the 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", "cell_type": "code",
"execution_count": 9, "execution_count": 10,
"metadata": {}, "metadata": {},
"outputs": [ "outputs": [
{ {
@@ -606,18 +613,18 @@
"\n", "\n",
"<br>\n", "<br>\n",
"\n", "\n",
"## Z Space\n", "## y Space\n",
"\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 in between 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. We will return to NNs soon, but for now let's focus on $\\mathbf{z}$. \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 in between that matter. In our case, we only have a single $\\mathbf{y}$ space, but for DL settings, we might have a large number of latent spaces, over which we have a certain amount of control. We will return to NNs soon, but for now let's focus on $\\mathbf{y}$. \n",
"\n", "\n",
"A 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 coincidentally, we already stored it in the `historyPGz` list above.\n", "A first thing to note is that for PG, we explicitly map from $L$ to $\\mathbf{y}$, and then continue with a mapping to $\\mathbf{x}$. Thus we already obtained the trajectory in $\\mathbf{y}$ space, and not coincidentally, we already stored it in the `historyPGy` list above.\n",
"\n", "\n",
"Let's directly take a look what PG did in $\\mathbf{z}$ space:" "Let's directly take a look what PG did in $\\mathbf{y}$ space:"
] ]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 10, "execution_count": 11,
"metadata": {}, "metadata": {},
"outputs": [ "outputs": [
{ {
@@ -626,13 +633,13 @@
"Text(0, 0.5, 'z1')" "Text(0, 0.5, 'z1')"
] ]
}, },
"execution_count": 10, "execution_count": 11,
"metadata": {}, "metadata": {},
"output_type": "execute_result" "output_type": "execute_result"
}, },
{ {
"data": { "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 "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAGDCAYAAAAyH0uIAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAc+klEQVR4nO3de5RdZZnn8e9DAhWXpArtQuju6KCFTcbbIATBTA9gFNBup9XltVHHmL4oLa4WZ8YWWwGdQYYZF6HpMC0daZHuRFg9beOIIzReclOCNN4K76QlQnMtIaeIQBHDM3/sU+RUUVWpqpxT+7yp72etvU6dd7/77OfsVH7ZefctMhNJUnkOqLsASdLsGOCSVCgDXJIKZYBLUqEMcEkqlAEuSYUywCWpUAa4JBXKAJekQhngklQoA1ySCmWAS1KhDHBJKpQBrmJExMsiIiPidRPMO70576VTLH9gRJwbET+NiEcj4hcRsSUiTmnpc0VE7IyI50TE9RHxy4i4KyLOiYgY93n/JSK+0fycRyLiloh4wyTrfltEfDMiHo6IByNiU0ScOq7PqyJic3OdD0XEFyPi+TPfUpovDHCVZANwB/DWCea9FdiWmTdOsfx5wLnA14AzgfOBnwPHjOu3ALgOuBf4AHAL8NHm1OpPgW8D5wAfAn4F/H1E/G5rp4g4F/hbYFez77nN77Gipc/bgS8CO4E/A/4b8DxgS0QcMcV30nyWmU5OxUzAx4FHgb6WtkOpwvG8vSz7HeDavfS5Akjgkpa2AK4FRoD+lvanjFv2QGAQ+EpL25HAbuBzwAHj+kfz9WDgQeCvx80/DNgxvt3JaXRyD1yluRLoAVqHKt4MLAT+bi/L7gCeHxHPncZ61oz+kJnZfH8Q8IqW9kdGf46IpwF9wGbG7tG/lup/uh/LzMdbV9D8XIBTgEOAz0ZE/+hEFfw3AS+bRr2ahxbWXYA0E5n5o4i4mWrI5PJm81uBrZl5214WPwf4PPCTiLiVapjkbzPze+P6PQ78y7i2nzRfjxhtiIhXAx8Gjqb6R+WJMlt+Hmh+3g+mqGv0H5SvTjJ/eIplNY8Z4CrRlcBfRMQSquA8gWpMe0qZuSkiBoDXAKcCfwicFRHvzsxPzaSAiPgPwP8FNgF/AtxNNYzzTuD0mXwWe45FvR24Z4L5v5rh52meMMBVoquAi4DfB55CFZxXT2fBzHwA+DTw6Yg4mCqAzwNaA/wA4Dns2esG+K3m6+3N19dTjcWflpkjo50i4p3jVrmt+XnPoxqDn8i25ut9mfnl6XwPCTwLRQXKzCHgS8DbqIZPrmu2TSkifm3c5+wEbmPs8MeoJ/bom6cPnkn1D8VXms27qYZKFrT0O4JqzLvVNVRDKOdExJi/by2nJV5PNUzyoYg4cIK6D53yi2necg9cpboS+D/Nnz8yzWV+EBEbqE4LfABYRnUwdM24fo8Cr4yIz1AdRHwV8LvAxzPz/mafLwLvB66LiPXAM4D3UP2D8KLRD8rM2yLi/GaNmyPic1RnsxwH3AWcnZnDEXEG1amG34qIq4D7gWc11/t1pjFEpPkn9hwIl8oREQdRjRcfAByemY9OY5k/B36PajikB9hOFZr/KzN3NftcQRXqLwL+Cvht4CHgk4w7kyQiVgEfpAranwEXUh3kPDczx1/0807gvVRDKQ8D3wP+e+uQSUSc3Py8E5r1/SvVWS1rMvOWaW4azSMGuIoUEQup9mC/kJl/0MbPvQJ4Q2Ye3K7PlDrFMXCV6rVUF/BcWXMdUm0cA1dRIuJ4quGNjwDfzsyNNZck1cY9cJXmDKqx6fuA/1RzLVKtHAOXpEK5By5JhTLAJalQRR/EbF7J9htU5+lK0v5iMXBX7mWMu+gApwrvO+suQpI6YAnVxVyTKj3AHwK444476O3trbsWSdpnw8PDPPOZz4RpjCyUHuAA9Pb2GuCS5h0PYkpSoQxwSSqUAS5JhTLAJalQBrgkFcoAl6RO6fC9pgxwSeqE3bth2bLqtUMMcEnqhM2bYXAQtmzp2Cr2iwt5JKkrDA1Bo1H9vHo17NpVvS5ZUrX19UF/f9tW5x64JLXL+efD0qWwfDls2lS1bdxYvV+6tJrfRga4JLXL6tVw2WWwYAHs2FG17dgBCxfC2rXV/DYywCWpnVatghUrxratWAErV7Z9VY6BS1I7ZVbDJv39MDAA27bBhg1Ve0RbV+UeuCS10+Ag7NwJ69fD1q2wbl31fnCw7asq+qHGEdELNBqNhreTldQdRkaqqTWTGg1YtAh6eva6+PDwMH19fQB9mTk8VV+HUCSpnXp6nhzUVSC3nUMoklQoA1ySCmWAS1KhDHBJKlTtAR4RvxkRfxcRv4iIRyJiMCKW1V2XJHW7Ws9CiYinAV8Hvga8CrgfeC7wYJ11SVIJ6j6N8M+AOzLznS1tP6urGEkqSd1DKL8H/HNE/H1E3BcR346IP5qsc0T0RETv6AQsnrtSJam71B3gzwHOAH4KnAb8FXBJRLxjkv5nA42W6c65KFKSulGtl9JHxGPAP2fm8pa2S4DjMvOlE/TvAVovcVoM3Oml9JL2FyVdSn838INxbT8EXj9R58wcAUZG30eb7+wlSSWpewjl68BR49p+C9heQy2SVJS6A3w1cEJEfCgijoyI04E/Bi6tuS5J6nq1Bnhm3gy8Dvh94FbgI8D7MnNdnXVJUgnqHgMnM68Frq27DkkqTd1DKJKkWTLAJalQBrgkFcoAl6RCGeCSVCgDXJIKZYBLUqEMcEkqlAEuSYUywCWpUAa4JBXKAJekQhngklQoA1ySCmWAS1KhDHBJKpQBLkmFMsAlqVAGuCQVygCXpEIZ4JJUKANckgplgEtSoQxwSSqUAS5JhTLAJalQBrgkFcoAl6RCGeCSVCgDXJIKZYBLUqEMcEkqVK0BHhHnRUSOm35UZ02SVIqFdRcAfB94Rcv7X9VViCSVpBsC/FeZeU/dRUhSabphDPy5EXFXRPxLRKyLiGdN1jEieiKid3QCFs9hnZLUVeoO8JuAlcArgTOAZwObI2KyYD4baLRMd85BjZLUlSIz667hCRFxCLAdeH9mXj7B/B6gp6VpMXBno9Ggt7d3boqUpA4aHh6mr68PoC8zh6fq2w1j4E/IzB0R8RPgyEnmjwAjo+8jYq5Kk6SuU/cQyhgRcTAwANxddy2S1O3qPg/8ExFxUkQcERHLgX8EdgOfrbMuSSpB3UMoS6jC+teA+4EtwAmZeX+tVUlSAWoN8Mx8S53rl6SSddUYuCRp+gxwSSqUAS5JhTLAJalQBrgkFcoAl6RCGeCSVCgDXJIKZYBLUqEMcEkqlAEuSYUywCWpUAa4JBXKAJekQhngkvboomfkau8McEmV3bth2bLqVUUwwCVVNm+GwUHYsqXuSjRNdT9STVKdhoag0ah+Xr0adu2qXpcsqdr6+qC/v776NCX3wKX57PzzYelSWL4cNm2q2jZurN4vXVrNV9cywKX5bPVquOwyWLAAduyo2nbsgIULYe3aar66lgEuzXerVsGKFWPbVqyAlStrKUfT5xi4NN9lVsMm/f0wMADbtsGGDVV7RN3VaQrugUvz3eAg7NwJ69fD1q2wbl31fnCw7sq0F5EFn7gfEb1Ao9Fo0NvbW3c5UplGRqqp9e9QowGLFkFPT311zVPDw8P09fUB9GXm8FR9HUKR5ruenicHdRUg6nIOoUhSoQxwSSqUAS5JhTLAJalQBrgkFcoAl6RCGeCSVCgDXJIK1TUBHhEfjIiMiIvrrkWSStAVAR4RxwHvAr5Xdy2SVIraAzwiDgbWAX8EPFhzOZJUjNoDHLgU+GJmfnlvHSOiJyJ6RydgcefLk6TuVOvNrCLiLcAxwHHTXORs4NzOVSRJ5ahtDzwingn8BfDWzHx0motdAPS1TEs6VJ4kdb0698CPBZ4BfCv2PPVjAXBiRJwJ9GTm7tYFMnMEGBl9Hz4tRNI8VmeAfwV44bi2TwM/Ai4cH96SpLFqC/DMfAi4tbUtIn4J/CIzb514KUnSqG44C0WSNAtd9Ui1zDy57hokqRTugUtSoQxwSSqUAS5JhTLAJalQBrgkFcoAl6RCGeCSVCgDXJIKZYBLUqEMcEkqlAEuSYUywCWpUAa4JBXKAJekQhngklQoA1ySCmWAS1KhDHBJKpQBLkmFMsAlqVAGuCQVygCXpEIZ4JJUKANckgplgEtSodoa4BHx7yJidzs/U5I0sU7sgUcHPlOSNM7CmXSOiM/tpUsfkLMvR+oymRDuk6g7zXQP/D8Ci4DGJNPOtlYn1Wn3bli2rHqVutCM9sCBHwL/kJmXTzQzIo4GXr2vRUldYfNmGByELVvgpJPqrkZ6kpkG+C3AMcCEAQ6MAD/fp4qkOg0NQaNR/bx6NezaVb0uWVK19fVBf3999UktInP6Q9YR0QMsyMyHO1fS9EVEL9BoNBr09vbWXY72B2edBWvWwNOfDo89Bjt2wCGHwEEHwQMPwJlnVoEudcjw8DB9fX0AfZk5PFXfGY2BZ+ZIZj4cEV+NiHPHz4+Ip0XEV6f7eRFxRkR8LyKGm9ONEfGqmdQktdXq1XDZZbBgQRXeUL0uXAhr1xre6iqzPY3wZODMiLgmIp7a0n4QMJPBwjuBDwLHAsuArwKfj4jnz7Iuad+tWgUrVoxtW7ECVq6spRxpMvtyHvgrgMOBrRFxxGw+IDO/kJn/LzN/mpk/ycw/pzqT5YR9qEvaN5mwcWM11n388dXrhg1Vu9RF9iXA76ba2x4Ebo6Ik/elkIhYEBFvAZ4K3DhJn56I6B2dgMX7sk5pQoODsHMnrF8PW7fCunXV+8HBuiuTxpjpWSijEqoxceD0iPgwcB1w4Uw/KCJeSBXYi6j2vl+XmT+YpPvZwJPG3qW2Ouoo2L4dRg+Mn3oq3H47LFpUa1nSeDM6C+WJhSIeBw7PzPta2l4PfAZ4SmYumMFnHQQ8i+oqzjcAfwicNFGIN8+C6WlpWgzc6VkLine truncated
"text/plain": [ "text/plain": [
"<Figure size 400x400 with 1 Axes>" "<Figure size 400x400 with 1 Axes>"
] ]
@@ -644,11 +651,11 @@
} }
], ],
"source": [ "source": [
"historyPGz = onp.asarray(historyPGz)\n", "historyPGy = onp.asarray(historyPGy)\n",
"\n", "\n",
"axes = plt.figure(figsize=(4, 4), dpi=100).gca()\n", "axes = plt.figure(figsize=(4, 4), dpi=100).gca()\n",
"axes.set_title('z space')\n", "axes.set_title('y space')\n",
"axes.scatter(historyPGz[:,0], historyPGz[:,1], lw=0.5, color='red', marker='*')\n", "axes.scatter(historyPGy[:,0], historyPGy[:,1], lw=0.5, color='red', marker='*')\n",
"axes.scatter([0], [0], lw=0.25, color='black', marker='*') \n", "axes.scatter([0], [0], lw=0.25, color='black', marker='*') \n",
"axes.set_xlabel('z0')\n", "axes.set_xlabel('z0')\n",
"axes.set_ylabel('z1')" "axes.set_ylabel('z1')"
@@ -658,7 +665,7 @@
"cell_type": "markdown", "cell_type": "markdown",
"metadata": {}, "metadata": {},
"source": [ "source": [
"For PG we're making explicit steps in $\\mathbf{z}$ space, which progress in a straight diagonal line to the origin (which is likewise the solution in $\\mathbf{z}$ space).\n", "For PG we're making explicit steps in $\\mathbf{y}$ space, which progress in a straight diagonal line to the origin (which is likewise the solution in $\\mathbf{y}$ space).\n",
"\n" "\n"
] ]
}, },
@@ -666,51 +673,51 @@
"cell_type": "markdown", "cell_type": "markdown",
"metadata": {}, "metadata": {},
"source": [ "source": [
"Interestingly, neither GD nor Newton's method give us information about progress in intermediate spaces (like the $\\mathbf{z}$ space). \n", "Interestingly, neither GD nor Newton's method give us information about progress in intermediate spaces (like the $\\mathbf{y}$ space). \n",
"\n", "\n",
"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 there. (For NNs in general we won't know at which latent-space points we end up after a GD update until we've actually computed all updated weights.)\n", "For **GD** we're concatenating the Jacobians, so we're moving in directions that locally should decrease the loss. However, the $\\mathbf{y}$ position is influenced by $\\mathbf{x}$, and hence we don't know where we end up in $\\mathbf{y}$ space until we have the definite point there. (For NNs in general we won't know at which latent-space points we end up after a GD update until we've actually computed all updated weights.)\n",
"\n", "\n",
"More specifically, we have an update $-\\eta \\frac{\\partial L}{\\partial \\mathbf{x}}$ for GD, which means we arrive at $\\mathbf{z}(\\mathbf{x} -\\eta \\frac{\\partial L}{\\partial \\mathbf{x}})$ in $\\mathbf{z}$ space. A Taylor expansion with \n", "More specifically, we have an update $-\\eta \\frac{\\partial L}{\\partial \\mathbf{x}}$ for GD, which means we arrive at $\\mathbf{y}(\\mathbf{x} -\\eta \\frac{\\partial L}{\\partial \\mathbf{x}})$ in $\\mathbf{y}$ space. A Taylor expansion with \n",
"$h = \\eta \\frac{\\partial L}{\\partial \\mathbf{x}}$ yields \n", "$h = \\eta \\frac{\\partial L}{\\partial \\mathbf{x}}$ yields \n",
"\n", "\n",
"$\n", "$\n",
"\\quad\n", "\\quad\n",
"\\mathbf{z}(\\mathbf{x} - h) = \n", "\\mathbf{y}(\\mathbf{x} - h) = \n",
"\\mathbf{z}(\\mathbf{x}) - h \\frac{\\partial \\mathbf{z}}{\\partial \\mathbf{x}} + \\mathcal{O}( h^2 )\n", "\\mathbf{y}(\\mathbf{x}) - h \\frac{\\partial \\mathbf{y}}{\\partial \\mathbf{x}} + \\mathcal{O}( h^2 )\n",
"= \\mathbf{z}(x) - \\eta \\frac{\\partial L}{\\partial \\mathbf{z}} (\\frac{\\partial \\mathbf{z}}{\\partial x})^2 + \\mathcal{O}( h^2 )\n", "= \\mathbf{y}(x) - \\eta \\frac{\\partial L}{\\partial \\mathbf{y}} (\\frac{\\partial \\mathbf{y}}{\\partial x})^2 + \\mathcal{O}( h^2 )\n",
"$.\n", "$.\n",
"\n", "\n",
"And $\\frac{\\partial L}{\\partial \\mathbf{z}} (\\frac{\\partial \\mathbf{z}}{\\partial \\mathbf{x}})^2$ clearly differs from the step $\\frac{\\partial L}{\\partial \\mathbf{z}}$ we would compute during the backpropagation pass in GD for $\\mathbf{z}$.\n", "And $\\frac{\\partial L}{\\partial \\mathbf{y}} (\\frac{\\partial \\mathbf{y}}{\\partial \\mathbf{x}})^2$ clearly differs from the step $\\frac{\\partial L}{\\partial \\mathbf{y}}$ we would compute during the backpropagation pass in GD for $\\mathbf{y}$.\n",
"\n", "\n",
"**Newton's method** does not fare much better: we compute first-order derivatives like for GD, and the second-order derivatives for the Hessian for the full process. But since both are approximations, the actual intermediate states resulting from an update step are unknown until the full chain is evaluated. In the _Consistency in function compositions_ paragraph for Newton's method in {doc}`physgrad` the squared $\\frac{\\partial \\mathbf{z}}{\\partial \\mathbf{x}}$ term for the Hessian already indicated this dependency.\n", "**Newton's method** does not fare much better: we compute first-order derivatives like for GD, and the second-order derivatives for the Hessian for the full process. But since both are approximations, the actual intermediate states resulting from an update step are unknown until the full chain is evaluated. In the _Consistency in function compositions_ paragraph for Newton's method in {doc}`physgrad` the squared $\\frac{\\partial \\mathbf{y}}{\\partial \\mathbf{x}}$ term for the Hessian already indicated this dependency.\n",
"\n", "\n",
"With **PGs** we do not have this problem: PGs can directly map points in $\\mathbf{z}$ to $\\mathbf{x}$ via the inverse function. Hence we know eactly where we started in $\\mathbf{z}$ space, as this position is crucial for evaluating the inverse.\n", "With **PGs** we do not have this problem: PGs can directly map points in $\\mathbf{y}$ to $\\mathbf{x}$ via the inverse function. Hence we know eactly where we started in $\\mathbf{y}$ space, as this position is crucial for evaluating the inverse.\n",
"\n", "\n",
"In the simple setting of this section, we only have a single latent space, and we already stored all values in $\\mathbf{x}$ space during the optimization (in the `history` lists). Hence, now we can go back and re-evaluate `fun_z` to obtain the positions in $\\mathbf{z}$ space." "In the simple setting of this section, we only have a single latent space, and we already stored all values in $\\mathbf{x}$ space during the optimization (in the `history` lists). Hence, now we can go back and re-evaluate `fun_y` to obtain the positions in $\\mathbf{y}$ space."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"x = np.asarray([3.,3.])\n",
"eta = 0.01\n",
"historyGDz = []\n",
"historyNtz = []\n",
"\n",
"for i in range(1,10):\n",
" historyGDz.append(fun_z(historyGD[i]))\n",
" historyNtz.append(fun_z(historyNt[i]))\n",
"\n",
"historyGDz = onp.asarray(historyGDz)\n",
"historyNtz = onp.asarray(historyNtz)\n"
] ]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 12, "execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"x = np.asarray([3.,3.])\n",
"eta = 0.01\n",
"historyGDy = []\n",
"historyNty = []\n",
"\n",
"for i in range(1,10):\n",
" historyGDy.append(fun_y(historyGD[i]))\n",
" historyNty.append(fun_y(historyNt[i]))\n",
"\n",
"historyGDy = onp.asarray(historyGDy)\n",
"historyNty = onp.asarray(historyNty)\n"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": { "metadata": {
"scrolled": true "scrolled": true
}, },
@@ -721,13 +728,13 @@
"Text(0, 0.5, 'z1')" "Text(0, 0.5, 'z1')"
] ]
}, },
"execution_count": 12, "execution_count": 13,
"metadata": {}, "metadata": {},
"output_type": "execute_result" "output_type": "execute_result"
}, },
{ {
"data": { "data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAGDCAYAAAA26pu1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAlkElEQVR4nO3df5zcVX3v8dcnm2SjkF2FhQQbKbhC0luwYH6Ae70EIpDa670Kiliu95rGtmrF22oftdLij/aay/W2dSkNahpBbZuo9NrSVisGCwmJEGjRyoISdGuAFEiyxp1lhewuu5/7x5kh311mktndmTnfM/N+Ph7fx+R75nx3PplMPnvmfM8Pc3dERCRdc2IHICIis6NELiKSOCVyEZHEKZGLiCROiVxEJHFK5CIiiVMiFxFJnBK5iEjilMhFRBKnRC4ikjglchGRxCmRi4gkTolcRCRxSuSSFDM7zcy80nGMaxea2fVmttfMRszsgJndbmavztTZbmYPmtlyM7vbzJ41sx+Z2bun/Kz5ZvaHZna/mRXM7KdmttPMLirzunPM7DfNrM/MDpvZQTO7zcxWTKn39uLPe9bMDpnZl8zs5bN9z6T5mZaxlZSY2XHAZVOK5wG9wKi7n3yUa7cAbwE2At8DTgReC3zZ3bcU62wHzgDmArcAjwBvLdZ7p7vfXKzXBTwAfBH4AbAQeCfwCmCVu/9r5nU/B6wDvg58o/iz/xPwTXffWKzz+8D/Kr7mDuAk4H3AMHCuuw9O422SVuPuOnQkfQA3As8BFx2j3iCw8Rh1tgMOfCBTNh/4DrAfmFcsawPmT7n2JcBTwE2ZsouKP+9Py7xWqSH1s8X4f2/K82cBY1PLdeiYeqhrRZJmZv8D+A3gg+5+5zGqDwLnmdnLjlHvOWBT6cTdR4vnJwPLi2XjxfJS18kJhJb2vwCvzvysNxMS+R9MfRF3L30dvpzQzXmLmXWVDsIvhR8QfhmIVKRELskys3OAzwBfdPdPVnHJBwmt3MfN7D4z+5iZvaJMvSfc/adTyh4pPp6Wef13mNkDwGHgx8BB4D8DnZnruos/79BR4joDMELSPjjl+DnCLxCRiubGDkBkJszspcBXCAn2V6u5xt1vMbOdhD72S4HfAX7XzC53969P8/XfDnweuBX4I+AAMA5cQ0je0zGH0Gp/ffFnTDU8zZ8nLUaJXJJjZnOALYQ+6Yvd/Zlqr3X3J4FPAZ8ys5OBbwO/T7gRWfIyMztuSqv8zOLj3uLjW4B/Ay7PdJFgZlO7UPqBtWZ2wlFa5f2EFvmP3P2RCnVEKlLXiqToo8Ba4Jfd/UfVXGBmbWaW7fLA3Q8ATwDtU6rPBd6VuXZ+8fwgcH+xuNRytky984DXTPlZXynW+WiZmErX/k3x5300U/Z8HTM78Vh/P2ltapFLUszsbODDwF3AycUujue5+19VuHQhsM/M/h/wXUJ3xcXASuC3p9R9gtDlchqh6+ZK4Bzg1919rFjnq4SblH9rZl8DTgfeTRjWeHwmnjvN7C+B/2lmZwC3ERpQ/wm4kzCKpt/MrgWuA04zs1uBp4s/8zLgz4E/rvItkhakceSSFDO7kJAAy3J3K1debFV/nNA3/gpCMv0hsMndP52ptx3oAt4B/BlwLmHY4R+5+42ZegZ8iNBSX0xI4B8GrgAudPfTMnXbgPdzZJx5gTC65Vp3/3am3uXFeucWix4H/gm4QV0ucjRK5CIZpUTu7mfFjkWkWuojFxFJnBK5iEjilMhFRBKnPnIRkcSpRS4ikjglchGRxCU9Iag4lvdlhMkTIiLNYiFhsbWq+r6TTuSEJL4vdhAiInWwBPj3aiqmnsifBnj88cfp6OiIHYuIyKwNDQ3x8pe/HKbR05B6Igego6NDiVxEWpZudoqIJE6JXEQkcUrkIiKJUyIXEUmcErmISOKUyEVakdZYaipK5CKtZmIcblsRHqUpKJGLtJqDO6HQBwd3xY5EaqQpJgSJyDEcHoCxQvjzw70wMQZ7euHFS0LZvE5Y0BUvPpkVJXKRVvDQBnhkI8w/ASZGQ9n+HbCtB0YPwZlXw/LeuDHKjKlrRaQVLO+FVZvA2mBsMJSNDYLNhVWblcQTFzWRm9leM/Myx40x4xJpSt3rYdGayWWL10D3uijhSO3E7lpZCbRlzs8Cbgf+Ok44Ik3MHQ7sgPYuOL4bhvth//ZQbhY7OpmFqC1ydz/o7k+VDuANQD+wI2ZcIk1psA+eG4aerbB2N/RsCeeDfbEjk1nKzebLZjYfeAL4pLv/7wp12oH2TNFCYF+hUNAytiLHMj4CEyMwL/N/ZbQAbQugrb3yddJQQ0NDdHZ2AnS6+1A11+TpZuebgJcAnz9KnWuAQubQ7kAi1Wprn5zEAeZ3Kok3gTwl8ncCX3f3J45S5zqgM3MsaURgIiJ5FvtmJwBm9rPAxcDlR6vn7iPASOa6OkcmIpJ/eWmR/wpwAPha7EBERFITPZGb2RxCIv+Cuz8XOx4RkdRET+SELpVTgZtjByIikqLofeTuvg1QZ7eIyAzloUUuIiKzoEQuIpI4JXIRkcQpkYuIJE6JXEQkcUrkIiKJUyIXEUmcErmISOKUyEVEEqdELiKSOCVyEZHEKZGLiCROiVxEJHFK5CIiiVMiFxFJnBK5iFTmHjsCqYISuYiUNz4OK1aER8k1JXIRKW/nTujrg127YkcixxB9qzcRyZGBASgUwp97e2FsLDwuWRLKOjuhqytefFKWWuQicsSGDbBsGfT0wF13hbIdO8L5smXheckdJXIROaK3FzZtgrY2GBwMZYODMHcubN4cnpfcUSIXkcnWr4c1ayaXrVkD69ZFCUeOTX3kIjKZe+hO6eqC7m7o74ft20O5WezopAy1yEVksr4+GB6GrVth927YsiWc9/XFjkwqME94wL+ZdQCFQqFAR0dH7HBEmsPISDiy/6cKBViwANrb48XVIoaGhujs7ATodPehaq5R14qITNbe/sKEHRKL5JS6VkREEhc9kZvZz5jZX5nZj83sWTPrM7MVseMSEUlF1K4VM3sp8C3gTuD1wEHgDOAnMeMSEUlJ7D7y3wUed/dfyZT9KFYwIiIpit218l+BfzGzvzazA2b2HTP7tUqVzazdzDpKB7CwcaGKiORT7ET+CuA9wA+AtcCngRvM7B0V6l8DFDLHvkYEKSKSZ1HHkZvZKPAv7t6TKbsBWOnurylTvx3IjotaCOzTOHIRaRYpjiN/EvjelLLvA28uV9ndR4CR0rlpurCISPSulW8BS6eUnQk8GiEWEZEkxU7kvcD5ZvZ7ZvZKM7sK+HXgxshxiYgkI2oid/d/Bi4Dfhl4EPgw8FvuviVmXCIiKYndR467fxX4auw4RERSFbtrRUREZkmJXEQkcUrkIiKJUyKX1pTwhioiUymRS+uZGIfbVoRHkSagRC6t5+BOKPTBwV2xIxGpiejDD0Ua4vAAjBXCnx/uhYkx2NMLL14SyuZ1woKuePGJzIISubSGhzbAIxth/gkwMRrK9u+AbT0wegjOvBqW98aNUWSG1LUirWF5L6zaBNYGY4OhbGwQbC6s2qwkLklTIpfW0b0eFq2ZXLZ4DXSvixKOSK2oa0Vahzsc2AHtXXB8Nwz3w/7toVxLIkvC1CKX1jHYB88NQ89WWLsberaE88G+2JGJzErUHYJmq7hvZ0E7BElVxkdgYgTmZT4rowVoWwBt7ZWvE2mgFHcIEmmctvYXJuz5nXFiEakhda2IiCROiVxEJHFK5CIiiVMiFxFJnBK5iEjilMhFRBKnRC4ikjglchGRxCmRi4gkTolcRCRxSuQiIolTIhcRSZwSuYhI4pTIRUQSFzWRm9nHzMynHA/HjElEJDV5WI/8IeDizPlzsQIREUlRHhL5c+7+VOwgRERSlYc+8jPM7Akz+zcz22Jmp1aqaGbtZtZROoCFDYxTRCSXYifye4F1wC8C7wFOB3aaWaUEfQ1QyBz7GhCjiEiu5WrzZTN7CfAo8AF3v6nM8+1AdtPFhcA+bb4sIs0i+c2X3X3QzB4BXlnh+RFgpHRuZo0KTUQkt2J3rUxiZscD3cCTsWMREUlF7HHkf2xmq83sNDPrAf4WGAe+GDMuEZGUxO5aWUJI2icCB4FdwPnufjBqVCIiCYmayN39bTFfX0SkGeSqj1xERKZPiVxEJHFK5CIiiVMiFxFJnBK5iEjilMhFRBKnRC4ikjglchGRxCmRi4gkTolcRCRxSuRSOzla216klSiRS21MjMNtK8Jjs9IvKskpJXKpjYM7odAHB3fFjqQ+xsdhxYrwKJIzSuQyc4cH4On+cDzcCxNjsKf3SNnhgdgR1s7OndDXB7ua9BeVJC32euSSsoc2wCMbYf4JMDEayvbvgG09MHoIzrwalvfGjXE2BgagUAh/7u2FsbHwuGRJKOvshK6uePGJFKlFLjO3vBdWbQJrg7HBUDY2CDYXVm1OO4kDbNgAy5ZBTw/cdVco27EjnC9bFp4XyQElcpmd7vWwaM3kssVroHtdlHBqqrcXNm2CtjYYHAxlg4Mwdy5s3hyeF8kBJXKZHXc4sAPau+DE88Lj/u3NM8Jj/XpYM+UX1Zo1sG5dlHBEylEil9kZ7IPnhqFnK6zdDT1bwvlgX+zIasM9dKd0dcF554XH7dub5xeVNAUlcpmdjqXwpkfhlEvC+SmXwhv3hvJm0NcHw8OwdSvs3g1btoTzvib5RSVNwTzhloWZdQCFQqFAR0dH7HCkGY2MhCP7+SoUYMECaG+PF5c0raGhITo7OwE63X2omms0/FDkaNrbX5iww38ykdxQ14qISOKUyEVEEqdELiKSOCVyEZHEKZGLiCROiVxEJHG5SeRm9iEzczO7PnYsIiIpyUUiN7OVwLuAB2LHIiKSmuiJ3MyOB7YAvwb8JHI4IiLJiZ7IgRuBr7n7N49V0czazayjdAAL6x+eiEi+RZ2ib2ZvA14NrKzykmuAj9YvIhGR9ERrkZvZy4ELine truncated "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAGDCAYAAAA26pu1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAlxElEQVR4nO3dfZycZX3v8c8vm2SDJruKiwk2UnCV5FSlah7APR6CkZDaegr4hEV7GmNbS42vVs9pKy0i9phyPOfUpTSoaQQRm/h02toWKwTFPAEBS21ZUELZmkAKJFnjzmYL2Sy7v/PHNUPuXWaS2d2Zve5r5vt+ve7X7H3Nfe/928nkN9dc9/Vg7o6IiKRrRuwARERkapTIRUQSp0QuIpI4JXIRkcQpkYuIJE6JXEQkcUrkIiKJUyIXEUmcErmISOKUyEVEEqdELiKSOCVyEZHEKZGLiCROiVySYWZvNjM3s0vLPHd58bk3nuD8WWb2CTP7VzM7amY/MbNdZrYqc8zNZjZoZq8ws9vN7D/M7Akzu9rMbNzv+x9mdnfx9zxjZveb2TsrXPt9ZnafmT1tZj81sx1mdtG4Y95qZjuL1zxiZt8ys1dP/JWSZqNELinZBjwOvLfMc+8Fet39nhOcfw3wCeB7wDpgPfAY8IZxx7UAtwEHgN8H7gc+Wdyyfgf4AXA18IfAs8A3zOyXsgeZ2SeALwPDxWM/Ufw7VmaO+VXgW8Ag8AfA/wR+DthlZmee4G8SAXfXpi2ZDfgT4CjQnik7jZAkrznJuf8M3HqSY24GHLg+U2bArcAQ0JEpP2XcubOAHuC7mbJXAiPAXwMzxh1vxce5wE+Bvxj3/Hygf3y5Nm3jN9XIJTW3AK1AtgnjMmAm8JcnObcfeLWZvaqK62wo/eDuXtyfDVyYKX+m9LOZvRhoB3YytoZ/CeGb7x+7+2j2AsXfC7AKeBHwFTPrKG2ED4B7gTdXEa80sZmxAxCZCHd/2My+T2hKubFY/F5gt7s/epLTrwb+FnjEzB4kNJ982d0fGHfcKPBv48oeKT6eWSows7cBVwGvI3y4PBdm5ufO4u/74QniKn2w3Fnh+YETnCuiRC5JugX4MzNbSEig5xHavE/I3XeYWSdwMXAR8OvAR8zst9z9CxMJwMz+C/B3wA7gt4EnCc077wcun8jv4vi9ql8Fnirz/LMT/H3SZJTIJUVfBT4D/ApwCiGBfq2aE939MPBF4ItmNpeQiK8Bsol8BvAKjtfCAc4uPu4tPr6D0Fa/2t2HSgeZ2fvHXbK3+Pt+jtBGX05v8fGgu3+nmr9DJEtt5JIcd+8Dvg28j9Csclux7ITM7CXjfs8g8Chjm0VKnqvhF7sdriN8YHy3WDxCaEJpyRx3JqFNPOubhKaVq81szP+3THfG2wnNJ39oZrPKxH3aCf8waXqqkUuqbgH+X/Hnj1d5zg/NbBuhO+FhYCnhpumGcccdBX7BzL5EuNn4VuCXgD9x90PFY74FfBS4zcy2AC8FPkT4YDin9Ivc/VEzW1+McaeZ/TWh98sy4AngSncfMLMrCF0U/8nMvgocAs4oXvcuqmg6kuZlx2+ci6TDzGYT2pNnAAvc/WgV5/wR8MuEZpJWYB8hef4fdx8uHnMzIbmfA3wOeBNwBPg843qemNla4GOEhPtj4NOEm6GfcPfxg4feD3yY0MTyNPAA8KlsU4qZXVD8fecV4/t3Qi+YDe5+f5UvjTQhJXJJkpnNJNRo/97dP1DD33sz8E53n1ur3ylSb2ojl1RdQhgIdEvkOESiUxu5JMXMziU0e3wc+IG7b48ckkh0qpFLaq4gtF0fBP5b5FhEckFt5CIiiVONXEQkcUrkIiKJS/pmZ3Fk3MsI/XxFRBrFPOAJr7LtO+lETkji+2MHISJSBwsJg8JOKvVEfgTg8ccfp62tLXYsIiJTNjAwwMtf/nKYQEtD6okcgLa2NiVyEWlautkpIpI4JXIRkcQpkYuIJE6JXEQkcUrkIiKJUyIXaUaaY6mhKJGLNJvREbhtaXiUhqBELtJsDu2EQg8c2hU7EqmRhhgQJCIncbQPhgvh54e7YXQY9nTDCxaGslntMKcjXnwyJUrkIs3gofXwyAaYfSqMHgtlB7bD1i44dhjOXgdLuuPGKJOmphWRZrCkG5ZvBGuB4f5QNtwPNhOWb1IST1zURG5me83My2w3xIxLpCF1roX5K8eWLVgJnWuihCO1E7tpZRnQktl/DXAH8I044Yg0MHc4uB1aO2BuJwz2woFtodwsdnQyBVFr5O5+yN2fKm3A24BeQCuji9Rafw88OwhdW2D1bujaHPb7e2JHJlOUm8WXzWw28ATwGXf/kwrHtAKtmaJ5wP5CoaBpbEVOZmQIRodgVub/yrECtMyBltbK58m0GhgYoL29HaDd3QeqOSdPNzsvAV4E3HyCY64ECplNqwOJVKuldWwSB5jdriTeAPKUyD8AfNvdnzjBMdcC7Zlt4XQEJiKSZ7FvdgJgZj8LXAi8/UTHufsQMJQ5r86RiYjkX15q5O8HDgLfih2IiEhqoidyM5tBSORfcvdnY8cjIpKa6Imc0KRyBnBT7EBERFIUvY3c3bcCauwWEZmkPNTIRURkCpTIRUQSp0QuIpI4JXIRkcQpkYuIJE6JXEQkcUrkIiKJUyIXEUmcErmISOKUyEVEEqdELiKSOCVyEZHEKZGLiCROiVxEJHFK5CIiiVMiF5HK3GNHIFVQIheR8kZGYOnS8Ci5pkQuIuXt3Ak9PbBrV+xI5CSiL/UmIjnS1weFQvi5uxuGh8PjwoWhrL0dOjrixSdlqUYuIsetXw+LF0NXF+zYEcq2bw/7ixeH5yV3lMhF5Ljubti4EVpaoL8/lPX3w8yZsGlTeF5yR4lcRMZauxZWrhxbtnIlrFkTJRw5ObWRi8hY7qE5paMDOjuhtxe2bQvlZrGjkzJUIxeRsXp6YHAQtmyB3bth8+aw39MTOzKpwDzhDv9m1gYUCoUCbW1tscMRaQxDQ2HL/p8qFGDOHGhtjRdXkxgYGKC9vR2g3d0HqjlHTSsiMlZr6/MTdkgsklNqWhERSVz0RG5mP2Nmf2lmPzGzZ8ysx8yWxo5LRCQVUZtWzOzFwF3A94C3AoeAVwE/jRmXiEhKYreR/wHwuLu/P1P241jBiIikKHbTyi8D/2hm3zCzg2b2AzP7jUoHm1mrmbWVNmDe9IUqIpJPsRP5K4ArgH8FVgOfA643s1+rcPyVQCGz7Z+OIEVE8ixqP3IzOwb8o7t3ZcquB5a5+xvLHN8KZPtFzQP2qx+5iDSKFPuRPwn8cFzZj4B3lDvY3YeAodK+abiwiEj0ppW7gEXjys4G9kWIRUQkSbETeTdwnpn9oZm90swuB34TuCFyXCIiyYiayN39+8ClwK8ADwIfB37X3TfHjEtEJCWx28hx91uBW2PHISKSqthNKyIiMkVK5CIiiVMiFxFJnBK5NKeEF1QRGU+JXJrP6AjctjQ8ijQAJXJpPod2QqEHDu2KHYlITUTvfigyLY72wXAh/PxwN4wOw55ueMHCUDarHeZ0xItPZAqUyKU5PLQeHtkAs0+F0WOh7MB22NoFxw7D2etgSXfcGEUmSU0r0hyWdMPyjWAtMNwfyob7wWbC8k1K4pI0JXJpHp1rYf7KsWULVkLnmijhiNSKmlakebjDwe3Q2gFzO2GwFw5sC+WaElkSphq5NI/+Hnh2ELq2wOrd0LU57Pf3xI5MZEqirhA0VcV1OwtaIUiqMjIEo0MwK/NeOVaAljnQ0lr5PJFplOIKQSLTp6X1+Ql7dnucWERqSE0rIiKJUyIXEUmcErmISOKUyEVEEqdELiKSOCVyEZHEKZGLiCROiVxEJHFK5CIiiVMiFxFJnBK5iEjilMhFRBKnRC4ikjglchGRxEVN5GZ2jZn5uO3hmDGJiKQmD/ORPwRcmNl/NlYgIiIpykMif9bdn4odhIhIqvLQRv4qM3vCzP7NzDab2RmVDjSzVjNrK23AvGmMU0Qkl2In8nuBNcAvAFcAZwE7zaxSgr4SKGS2/dMQo4hIruVq8WUzexGwD/iou99Y5vlWILvo4jxgvxZfFpFGkfziy+7eb2aPAK+s8PwQMFTaN7PpCk1EJLdiN62MYWZzgU7gydixiIikInY/8v9rZivM7Ewz6wL+BhgBvhIzLhGRlMRuWllISNovAQ4Bu4Dz3P1Q1KhERBISNZG7+3tiXl9EpBHkqo1cREQmTolcRCRxSuQiIolTIhcRSZwSuYhI4pTIRUQSp0QuIpI4JXIRkcQpkYuIJE6JXEQkcUrkUjs5mttepJkokUttjI7AbUvDY6PSB5XklBK51MahnVDogUO7YkdSHyMjsHRpeBTJGSVymbyjfXCkN2wPd8PoMOzpPl52tC92hLWzcyf09MCuBv2gkqTFno9cUvbQenhkA8w+FUaPhbID22FrFxw7DGevgyXdcWOcir4+KBTCz93dMDwcHhcuDGXt7dDRES8+kSLVyGXylnTD8o1gLTDcH8qG+8FmwvJNaSdxgPXrYfFi6OqCHTtC2fbtYX/x4vC8SA4okcvUdK6F+SvHli1YCZ1rooRTU93dsHEjtLRAf38o6++HmTNh06bwvEgOKJHL1LjDwe3Q2gEvOTc8HtjWOD081q6FleM+qFauhDVrooQjUo4SuUxNfw88OwhdW2D1bujaHPb7e2JHVhvuoTmlowPOPTc8btvWOB9U0hCUyGVq2hbBJfvg9FVh//SL4OK9obwR9PTA4CBs2QK7d8PmzWG/p0E+qKQhmCdcszCzNqBQKBRoa2uLHY40oqGhsGXfX4UCzJkDra3x4pKGNTAwQHt7O0C7uw9Uc466H4qcSGvr8xN2+E8mkhtqWhERSZwSuYhI4pTIRUQSp0QuIpI4JXIRkcQpkYuIJC43idzMPmZmbmbXxY5FRCQluUjkZrYM+CDwQOxYRERSEz2Rm9lcYDPwG8BPI4cjIpKc6IkcuAHLine truncated
"text/plain": [ "text/plain": [
"<Figure size 400x400 with 1 Axes>" "<Figure size 400x400 with 1 Axes>"
] ]
@@ -740,10 +747,10 @@
], ],
"source": [ "source": [
"axes = plt.figure(figsize=(4, 4), dpi=100).gca()\n", "axes = plt.figure(figsize=(4, 4), dpi=100).gca()\n",
"axes.set_title('z space')\n", "axes.set_title('y space')\n",
"axes.scatter(historyGDz[:,0], historyGDz[:,1], lw=0.5, marker='*', color='blue')\n", "axes.scatter(historyGDy[:,0], historyGDy[:,1], lw=0.5, marker='*', color='blue')\n",
"axes.scatter(historyNtz[:,0], historyNtz[:,1], lw=0.5, marker='*', color='orange')\n", "axes.scatter(historyNty[:,0], historyNty[:,1], lw=0.5, marker='*', color='orange')\n",
"axes.scatter(historyPGz[:,0], historyPGz[:,1], lw=0.5, marker='*', color='red')\n", "axes.scatter(historyPGy[:,0], historyPGy[:,1], lw=0.5, marker='*', color='red')\n",
"axes.scatter([0], [0], lw=0.25, color='black', marker='*') \n", "axes.scatter([0], [0], lw=0.25, color='black', marker='*') \n",
"axes.set_xlabel('z0')\n", "axes.set_xlabel('z0')\n",
"axes.set_ylabel('z1')" "axes.set_ylabel('z1')"
@@ -753,7 +760,7 @@
"cell_type": "markdown", "cell_type": "markdown",
"metadata": {}, "metadata": {},
"source": [ "source": [
"These trajectories confirm the intuition outlined in the previous sections: GD in blue gives a very sub-optimal trajectory in $\\mathbf{z}$. Newton (in orange) does better, but is still clearly curved, in contrast to the straight, and diagonal red trajectory for the PG-based optimization.\n", "These trajectories confirm the intuition outlined in the previous sections: GD in blue gives a very sub-optimal trajectory in $\\mathbf{y}$. Newton (in orange) does better, but is still clearly curved, in contrast to the straight, and diagonal red trajectory for the PG-based optimization.\n",
"\n", "\n",
"The behavior in intermediate spaces becomes especially important when they're not only abstract latent spaces as in this example, but when they have actual physical meanings." "The behavior in intermediate spaces becomes especially important when they're not only abstract latent spaces as in this example, but when they have actual physical meanings."
] ]
@@ -784,14 +791,14 @@
"\n", "\n",
"## Approximate inversions\n", "## Approximate inversions\n",
"\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", "If an analytic inverse like the `fun_y_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", "\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." "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", "cell_type": "code",
"execution_count": 13, "execution_count": 14,
"metadata": {}, "metadata": {},
"outputs": [ "outputs": [
{ {
@@ -807,42 +814,42 @@
"array([2.00000003, 1.41421353])" "array([2.00000003, 1.41421353])"
] ]
}, },
"execution_count": 13, "execution_count": 14,
"metadata": {}, "metadata": {},
"output_type": "execute_result" "output_type": "execute_result"
} }
], ],
"source": [ "source": [
"def fun_z_inv_opt(target_y, x_ini):\n", "def fun_y_inv_opt(target_y, x_ini):\n",
" # a bit ugly, we switch to pure scipy here inside each iteration for BFGS\n", " # a bit ugly, we switch to pure scipy here inside each iteration for BFGS\n",
" import numpy as np\n", " import numpy as np\n",
" from scipy.optimize import fmin_l_bfgs_b\n", " from scipy.optimize import fmin_l_bfgs_b\n",
" target_y = onp.array(target_y)\n", " target_y = onp.array(target_y)\n",
" x_ini = onp.array(x_ini)\n", " x_ini = onp.array(x_ini)\n",
"\n", "\n",
" def fun_z_opt(x,target_y=[2,2]):\n", " def fun_y_opt(x,target_y=[2,2]):\n",
" y = onp.array( [x[0], x[1]*x[1]] ) # we cant use fun_z from JAX here\n", " y = onp.array( [x[0], x[1]*x[1]] ) # we cant use fun_y from JAX here\n",
" ret = onp.sum( onp.square(y-target_y) )\n", " ret = onp.sum( onp.square(y-target_y) )\n",
" return ret\n", " return ret\n",
" \n", " \n",
" ret = fmin_l_bfgs_b(lambda x: fun_z_opt(x,target_y), x_ini, approx_grad=True )\n", " ret = fmin_l_bfgs_b(lambda x: fun_y_opt(x,target_y), x_ini, approx_grad=True )\n",
" #print( ret ) # return full BFGS details\n", " #print( ret ) # return full BFGS details\n",
" return ret[0]\n", " return ret[0]\n",
"\n", "\n",
"print(\"BFGS optimization test run, find x such that y=[2,2]:\")\n", "print(\"BFGS optimization test run, find x such that y=[2,2]:\")\n",
"fun_z_inv_opt([2,2], [3,3])\n" "fun_y_inv_opt([2,2], [3,3])\n"
] ]
}, },
{ {
"cell_type": "markdown", "cell_type": "markdown",
"metadata": {}, "metadata": {},
"source": [ "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." "Nonetheless, we can now use this numerically inverted $\\mathbf{y}$ function to perform the PG optimization. Apart from calling `fun_y_inv_opt`, the rest of the code is unchanged."
] ]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 14, "execution_count": 15,
"metadata": {}, "metadata": {},
"outputs": [ "outputs": [
{ {
@@ -869,12 +876,12 @@
"\n", "\n",
"for i in range(10): \n", "for i in range(10): \n",
" # same as before, Newton step for L(y)\n", " # same as before, Newton step for L(y)\n",
" y = fun_z(x)\n", " y = fun_y(x)\n",
" GL = jax.grad(fun_L)(y)\n", " GL = jax.grad(fun_L)(y)\n",
" y += -eta * np.matmul( np.linalg.inv( jax.jacobian(jax.jacobian(fun_L))(y) ) , GL)\n", " y += -eta * np.matmul( np.linalg.inv( jax.jacobian(jax.jacobian(fun_L))(y) ) , GL)\n",
"\n", "\n",
" # optimize for inverse physics, assuming we dont have access to an inverse for fun_z\n", " # optimize for inverse physics, assuming we dont have access to an inverse for fun_y\n",
" x = fun_z_inv_opt(y,x)\n", " x = fun_y_inv_opt(y,x)\n",
" history.append(x)\n", " history.append(x)\n",
" updates.append( history[-2] - history[-1] )\n", " updates.append( history[-2] - history[-1] )\n",
" print( \"PG iter %d: \"%i + format(x) )\n" " print( \"PG iter %d: \"%i + format(x) )\n"
@@ -889,7 +896,7 @@
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 15, "execution_count": 16,
"metadata": {}, "metadata": {},
"outputs": [ "outputs": [
{ {
@@ -919,7 +926,7 @@
"\n", "\n",
"Based on this code example you can try the following modifications:\n", "Based on this code example you can try the following modifications:\n",
"\n", "\n",
"- Instead of the simple L(z(x)) function above, try other, more complicated functions.\n", "- Instead of the simple L(y(x)) function above, try other, more complicated functions.\n",
"\n", "\n",
"- Replace the simple \"regular\" gradient descent with another optimizer, e.g., commonly used DL optimizers such as AdaGrad, RmsProp or Adam. Compare the versions above with the new trajectories." "- Replace the simple \"regular\" gradient descent with another optimizer, e.g., commonly used DL optimizers such as AdaGrad, RmsProp or Adam. Compare the versions above with the new trajectories."
] ]
physgrad-nn.md
+20 -20
View File
@@ -8,7 +8,7 @@ In this section, we will first show how PGs can be integrated into the optimizat
## Physical Gradients and loss functions ## Physical Gradients and loss functions
As before, we consider a scalar objective function $L(z)$ that depends on the result of an invertible simulator $z = \mathcal P(x)$. In {doc}`physgrad` we've outlined the inverse gradient (IG) update $\Delta x = \frac{\partial x}{\partial L} \cdot \Delta L$, where $\Delta L$ denotes a step to take in terms of the loss. As before, we consider a scalar objective function $L(y)$ that depends on the result of an invertible simulator $y = \mathcal P(x)$. In {doc}`physgrad` we've outlined the inverse gradient (IG) update $\Delta x = \frac{\partial x}{\partial L} \cdot \Delta L$, where $\Delta L$ denotes a step to take in terms of the loss.
By applying the chain rule and substituting the IG $\frac{\partial x}{\partial L}$ for the PG, we obtain By applying the chain rule and substituting the IG $\frac{\partial x}{\partial L}$ for the PG, we obtain
@@ -17,27 +17,27 @@ $$
\Delta x \Delta x
&= \frac{\partial x}{\partial L} \cdot \Delta L &= \frac{\partial x}{\partial L} \cdot \Delta L
\\ \\
&= \frac{\partial x}{\partial z} \left( \frac{\partial z}{\partial L} \cdot \Delta L \right) &= \frac{\partial x}{\partial y} \left( \frac{\partial y}{\partial L} \cdot \Delta L \right)
\\ \\
&= \frac{\partial x}{\partial z} \cdot \Delta z &= \frac{\partial x}{\partial y} \cdot \Delta y
\\ \\
&= \mathcal P^{-1}_{(x_0,z_0)}(z_0 + \Delta z) - x_0 + \mathcal O(\Delta z^2) &= \mathcal P^{-1}_{(x_0,y_0)}(y_0 + \Delta y) - x_0 + \mathcal O(\Delta y^2)
. .
\end{aligned} \end{aligned}
$$ $$
This equation has turned the step w.r.t. $L$ into a step in $z$ space: $\Delta z$. This equation has turned the step w.r.t. $L$ into a step in $y$ space: $\Delta y$.
However, it does not prescribe a unique way to compute $\Delta z$ since the derivative $\frac{\partial z}{\partial L}$ as the right-inverse of the row-vector $\frac{\partial L}{\partial z}$ puts almost no restrictions on $\Delta z$. However, it does not prescribe a unique way to compute $\Delta y$ since the derivative $\frac{\partial y}{\partial L}$ as the right-inverse of the row-vector $\frac{\partial L}{\partial y}$ puts almost no restrictions on $\Delta y$.
Instead, we use a Newton step (equation {eq}`quasi-newton-update`) to determine $\Delta z$ where $\eta$ controls the step size of the optimization steps. Instead, we use a Newton step (equation {eq}`quasi-newton-update`) to determine $\Delta y$ where $\eta$ controls the step size of the optimization steps.
Here an obvious questions is: Doesn't this leave us with the disadvantage of having to compute the inverse Hessian, as discussed before? Here an obvious questions is: Doesn't this leave us with the disadvantage of having to compute the inverse Hessian, as discussed before?
Luckily, unlike with regular Newton or quasi-Newton methods, where the Hessian of the full system is required, here, the Hessian is needed only for $L(z)$. Even better, for many typical $L$ its computation can be completely forgone. Luckily, unlike with regular Newton or quasi-Newton methods, where the Hessian of the full system is required, here, the Hessian is needed only for $L(y)$. Even better, for many typical $L$ its computation can be completely forgone.
E.g., consider the case $L(z) = \frac 1 2 || z^\textrm{predicted} - z^\textrm{target}||_2^2$ which is the most common supervised objective function. E.g., consider the case $L(y) = \frac 1 2 || y^\textrm{predicted} - y^\textrm{target}||_2^2$ which is the most common supervised objective function.
Here $\frac{\partial L}{\partial z} = z^\textrm{predicted} - z^\textrm{target}$ and $\frac{\partial^2 L}{\partial z^2} = 1$. Here $\frac{\partial L}{\partial y} = y^\textrm{predicted} - y^\textrm{target}$ and $\frac{\partial^2 L}{\partial y^2} = 1$.
Using equation {eq}`quasi-newton-update`, we get $\Delta z = \eta \cdot (z^\textrm{target} - z^\textrm{predicted})$ which can be computed without evaluating the Hessian. Using equation {eq}`quasi-newton-update`, we get $\Delta y = \eta \cdot (y^\textrm{target} - y^\textrm{predicted})$ which can be computed without evaluating the Hessian.
Once $\Delta z$ is determined, the gradient can be backpropagated to earlier time steps using the inverse simulator $\mathcal P^{-1}$. We've already used this combination of a Newton step for the loss and PGs for the PDE in {doc}`physgrad-comparison`. Once $\Delta y$ is determined, the gradient can be backpropagated to earlier time steps using the inverse simulator $\mathcal P^{-1}$. We've already used this combination of a Newton step for the loss and PGs for the PDE in {doc}`physgrad-comparison`.
## NN training ## NN training
@@ -50,8 +50,8 @@ we aim for using PGs to accurately optimize through the simulation.
Consider the following setup: Consider the following setup:
A neural network makes a prediction $x = \mathrm{NN}(a \,;\, \theta)$ about a physical state based on some input $a$ and the network weights $\theta$. A neural network makes a prediction $x = \mathrm{NN}(a \,;\, \theta)$ about a physical state based on some input $a$ and the network weights $\theta$.
The prediction is passed to a physics simulation that computes a later state $z = \mathcal P(x)$, and hence The prediction is passed to a physics simulation that computes a later state $y = \mathcal P(x)$, and hence
the objective $L(z)$ depends on the result of the simulation. the objective $L(y)$ depends on the result of the simulation.
```{admonition} Combined training algorithm ```{admonition} Combined training algorithm
@@ -59,27 +59,27 @@ the objective $L(z)$ depends on the result of the simulation.
To train the weights $\theta$ of the NN, we then perform the following updates: To train the weights $\theta$ of the NN, we then perform the following updates:
* Evaluate $\Delta z$ via a Newton step as outlined above * Evaluate $\Delta y$ via a Newton step as outlined above
* Compute the PG $\Delta x = \mathcal P^{-1}_{(x, z)}(z + \Delta z) - x$ using an inverse simulator * Compute the PG $\Delta x = \mathcal P^{-1}_{(x, y)}(y + \Delta y) - x$ using an inverse simulator
* Use GD or a GD-based optimizer to compute the updates to the network weights, $\Delta\theta = \eta_\textrm{NN} \cdot \frac{\partial y}{\partial\theta} \cdot \Delta y$ * Use GD or a GD-based optimizer to compute the updates to the network weights, $\Delta\theta = \eta_\textrm{NN} \cdot \frac{\partial y}{\partial\theta} \cdot \Delta y$
``` ```
The combined optimization algorithm depends on both the **learning rate** $\eta_\textrm{NN}$ for the network as well as the step size $\eta$ from above, which factors into $\Delta z$. The combined optimization algorithm depends on both the **learning rate** $\eta_\textrm{NN}$ for the network as well as the step size $\eta$ from above, which factors into $\Delta y$.
To first order, the effective learning rate of the network weights is $\eta_\textrm{eff} = \eta \cdot \eta_\textrm{NN}$. To first order, the effective learning rate of the network weights is $\eta_\textrm{eff} = \eta \cdot \eta_\textrm{NN}$.
We recommend setting $\eta$ as large as the accuracy of the inverse simulator allows, before choosing $\eta_\textrm{NN} = \eta_\textrm{eff} / \eta$ to achieve the target network learning rate. We recommend setting $\eta$ as large as the accuracy of the inverse simulator allows, before choosing $\eta_\textrm{NN} = \eta_\textrm{eff} / \eta$ to achieve the target network learning rate.
This allows for nonlinearities of the simulator to be maximally helpful in adjusting the optimization direction. This allows for nonlinearities of the simulator to be maximally helpful in adjusting the optimization direction.
**Note:** **Note:**
For simple objectives like a loss of the form $L=|z - z^*|^2$, this procedure can be easily integrated into an GD autodiff pipeline by replacing the gradient of the simulator only. For simple objectives like a loss of the form $L=|y - y^*|^2$, this procedure can be easily integrated into an GD autodiff pipeline by replacing the gradient of the simulator only.
This gives an effective objective function for the network This gives an effective objective function for the network
$$ $$
L_\mathrm{NN} = \frac 1 2 | x - \mathcal P_{(x,z)}^{-1}(z + \Delta z) |^2 L_\mathrm{NN} = \frac 1 2 | x - \mathcal P_{(x,y)}^{-1}(y + \Delta y) |^2
$$ $$
where $\mathcal P_{(x,z)}^{-1}(z + \Delta z)$ is treated as a constant. where $\mathcal P_{(x,y)}^{-1}(y + \Delta y)$ is treated as a constant.
## Iterations and time dependence ## Iterations and time dependence
physgrad.md
+107 -100
View File
@@ -10,17 +10,20 @@ The latter two methods share similarities, but in the loss term case, the evalua
This is a natural choice from a deep learning perspective, but we haven't questioned at all whether this is actually a good choice. This is a natural choice from a deep learning perspective, but we haven't questioned at all whether this is actually a good choice.
Not too surprising after this introduction: A central insight of the following chapter will be that regular gradients are often a _sub-optimal choice_ for learning problems involving physical quantities. Not too surprising after this introduction: A central insight of the following chapter will be that regular gradients are often a _sub-optimal choice_ for learning problems involving physical quantities.
It turns out that both supervised and DP gradients have their pros and cons. In the following, we'll analyze this in more detail. In particular, we'll illustrate how the multi-modal problems (as hinted at in {doc}`intro-teaser`) negatively influence NNs. Then we'll show how scaling problems of DP gradients affect NN training. Finally, we'll explain several alternatives to prevent these problems. It turns out that a key property that is missing in regular gradients is a proper _inversion_ of the Jacobian matrix. It turns out that both supervised and DP gradients have their pros and cons. In the following, we'll analyze this in more detail. In particular,
we'll show how scaling problems of DP gradients affect NN training.
Then, we'll also illustrate how the multi-modal problems (as hinted at in {doc}`intro-teaser`) negatively influence NNs.
Finally, we'll explain several alternatives to prevent these problems. It turns out that a key property that is missing in regular gradients is a proper _inversion_ of the Jacobian matrix.
```{admonition} A preview of this chapter ```{admonition} A preview of this chapter
:class: tip :class: tip
Below, we'll proceed in the following steps: Below, we'll proceed in the following steps:
- We'll illustrate how the multi-modal problems (as hinted at in {doc}`intro-teaser`) negatively influence NNs - We'll show how scaling problems of DP gradients affect NNs,
- Then we'll show how scaling problems of DP gradients affect NN training. - and how multi-modal problems (cf. {doc}`intro-teaser`) deteriorate NN training
- Finally we'll explain several alternatives to prevent these problems. - Finally we'll explain several alternatives to prevent these problems.
- It turns out that a key property that is missing in regular gradients is a proper _inversion_ of the Jacobian matrix. - What's missing in our GD/Adam&Co. runs so far is a proper _inversion_ of the Jacobian matrix.
``` ```
@@ -31,70 +34,58 @@ Below, we'll proceed in the following steps:
% [analytical, more practical approach] % [analytical, more practical approach]
XXX PG physgrad chapter notes from dec 23 XXX XXX PG physgrad chapter notes from dec 23 XXX
- recap formulation P(x)=z , L() ... etc. rename z to y? [? recap formulation P(x)=z , L() ?]
- intro after dL/dx bad, Newton? discussion is repetitive - intro after dL/dx bad, Newton? discussion is repetitive
[older commment - more intro to quasi newton?] [older commment - more intro to quasi newton?]
- GD - is "diff. phys." , rename? add supervised before? - GD - is "diff. phys." , rename? add supervised before?
comparison notebook: comparison notebook:
- why z, rename to y? (see above)
- add legends to plot - add legends to plot
- summary "tighest possible" bad -> rather, illustrates what ideal direction can do - summary "tighest possible" bad -> rather, illustrates what ideal direction can do
%```{admonition} Looking ahead
%:class: tip
%Below, we'll proceed in the following steps:
%- we'll first show the problems with regular gradient descent, especially for functions that combine small and large scales,
%- a central insight will be that an _inverse gradient_ is a lot more meaningful than the regular one,
%- finally, we'll show how to use inverse functions (and especially inverse PDE solvers) to compute a very accurate update that includes higher-order terms.
%```
```{admonition} Looking ahead
:class: tip
Below, we'll proceed in the following steps:
- we'll first show the problems with regular gradient descent, especially for functions that combine small and large scales,
- a central insight will be that an _inverse gradient_ is a lot more meaningful than the regular one,
- finally, we'll show how to use inverse functions (and especially inverse PDE solvers) to compute a very accurate update that includes higher-order terms.
```
## Overview
All NNs of the previous chapters were trained with gradient descent (GD) via backpropagation, GD and hence backpropagation was also employed for the PDE solver (_simulator_) $\mathcal P$, computing the composite gradient
$(\partial L / \partial x)^T$ for the loss function $L$:
$$
\Big( \frac{\partial L}{\partial x} \Big)^T = \Big( \frac{\partial L}{\partial \mathcal P(x)} \Big)^T
\Big( \frac{\partial \mathcal P(x)}{\partial x} \Big)^T
$$
In the field of classical optimization, techniques such as Newton's method or BFGS variants are commonly used to optimize numerical processes since they can offer better convergence speed and stability.
These methods likewise employ gradient information, but substantially differ from GD in the way they
compute the update step, typically via higher order derivatives.
The PG which we'll derive below can take into account nonlinearities to produce better optimization updates when an (full or approximate) inverse simulator is available.
In contrast to classic optimization techniques, we show how a differentiable or invertible physics
simulator can be leveraged to compute the PG without requiring higher-order derivatives of the simulator.
In the following, we will stop using GD for everything, and instead use the aforementioned PGs for the simulator.
This update is combined with a GD based step for updating the weights in the NNs.
This setup, consisting of two fundamentally different optimization schemes, will result in an improved end-to-end training.
```{figure} resources/placeholder.png
---
height: 220px
name: pg-training
---
TODO, visual overview of PG training
```
![Divider](resources/divider3.jpg) ![Divider](resources/divider3.jpg)
## Traditional optimization methods ## Overview - Traditional optimization methods
We'll start by revisiting the most commonly used optimization methods -- gradient descent (GD) and quasi-Newton methods -- and describe their fundamental limits and drawbacks on a theoretical level. As before, let $L(x)$ be a scalar loss function, subject to minimization. The goal is to compute a step in terms of the input parameters $x$ , denoted by $\Delta x$. The different versions of $\Delta x$ will be denoted by a subscript.
All NNs of the previous chapters were trained with gradient descent (GD) via backpropagation. GD with backpropagation was also employed for the PDE solver (_simulator_) $\mathcal P$. As a central quantitiy, this gives the composite gradient
$(\partial L / \partial x)^T$ of the loss function $L$:
$$
\Big( \frac{\partial L}{\partial x} \Big)^T =
\Big( \frac{\partial \mathcal P(x)}{\partial x} \Big)^T
\Big( \frac{\partial L}{\partial \mathcal P(x)} \Big)^T
$$ (loss-deriv)
We've shown that using $\partial L/\partial x$ works, but
in the field of classical optimization, other algorithms are more widely used than GD (so-called quasi-Newton methods), and they use different updates.
Hence in the following we'll revisit GD, and discuss the pros and cons of the different methods on a theoretical level. Among others, it's interesting to discuss why classical optimization algorithms aren't widely used for NN training despite having some obvious advantages.
As before, let $L(x)$ be a scalar loss function, subject to minimization. The goal is to compute a step in terms of the input parameters $x$ , denoted by $\Delta x$. The different version of $\Delta x$ will be denoted by a subscript.
Note that we exclusively consider multivariate functions, and hence all symbols represent vector-valued expressions unless specified otherwise. Note that we exclusively consider multivariate functions, and hence all symbols represent vector-valued expressions unless specified otherwise.
%techniques such as Newton's method or BFGS variants are commonly used to optimize numerical processes since they can offer better convergence speed and stability. These methods likewise employ gradient information, but substantially differ from GD in the way they compute the update step, typically via higher order derivatives.
%```{figure} resources/placeholder.png
%---
%height: 220px
%name: pg-training
%---
%TODO, visual overview of PG training
%```
### Gradient descent ### Gradient descent
@@ -106,39 +97,53 @@ $$ (GD-update)
where $\eta$ is the scalar learning rate. where $\eta$ is the scalar learning rate.
The Jacobian $\frac{\partial L}{\partial x}$ describes how the loss reacts to small changes of the input. The Jacobian $\frac{\partial L}{\partial x}$ describes how the loss reacts to small changes of the input.
% Surprisingly, this very widely used update has a number of undesirable properties that we'll highlight in the following. Note that we've naturally applied this update in supervised settings such as {doc}`supervised-airfoils`, but we've also used it in the differentiable physics approaches. E.g., in {doc}`diffphys-code-sol` we've computed the derivative of the fluid solver. In the latter case, we've still only updated the NN parameters, but the fluid solver Jacobian was part of {eq}`GD-update`, as shown in {eq}`loss-deriv`.
Surprisingly, this very widely used construction has a number of undesirable properties.
**Units** 📏 **Units** 📏
A first indicator that something is amiss with GD is that it inherently misrepresents dimensions. A first indicator that something is amiss with GD is that it inherently misrepresents dimensions.
Assume two parameters $x_1$ and $x_2$ have different physical units. Assume two parameters $x_1$ and $x_2$ have different physical units.
Then the GD parameter updates scale with the inverse of these units because the parameters appear in the denominator for the GD update above. Then the GD parameter updates scale with the inverse of these units because the parameters appear in the denominator for the GD update above ($\cdots / \partial x$).
The learning rate $\eta$ could compensate for this discrepancy but since $x_1$ and $x_2$ have different units, there exists no single $\eta$ to produce the correct units for both parameters. The learning rate $\eta$ could compensate for this discrepancy but since $x_1$ and $x_2$ have different units, there exists no single $\eta$ to produce the correct units for both parameters.
One could argue that units aren't very important for the parameters of NNs, but nonetheless it's unnerving from a physics perspective that they're wrong.
**Function sensitivity** 🔍 **Function sensitivity** 🔍
GD has inherent problems when functions are not normalized. GD has also inherent problems when functions are not _normalized_.
Assume the range of $L(x)$ lies on a different scale than $x$. This can be illustrated with a very simple example:
Consider the function $L(x) = c \cdot x$ for example where $c \ll 1$ or $c \gg 1$. consider the function $L(x) = c \cdot x$.
Then the parameter updates of GD scale with $c$, i.e. $\Delta x_{\text{GD}} = \eta \cdot c$. Then the parameter updates of GD scale with $c$, i.e. $\Delta x_{\text{GD}} = \eta \cdot c$, and
Such behavior can occur easily in complex functions such as deep neural networks if the layers are not normalized correctly. $L(x+\Delta x_{\text{GD}})$ will even have terms on the order of $c^2$.
% If $L$ is normalized via $c=1$, everything's fine. But in practice, we'll often
For sensitive functions, i.e. _small changes_ in $x$ cause **large** changes in $L$, GD counter-intuitively produces large $\Delta x_{\text{GD}}$, causing even larger steps in $L$ (exploding gradients). have $c \ll 1$, or even worse $c \gg 1$, and then our optimization will be in trouble.
For insensitive functions where _large changes_ in the input don't change the output $L$ much, GD produces **small** updates, which can lead to the optimization coming to a standstill (that's the classic _vanishing gradients_ problem).
% More specifically, if we look at how the loss changes, the expansion around $x$ for
While normalization in combination with correct setting of the learning rate $\eta$ can be used to counteract this behavior in neural networks, these tools are not available when optimizing simulations. the update step of GD gives:
Applying normalization to a simulation anywhere but after the last solver step would destroy the simulation. $L(x+\Delta x) = L(x) + \Delta x \frac{\partial L}{\partial x} + \cdots $.
Setting the learning rate is also difficult when simulation parameters at different time steps are optimized simultaneously or when the magnitude of the simulation output varies w.r.t. the initial state. This first-order step causes a change in the loss of
$\big( L(x) - L(x+\Delta x) \big) = -\eta \cdot (\frac{\partial L}{\partial x})^2 + \mathcal O(\Delta x^2)$ . Hence the loss changes by the squared derivative, instead of being proportional to it, as one
might expect when applying SGD without much thought.
This demonstrates that
for sensitive functions, i.e. functions where _small changes_ in $x$ cause _large_ changes in $L$, GD counter-intuitively produces large $\Delta x_{\text{GD}}$. This causes even larger steps in $L$, and leads to exploding gradients.
For insensitive functions where _large changes_ in the input don't change the output $L$ much, GD produces _small_ updates, which can lead to the optimization coming to a halt. That's the classic _vanishing gradients_ problem.
Such sensitivity problems can occur easily in complex functions such as deep neural networks where the layers are typically not fully normalized.
Normalization in combination with correct setting of the learning rate $\eta$ can be used to counteract this behavior in NNs to some extent, but these tools are not available when optimizing physics simulations.
Applying normalization to a simulation anywhere but after the last solver step would destroy the state of the simulation.
Adjusting the learning rate is also difficult in practice, e.g. when simulation parameters at different time steps are optimized simultaneously or when the magnitude of the simulation output varies w.r.t. the initial state.
**Convergence near optimum** 💎 **Convergence near optimum** 💎
The loss landscape of any differentiable function necessarily becomes flat close to an optimum Finally, the loss landscape of any differentiable function necessarily becomes flat close to an optimum,
(the gradient approaches zero upon convergence). as the gradient approaches zero upon convergence.
Therefore $\Delta x_{\text{GD}} \rightarrow 0$ as the optimum is approached, resulting in slow convergence. Therefore $\Delta x_{\text{GD}} \rightarrow 0$ as the optimum is approached, resulting in slow convergence.
This is an important point, and we will revisit it below. It's also somewhat surprising at first, but it can actually This is an important point, and we will revisit it below. It's also somewhat surprising at first, but it can actually
stabilize the training. On the other hand, it also makes the learning process difficult to control. stabilize the training. On the other hand, it makes the learning process difficult to control.
@@ -158,7 +163,7 @@ This construction solves some of the problems of gradient descent from above, bu
Quasi-Newton methods definitely provide a much better handling of physical units than GD. Quasi-Newton methods definitely provide a much better handling of physical units than GD.
The quasi-Newton update from equation {eq}`quasi-newton-update` The quasi-Newton update from equation {eq}`quasi-newton-update`
produces the correct units for all parameters to be optimized, $\eta$ can stay dimensionless. produces the correct units for all parameters to be optimized. As a consequence, $\eta$ can stay dimensionless.
**Convergence near optimum** 💎 **Convergence near optimum** 💎
@@ -173,8 +178,8 @@ So far, quasi-Newton methods address both shortcomings of GD.
However, similar to GD, the update of an intermediate space still depends on all functions before that. However, similar to GD, the update of an intermediate space still depends on all functions before that.
This behavior stems from the fact that the Hessian of a function composition carries non-linear terms of the gradient. This behavior stems from the fact that the Hessian of a function composition carries non-linear terms of the gradient.
Consider a function composition $L(z(x))$, with $L$ as above, and an additional function $z(x)$. Consider a function composition $L(y(x))$, with $L$ as above, and an additional function $y(x)$.
Then the Hessian $\frac{d^2L}{dx^2} = \frac{\partial^2L}{\partial z^2} \left( \frac{\partial z}{\partial x} \right)^2 + \frac{\partial L}{\partial z} \cdot \frac{\partial^2 z}{\partial x^2}$ depends on the square of the inner gradient $\frac{\partial z}{\partial x}$. Then the Hessian $\frac{d^2L}{dx^2} = \frac{\partial^2L}{\partial y^2} \left( \frac{\partial y}{\partial x} \right)^2 + \frac{\partial L}{\partial y} \cdot \frac{\partial^2 y}{\partial x^2}$ depends on the square of the inner gradient $\frac{\partial y}{\partial x}$.
This means that the Hessian is influenced by the _later_ derivatives of a backpropagation pass, This means that the Hessian is influenced by the _later_ derivatives of a backpropagation pass,
and as a consequence, the update of any latent space is unknown during the computation of the gradients. and as a consequence, the update of any latent space is unknown during the computation of the gradients.
@@ -216,18 +221,20 @@ are still a very active research topic, and hence many extensions have been prop
As a first step towards _physical_ gradients, we introduce _inverse_ gradients (IGs), As a first step towards _physical_ gradients, we introduce _inverse_ gradients (IGs),
which already solve many of the aforementioned problems. Unfortunately, they come with their own set of problems, which is why they only represent an intermediate step. which already solve many of the aforementioned problems. Unfortunately, they come with their own set of problems, which is why they only represent an intermediate step.
Instead of $L$ (which is scalar), let's consider a general, potentially non-scalar function $z(x)$. We define the update Instead of $L$ (which is scalar), let's consider a general, potentially non-scalar function $y(x)$.
This will typically be the physical simulator later on, but to keep things general we'll call it $y$ for now.
We define the update
$$ $$
\Delta x_{\text{IG}} = \frac{\partial x}{\partial z} \cdot \Delta z. \Delta x_{\text{IG}} = \frac{\partial x}{\partial y} \cdot \Delta y.
$$ (IG-def) $$ (IG-def)
to be the IG update. to be the IG update.
Here, the Jacobian $\frac{\partial x}{\partial z}$, which is similar to the inverse of the GD update above, encodes how the inputs must change in order to obtain a small change $\Delta z$ in the output. Here, the Jacobian $\frac{\partial x}{\partial y}$, which is similar to the inverse of the GD update above, encodes how the inputs must change in order to obtain a small change $\Delta y$ in the output.
% %
The crucial step is the inversion, which of course requires the Jacobian matrix to be invertible (a drawback we'll get back to below). However, if we can invert it, this has some very nice properties. The crucial step is the inversion, which of course requires the Jacobian matrix to be invertible (a drawback we'll get back to below). However, if we can invert it, this has some very nice properties.
Note that instead of using a learning rate, here the step size is determined by the desired increase or decrease of the value of the output, $\Delta z$. Thus, we need to choose a $\Delta z$ instead of an $\eta$. This $\Delta z$ will show up frequently in the following equations, and make them look quite different to the ones above at first sight. Effectively, $\Delta z$ plays the same role as the learning rate, i.e., it controls the step size of the optimization. Note that instead of using a learning rate, here the step size is determined by the desired increase or decrease of the value of the output, $\Delta y$. Thus, we need to choose a $\Delta y$ instead of an $\eta$. This $\Delta y$ will show up frequently in the following equations, and make them look quite different to the ones above at first sight. Effectively, $\Delta y$ plays the same role as the learning rate, i.e., it controls the step size of the optimization.
@@ -235,7 +242,7 @@ Note that instead of using a learning rate, here the step size is determined by
**Positive Aspects** **Positive Aspects**
IGs scale with the inverse derivative. Hence the updates are automatically of the same units as the parameters without requiring an arbitrary learning rate: $\frac{\partial x}{\partial z}$ times $\Delta z$ has the units of $x$. IGs scale with the inverse derivative. Hence the updates are automatically of the same units as the parameters without requiring an arbitrary learning rate: $\frac{\partial x}{\partial y}$ times $\Delta y$ has the units of $x$.
% **Function sensitivity** % **Function sensitivity**
@@ -249,8 +256,8 @@ IGs show the opposite behavior of GD close to an optimum: they typically produce
% **Consistency in function compositions** % **Consistency in function compositions**
Additionally, IGs are consistent in function composition. Additionally, IGs are consistent in function composition.
The change in $x$ is $\Delta x_{\text{IG}} = \Delta L \cdot \frac{\partial x}{\partial z} \frac{\partial z}{\partial L}$ and the approximate change in $z$ is $\Delta z = \Delta L \cdot \frac{\partial z}{\partial x} \frac{\partial x}{\partial z} \frac{\partial z}{\partial L} = \Delta L \frac{\partial z}{\partial L}$. The change in $x$ is $\Delta x_{\text{IG}} = \Delta L \cdot \frac{\partial x}{\partial y} \frac{\partial y}{\partial L}$ and the approximate change in $y$ is $\Delta y = \Delta L \cdot \frac{\partial y}{\partial x} \frac{\partial x}{\partial y} \frac{\partial y}{\partial L} = \Delta L \frac{\partial y}{\partial L}$.
% In the example in table~\ref{tab:function-composition-example}, the change $\Delta z$ is the same no matter what space is used as optimization target. % In the example in table~\ref{tab:function-composition-example}, the change $\Delta y$ is the same no matter what space is used as optimization target.
The change in intermediate spaces is independent of their respective dependencies, at least up to first order. The change in intermediate spaces is independent of their respective dependencies, at least up to first order.
Consequently, the change to these spaces can be estimated during backpropagation, before all gradients have been computed. Consequently, the change to these spaces can be estimated during backpropagation, before all gradients have been computed.
@@ -261,11 +268,11 @@ Note that even Newton's method with its inverse Hessian didn't fully get this ri
So far so good. So far so good.
The above properties make the advantages of IGs clear, but we're not done, unfortunately. There are strong limitations to their applicability. The above properties make the advantages of IGs clear, but we're not done, unfortunately. There are strong limitations to their applicability.
% %
The IG $\frac{\partial x}{\partial z}$ is only well-defined for square Jacobians, i.e. for functions $z$ with the same inputs and output dimensions. The IG $\frac{\partial x}{\partial y}$ is only well-defined for square Jacobians, i.e. for functions $y$ with the same inputs and output dimensions.
In optimization, however, the input is typically high-dimensional while the output is a scalar objective function. In optimization, however, the input is typically high-dimensional while the output is a scalar objective function.
% %
And, somewhat similar to the Hessians of quasi-Newton methods, And, somewhat similar to the Hessians of quasi-Newton methods,
even when the $\frac{\partial z}{\partial x}$ is square, it may not be invertible. even when the $\frac{\partial y}{\partial x}$ is square, it may not be invertible.
Thus, we now consider the fact that inverse gradients are linearizations of inverse functions and show that using inverse functions provides additional advantages while retaining the same benefits. Thus, we now consider the fact that inverse gradients are linearizations of inverse functions and show that using inverse functions provides additional advantages while retaining the same benefits.
@@ -281,11 +288,11 @@ As long as the physical process does _not destroy_ information, the Jacobian is
In fact, it is believed that information in our universe cannot be destroyed so any physical process could in theory be inverted as long as we have perfect knowledge of the state. In fact, it is believed that information in our universe cannot be destroyed so any physical process could in theory be inverted as long as we have perfect knowledge of the state.
While evaluating the IGs directly can be done through matrix inversion or taking the derivative of an inverse simulator, we now consider what happens if we use the inverse simulator directly in backpropagation. While evaluating the IGs directly can be done through matrix inversion or taking the derivative of an inverse simulator, we now consider what happens if we use the inverse simulator directly in backpropagation.
Let $z = \mathcal P(x)$ be a forward simulation, and $\mathcal P(z)^{-1}=x$ its inverse (we assume it exists for now, but below we'll relax that assumption). Let $y = \mathcal P(x)$ be a forward simulation, and $\mathcal P(y)^{-1}=x$ its inverse (we assume it exists for now, but below we'll relax that assumption).
Equipped with the inverse we now define an update that we'll call the **physical gradient** (PG) {cite}`holl2021pg` in the following as Equipped with the inverse we now define an update that we'll call the **physical gradient** (PG) {cite}`holl2021pg` in the following as
$$ $$
\frac{\Delta x_{\text{PG}} }{\Delta z} \equiv \big( \mathcal P^{-1} (z_0 + \Delta z) - x_0 \big) \frac{\Delta x_{\text{PG}} }{\Delta y} \equiv \big( \mathcal P^{-1} (y_0 + \Delta y) - x_0 \big)
$$ (PG-def) $$ (PG-def)
@@ -294,24 +301,24 @@ $$ (PG-def)
% add ? $ / \Delta z $ on the right!? the above only gives $\Delta x$, see below % add ? $ / \Delta z $ on the right!? the above only gives $\Delta x$, see below
Note that this PG is equal to the IG from the section above up to first order, but contains nonlinear terms, i.e. Note that this PG is equal to the IG from the section above up to first order, but contains nonlinear terms, i.e.
$ \Delta x_{\text{PG}} / \Delta z = \frac{\partial x}{\partial z} + \mathcal O(\Delta z^2) $. $ \Delta x_{\text{PG}} / \Delta y = \frac{\partial x}{\partial y} + \mathcal O(\Delta y^2) $.
% %
The accuracy of the update also depends on the fidelity of the inverse function $\mathcal P^{-1}$. The accuracy of the update also depends on the fidelity of the inverse function $\mathcal P^{-1}$.
We can define an upper limit to the error of the local inverse using the local gradient $\frac{\partial x}{\partial z}$. We can define an upper limit to the error of the local inverse using the local gradient $\frac{\partial x}{\partial y}$.
In the worst case, we can therefore fall back to the regular gradient. In the worst case, we can therefore fall back to the regular gradient.
% We now show that these terms can help produce more stable updates than the IG alone, provided that $\mathcal P_{(x_0,z_0)}^{-1}$ is a sufficiently good approximation of the true inverse. % We now show that these terms can help produce more stable updates than the IG alone, provided that $\mathcal P_{(x_0,z_0)}^{-1}$ is a sufficiently good approximation of the true inverse.
% Let $\mathcal P^{-1}(z)$ be the true inverse function to $\mathcal P(x)$, assuming that $\mathcal P$ is fully invertible. % Let $\mathcal P^{-1}(z)$ be the true inverse function to $\mathcal P(x)$, assuming that $\mathcal P$ is fully invertible.
The intuition for why the PG update is a good one is that when The intuition for why the PG update is a good one is that when
applying the update $\Delta x_{\text{PG}} = \mathcal P^{-1}(z_0 + \Delta z) - x_0$ it will produce $\mathcal P(x_0 + \Delta x) = z_0 + \Delta z$ exactly, despite $\mathcal P$ being a potentially highly nonlinear function. applying the update $\Delta x_{\text{PG}} = \mathcal P^{-1}(y_0 + \Delta y) - x_0$ it will produce $\mathcal P(x_0 + \Delta x) = y_0 + \Delta y$ exactly, despite $\mathcal P$ being a potentially highly nonlinear function.
When rewriting this update in the typical gradient format, $\frac{\Delta x_{\text{PG}}}{\Delta z}$ replaces the gradient from the IG update above, and gives $\Delta x$. When rewriting this update in the typical gradient format, $\frac{\Delta x_{\text{PG}}}{\Delta y}$ replaces the gradient from the IG update above, and gives $\Delta x$.
**Fundamental theorem of calculus** **Fundamental theorem of calculus**
To more clearly illustrate the advantages in non-linear settings, we To more clearly illustrate the advantages in non-linear settings, we
apply the fundamental theorem of calculus to rewrite the ratio $\Delta x_{\text{PG}} / \Delta z$ from above. This gives, apply the fundamental theorem of calculus to rewrite the ratio $\Delta x_{\text{PG}} / \Delta y$ from above. This gives,
% \begin{equation} \label{eq:avg-grad} % \begin{equation} \label{eq:avg-grad}
@@ -323,10 +330,10 @@ apply the fundamental theorem of calculus to rewrite the ratio $\Delta x_{\text{
% focused on 1D for simplicity. Likewise, by integrating over $z$ we can obtain: % focused on 1D for simplicity. Likewise, by integrating over $z$ we can obtain:
$\begin{aligned} $\begin{aligned}
\frac{\Delta x_{\text{PG}}}{\Delta z} = \frac{\int_{z_0}^{z_0+\Delta z} \frac{\partial x}{\partial z} \, dz}{\Delta z} \frac{\Delta x_{\text{PG}}}{\Delta y} = \frac{\int_{y_0}^{y_0+\Delta y} \frac{\partial x}{\partial y} \, dy}{\Delta y}
\end{aligned}$ \end{aligned}$
Here the expressions inside the integral is the local gradient, and we assume it exists at all points between $z_0$ and $z_0+\Delta z_0$. Here the expressions inside the integral is the local gradient, and we assume it exists at all points between $y_0$ and $y_0+\Delta y_0$.
The local gradients are averaged along the path connecting the state before the update with the state after the update. The local gradients are averaged along the path connecting the state before the update with the state after the update.
The whole expression is therefore equal to the average gradient of $\mathcal P$ between the current $x$ and the estimate for the next optimization step $x_0 + \Delta x_{\text{PG}}$. The whole expression is therefore equal to the average gradient of $\mathcal P$ between the current $x$ and the estimate for the next optimization step $x_0 + \Delta x_{\text{PG}}$.
This effectively amounts to _smoothing the objective landscape_ of an optimization by computing updates that can take nonlinearities of $\mathcal P$ into account. This effectively amounts to _smoothing the objective landscape_ of an optimization by computing updates that can take nonlinearities of $\mathcal P$ into account.
@@ -340,29 +347,29 @@ The equations naturally generalize to higher dimensions by replacing the integra
### Global and local inverse functions ### Global and local inverse functions
Let $\mathcal P$ be a function with a square Jacobian and $z = \mathcal P(x)$. Let $\mathcal P$ be a function with a square Jacobian and $y = \mathcal P(x)$.
A global inverse function $\mathcal P^{-1}$ is defined only for bijective $\mathcal P$. A global inverse function $\mathcal P^{-1}$ is defined only for bijective $\mathcal P$.
If the inverse exists, it can find $x$ for any $z$ such that $z = \mathcal P(x)$. If the inverse exists, it can find $x$ for any $y$ such that $y = \mathcal P(x)$.
Instead of using this "perfect" inverse $\mathcal P^{-1}$ directly, we'll in practice often use a local inverse Instead of using this "perfect" inverse $\mathcal P^{-1}$ directly, we'll in practice often use a local inverse
$\mathcal P_{(x_0,z_0)}^{-1}(z)$, defined at the point $(x_0, z_0)$. This local inverse can be $\mathcal P_{(x_0,y_0)}^{-1}(y)$, defined at the point $(x_0, y_0)$. This local inverse can be
easier to obtain, as it only needs to exist near a given $z_0$, and not for all $z$. easier to obtain, as it only needs to exist near a given $y_0$, and not for all $y$.
For $\mathcal P^{-1}$ to exist $\mathcal P$ would need to be globally invertible. For $\mathcal P^{-1}$ to exist $\mathcal P$ would need to be globally invertible.
By contrast, a \emph{local inverse}, defined at point $(x_0, z_0)$, only needs to be accurate in the vicinity of that point. By contrast, a \emph{local inverse}, defined at point $(x_0, y_0)$, only needs to be accurate in the vicinity of that point.
If a global inverse $\mathcal P^{-1}(z)$ exists, the local inverse approximates it and matches it exactly as $z \rightarrow z_0$. If a global inverse $\mathcal P^{-1}(y)$ exists, the local inverse approximates it and matches it exactly as $y \rightarrow y_0$.
More formally, $\lim_{z \rightarrow z_0} \frac{\mathcal P^{-1}_{(x_0, z_0)}(z) - P^{-1}(z)}{|z - z_0|} = 0$. More formally, $\lim_{y \rightarrow y_0} \frac{\mathcal P^{-1}_{(x_0, y_0)}(y) - P^{-1}(y)}{|y - y_0|} = 0$.
Local inverse functions can exist, even when a global inverse does not. Local inverse functions can exist, even when a global inverse does not.
Non-injective functions can be inverted, for example, by choosing the closest $x$ to $x_0$ such that $\mathcal P(x) = z$. Non-injective functions can be inverted, for example, by choosing the closest $x$ to $x_0$ such that $\mathcal P(x) = y$.
With the local inverse, the PG is defined as With the local inverse, the PG is defined as
$$ $$
\frac{\Delta x_{\text{PG}}}{\Delta z} \equiv \big( \mathcal P_{(x_0,z_0)}^{-1} (z_0 + \Delta z) - x_0 \big) / \Delta z \frac{\Delta x_{\text{PG}}}{\Delta y} \equiv \big( \mathcal P_{(x_0,y_0)}^{-1} (y_0 + \Delta y) - x_0 \big) / \Delta y
$$ (local-PG-def) $$ (local-PG-def)
For differentiable functions, a local inverse is guaranteed to exist by the inverse function theorem as long as the Jacobian is non-singular. For differentiable functions, a local inverse is guaranteed to exist by the inverse function theorem as long as the Jacobian is non-singular.
That is because the inverse Jacobian $\frac{\partial x}{\partial z}$ itself is a local inverse function, albeit not the most accurate one. That is because the inverse Jacobian $\frac{\partial x}{\partial y}$ itself is a local inverse function, albeit not the most accurate one.
Even when the Jacobian is singular (because the function is not injective, chaotic or noisy), we can usually find good local inverse functions. Even when the Jacobian is singular (because the function is not injective, chaotic or noisy), we can usually find good local inverse functions.