diff --git a/old-phiflow1.md b/old-phiflow1.md new file mode 100644 index 0000000..00515e3 --- /dev/null +++ b/old-phiflow1.md @@ -0,0 +1,4 @@ +# Old Phiflow1 Examples + +Remove sometime... + diff --git a/physgrad-comparison.ipynb b/physgrad-comparison.ipynb new file mode 100644 index 0000000..a5e7cb2 --- /dev/null +++ b/physgrad-comparison.ipynb @@ -0,0 +1,675 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Optimizations with Physical Gradients\n", + "\n", + "Test equations: properly define z in R2 to R2 etc.
\n", + "
\n", + " \n", + "$\\mathbf{z}(\\mathbf{x}) = \\mathbf{z}(x_0,x_1) = [x_0 \\ x_1^2]^T$
\n", + "$L(\\mathbf{z}) = |\\mathbf{z}|^2 = z_0^2 + z_1^2$
\n", + "\n", + "(Note, to be in sync with python arrays, indices start at 0 here.)\n", + "\n", + "..." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "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", + "output_type": "stream", + "text": [ + "Some test calls:\n" + ] + }, + { + "data": { + "text/plain": [ + "(DeviceArray([3., 9.], dtype=float32),\n", + " DeviceArray(90., dtype=float32),\n", + " DeviceArray(90., dtype=float32),\n", + " DeviceArray(0, dtype=int32))" + ] + }, + "execution_count": 1, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "#import numpy as np\n", + "\n", + "import jax\n", + "import jax.numpy as np\n", + "import numpy as onp\n", + "\n", + "\n", + "# \"physics\" function z\n", + "def fun_z(x):\n", + " return np.array( [x[0], x[1]*x[1]] )\n", + "\n", + "# simple L2 loss\n", + "def fun_L(z):\n", + " #return z[0]*z[0] + z[1]*z[1] # \"manual version\"\n", + " return np.sum( np.square(z) )\n", + "\n", + "# composite function with L & z\n", + "def fun(x):\n", + " return fun_L(fun_z(x))\n", + "\n", + "\n", + "x = np.asarray([3,3], dtype=np.float32)\n", + "\n", + "print(\"Some test calls:\") \n", + "fun_z(x) , fun_L( fun_z(x) ), fun(x), fun([0,0])" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Jacobian L(z): [ 6. 18.]\n", + "Jacobian z(x): \n", + "[[1. 0.]\n", + " [0. 6.]]\n", + "Gradient for full L(x): [ 6. 108.]\n", + "\n", + "Sanity check with inverse Jacobian of z, this should give x again: [3. 3.]\n" + ] + } + ], + "source": [ + "\n", + "# this works:\n", + "print(\"Jacobian L(z): \" + format(jax.grad(fun_L)(fun_z(x))) )\n", + "\n", + "# the following would give an error as z (and hence fun_z) is not scalar\n", + "#jax.grad(fun_z)(x) \n", + "\n", + "# computing the jacobian of z is a valid operation:\n", + "J = jax.jacobian(fun_z)(x)\n", + "print( \"Jacobian z(x): \\n\" + format(J) ) \n", + "\n", + "# the following also gives error, JAX grad needs a single function object\n", + "#jax.grad( fun_L(fun_z) )(x) \n", + "\n", + "# instead use composite 'fun' from above\n", + "print(\"Gradient for full L(x): \" + format( jax.grad(fun)(x) ) )\n", + "\n", + "print( \"\\nSanity check with inverse Jacobian of z, this should give x again: \" + format(np.linalg.solve(J, np.matmul(J,x) )) )\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This already shows a problem, the gradient mostly points along $x_2$ with 108! \n", + "We know both dimensions should go towards zero, this points to problems..." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Gradient descent\n", + "\n", + "Update step:
\n", + "\n", + "$\\Delta \\mathbf{x} = \n", + "- \\eta ( J_{L} J_{\\mathbf{z}} )^T\n", + "=\n", + "- \\eta ( \\frac{\\partial L }{ \\partial \\mathbf{z} } \\frac{\\partial \\mathbf{z} }{ \\partial \\mathbf{x} } )^T$\n", + "\n", + "with step size $\\eta$, \n", + "start at $x=[3,3]$\n", + "\n", + "..." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "GD iter 0: [2.94 1.9200001]\n", + "GD iter 1: [2.8812 1.6368846]\n", + "GD iter 2: [2.823576 1.4614503]\n", + "GD iter 3: [2.7671044 1.3365935]\n", + "GD iter 4: [2.7117622 1.2410815]\n", + "GD iter 5: [2.657527 1.1646168]\n", + "GD iter 6: [2.6043763 1.1014326]\n", + "GD iter 7: [2.5522888 1.0479842]\n", + "GD iter 8: [2.501243 1.0019454]\n", + "GD iter 9: [2.4512184 0.96171147]\n" + ] + } + ], + "source": [ + "x = np.asarray([3.,3.])\n", + "eta = 0.01\n", + "history = [x]; updates = []\n", + "\n", + "for i in range(10):\n", + " G = jax.grad(fun)(x)\n", + " x += -eta * G\n", + " history.append(x); updates.append(G)\n", + " print( \"GD iter %d: \"%i + format(x) )\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let's look at the progression of this...\n", + "Target shown in black" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Text(0, 0.5, 'x1')" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAFtCAYAAADrr7rKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAcRElEQVR4nO3df5RcZZ3n8fc3TchxId16kEkQgsnIWUThDCBHDesAyrAqkmh2ZfWIcyawnhHHnKjr6poFXWcWJjN6FnEWRlZ0gqxzHP9ZEBRFGSd4ViIqyq81srKCA0iCona3uxIz5Lt/3Gotiqrqru6qvlX9vF/n3NPUc59b95tr+enbz33q3shMJEllWFZ3AZKkxWPoS1JBDH1JKoihL0kFMfQlqSCGviQVxNCXpIIY+pJUkIPqLmCxRUQAzwGm665FkvpoJfDjnOUbt8WFPlXgP1x3EZI0AEcBj3TrUGLoTwM89NBDjI+P112LJC3Y1NQUa9asgTmMYJQY+gCMj48b+pKK44VcSSqIoS9JBTH0Jakghr4kFcTQl6SCGPqSNEQG/TDDWkM/It4WEXdHxFRj2RURr55lm3Mj4vsR8URE3BMRZy9WvZI0CNPTsHUrrFsHa9ZUP7durdr7Lep8Rm5EbACeBH4ABPBHwHuAkzLzf7XpfyrwNWAb8HngTcB/AE7OzHvnuM9xYHJyctJ5+pJqNz0N69fD7t1w4MBv25ctg+OOg127YOXK7u8xNTXFxMQEwERmTnXrW+uZfmbemJk3ZeYPMvN/Z+ZFwC+Bl3bY5B3AlzLzw5m5OzPfD3wH2LJYNUtSP1100dMDH6rXu3fDxRf3d39DM6YfEWMR8UbgEGBXh27rgVta2m5utHd63xURMT6zUN2USJKGwo03Pj3wZxw4ADfc0N/91R76EXFCRPwS2AdcBWzKzO916L4a2NvStrfR3sk2YLJp8WZrkoZCJuzf373P/v39vbhbe+gD9wEnAi8BPgZ8KiJe0Mf33w5MNC1H9fG9JWneImD58u59li+v+vVL7aGfmb/OzPsz847M3AbcRTV2384eYFVL26pGe6f335eZUzML3kdf0hDZsKG6aNvOsmWwcWN/91d76LexDFjRYd0u4MyWtrPofA1AkobapZdWs3Rag39m9s4ll/R3f3XP098eEadFxNrG2P524Azgbxvrr220zfgo8KqIeHdEPD8iPgicAlyx2LVLUj+sXFlNy9yyBdauhSOPrH5u2TK36Zq9qnue/iepztyPoLrIejfwl5n5lcb6ncCDmbm5aZtzgUuAtVTz+9+bmTf1sE/n6UsaWpm9j+H3Mk+/1tCvg6EvaakZmS9nSZIWl6EvSQUx9CWpIIa+JBXE0Jekghj6klQQQ1+SCmLoS1JBDH1JKoihL0kFMfQlqSCGviQVxNCXpIIY+pJUEENfkgpi6EtSQQx9SSqIoS9JBTH0Jakghr4kFcTQl6SCGPqSVBBDX5IKYuhLUkEMfUkqiKEvSQUx9CWpIIa+JBXE0Jekghj6klQQQ1+SCmLoS1JBDH1JKoihL0kFMfQlqSC1hn5EbIuIb0XEdEQ8FhHXR8Sxs2yzOSKyZXlisWqWpFFW95n+6cCVwEuBs4DlwJcj4pBZtpsCjmhanjvIIiVpqTiozp1n5quaX0fEZuAx4EXA17pvmnsGWJokLUl1n+m3mmj8/Nks/Q6NiB9FxEMR8bmIeGGnjhGxIiLGZxZgZd+qlaQRMzShHxHLgMuBr2fmvV263gdcALwWeDPVv+G2iDiqQ/9twGTT8nC/apakUROZWXcNAETEx4BXAy/LzDkHc0QsB3YDn8nM97dZvwJY0dS0Enh4cnKS8fHxBVYtSfWbmppiYmICYCIzp7r1rXVMf0ZEXAGcA5zWS+ADZOb+iPgucEyH9fuAfU37WkipkjTS6p6yGY3A3wS8IjMfmMd7jAEnAI/2uz5JWmrqPtO/EngT1fj8dESsbrRPZuavACLiWuCRzNzWeP0B4BvA/cAzgfdQTdn8xOKWLkmjp+7Qf1vj586W9vOBaxr/fTRwoGnds4CrgdXAz4E7gFMz83sDq1KSloihuZC7WBrTNie9kCtpqejlQu7QTNmUJA2eoS9JBTH0Jakghr4kFcTQl6SCGPqSVBBDX5IKYuhLUkEMfUkqiKEvSQUx9CWpIIa+JBXE0Jekghj6klQQQ1+SCmLoS1JBDH1JKoihL0kFMfQlqSCGviQVxNCXpIIY+pJUEENfkgpi6EtSQQx9SSqIoS9JBTH0Jakghr4kFcTQl6SCGPqSVBBDX5IKYuhLUkEMfUkqiKEvSQUx9CWpILWGfkRsi4hvRcR0RDwWEddHxLFz2O7ciPh+RDwREfdExNmLUa8kjbq6z/RPB64EXgqcBSwHvhwRh3TaICJOBT4DfBI4CbgeuD4ijh94tZI04iIz667hNyLicOAx4PTM/FqHPp8FDsnMc5ravgHcmZkXzmEf48Dk5OQk4+PjfapckuozNTXFxMQEwERmTnXrW/eZfquJxs+fdemzHrilpe3mRvvTRMSKiBifWYCVCy9TkkbT0IR+RCwDLge+npn3dum6Gtjb0ra30d7ONmCyaXl4YZVK0ugamtCnGts/Hnhjn993O9VfEDPLUX1+f0kaGQfVXQBARFwBnAOclpmznYnvAVa1tK1qtD9NZu4D9jXtawGVSlL/ZMJiR1LdUzajEfibgFdk5gNz2GwXcGZL21mNdkkaatPTsHUrrFsHa9ZUP7durdoXQ62zdyLir4E3Aa8F7mtaNZmZv2r0uRZ4JDO3NV6fCtwKvA/4AtVw0H8ETp7lWsDMPp29I6kW09Owfj3s3g0HDvy2fdkyOO442LULVs5jqskozd55G9U4+07g0ablDU19jgaOmHmRmbdR/aL4Y+Au4PXA6+YS+JJUp4suenrgQ/V69264+OLB1zBU8/QXg2f6kuqybh08+GDn9WvXwgNzGeRuMUpn+pJUhEzYv797n/37q36DZOhL0iKIgOXLu/dZvnzws3kMfUlaJBs2VBdt21m2DDZuHHwNhr4kLZJLL61m6bQG/8zsnUsuGXwNhr4kLZKVK6tpmVu2VBdtjzyy+rlly/yna/bK2TuSVJN+fSPX2TuSNALquCuMoS9JBTH0Jakghr4kFcTQl6SCGPqSVBBDX5IKYuhLUkEMfUkqiKEvSQUx9CWpIIa+JBXE0Jekghj6klQQQ1+SCmLoS9KADdNjSwx9SRqA6WnYuhXWrYM1a6qfW7dW7XXyyVmS1GfT07B+PezeDQcO/LZ95lm4/X40ok/OkqQaXXTR0wMfqte7d8PFF9dTFxj6ktR3N9749MCfceAA3HDD4tbTzNCXpD7KhP37u/fZv7++i7uGviT1UQQsX969z/Ll9TwUHQx9Seq7DRuqi7btLFsGGzcubj1P2X99u5akpenSS6tZOq3BPzN755JL6qkLDH1J6ruVK6tpmVu2wNq1cOSR1c8tW/o/XbNXztOXpAHLHOwYvvP0JWmI1HXRth1DX5IK0rfQj4jjIuKHPW5zWkTcGBE/joiMiNfN0v+MRr/WZfWCipekQvTzTP9g4Lk9bnMIcBfw9h63OxY4oml5rMftJalIB821Y0RcNkuXw3vdeWZ+Efhi4/172fSxzPxFr/uTpNLNOfSBdwB3Ap2uDB+64Grm7s6IWAHcC3wwM7++iPuWpI4GPVNnoXoJ/fuBj2Tmp9utjIgTgTv6UVQXjwIXAt8GVgBvAXZGxEsy8zsd6lrR6Dujxhmykpai6enqzpo33ljdV2f58upbuZdeWu+c/HZ6Cf1vAy8C2oY+kMBAf79l5n3AfU1Nt0XE84B3AX/YYbNtwH8aZF2SytXp3vlXXglf/Wr9X8Zq1cuF3HcDl3damZl3ZWYdU0C/CRzTZf12YKJpOWoxipJUhmG+d347cw7pzNyTmT+KiJd36hMRb+1PWT05kWrYp63M3JeZUzMLUPPDyiQtJcN87/x25nNm/qWI+HBE/ObmoRHx7Ii4EfiLXt4oIg6NiBMb1wMA1jVeH91Yvz0irm3q/86IeG1EHBMRx0fE5cArgCvn8e+QpAUZ9nvntzOf0H85sAn4VkS8ICJeQzWLZpzqrLsXpwDfbSwAlzX++88ar48Ajm7qfzDwX4B7gFuB3wP+IDP/vvd/hiQtzLDfO7+dXi7kApCZtzXOzK8CvkP1i+P9wIeyx7u3ZeZOulz8zczNLa8/BHyot4olaXA2bKgu2rYb4qn73vntzPfC6z+nOkt/GPgnqm/I/rN+FSVJo2KY753fTs+hHxHvA3YBXwGOB14MnATcHRHr+1ueJA23Yb53fjs9308/Ih4FLmjcQmGmbTnw58DWzFzRceMh4P30JQ1SHd/I7eV++j2P6QMnZOZPmxsycz/wnoj4/DzeT5KWjGG6aNtOz8M7rYHfsu7WhZUjSaNhmKZh9sKHqEjSHE1Pw9atsG4drFlT/dy6tWofFT4jV5LmoNM9dmZm6dR50dZn5EpSn43aPXY6MfQlaQ5G7R47nRj6kjSLUbzHTieGviTNYhTvsdOJoS9Jc7Bhw9NvtTBjGO+x04mhL0mzyBy9e+x0Mp9v5ErSktfuubevfCX8/u/Dl77027aNG6vAH7Z77HRi6EtSi05z8q++ujqrv/tuOPTQ0RjDb+XwjiS1mMuc/FEMfDD0Jelplsqc/HYMfUlqspTm5Ldj6EtSk6U0J78dQ1+SmmQunTn57Th7R1LxWqdnjo3BM58JP//5U4dxRm1OfjuGvqSidbtl8rOeVU3NfPLJ0ZyT346hL6lo3aZn/uIX8OY3w+WXj+4YfivH9CUVbS7TM5dK4IOhL6lgS316ZjuGvqSiLeXpme0Y+pKK0vpw88cf79x31KdntuOFXEnF6DRTp52lMD2zHUNfUjE6zdSBagjn0ENhfHzpTM9sx9CXVIxuM3Uy4bDD4Ic/XFpj+K0c05dUhAMH5jZTZ6kz9CUtWc0XbY8+Gvbs6d5/qc3UacfhHUlLUi8XbWFpztRpxzN9SUtSt4u2rZbqTJ12DH1JS05m94u2AAcdBEceCWvXwpYtsGvX0pup006toR8Rp0XEjRHx44jIiHjdHLY5IyK+ExH7IuL+iNg8+EolDbvWL1099FD3/qtWwT/+IzzwAHz0o2UEPtQ/pn8IcBfwN8D/mK1zRKwDvgBcBZwHnAl8IiIezcybB1mopOHV6/g9VBdtOz0oZSmrNfQz84vAFwFibpfMLwQeyMx3N17vjoiXAe8CDH2pUL2M30M5F23bGbXfc+uBW1rabm60txURKyJifGYBCvkjTirHDTf0FvilXLRtZ9RCfzWwt6VtLzAeEc/osM02YLJpeXhw5UlaLM1j+LON34+NwXOeU95F23bqHtNfDNuBy5per8Tgl0Zar2P4a9Ys/dsrzNWohf4eYFVL2ypgKjN/1W6DzNwH7Jt5PcdrB5KGWK9z8DduNPBnjNrwzi6qGTvNzmq0SyrEbHPwZ5Q+ft9O3fP0D42IEyPixEbTusbroxvrt0fEtU2bXAX8bkR8KCKeHxF/Avwb4COLW7mkuszlEYfLlsFzn+v4fTt1D++cAvxD0+uZsfdPAZuBI4CjZ1Zm5gMR8RqqkH8H1dj8W5yjL5UjYvZHHB59dPWlKz1d3fP0dwIdR9oyc3OHbU4aWFGSht6GDXDlle2HeEqegz8XozamL0lcemk1Vt/6jVrH8Gdn6EsaOStXVmP1W7ZUc+9LvHHafEVm1l3Domp8K3dycnKS8fHxusuR1AeZZU/JnJqaYmJiAmAiM6e69fVMX9LIKznwe2XoS1JBDH1JKoihL0kFMfQlqSCGviQVxNCXpIIY+pJUEENfkgpi6EtSQQx9SSqIoS9JBTH0Jakghr4kFcTQl6SCGPqSVBBDX5IKYuhLUkEMfUkqiKEvSQUx9CWpIIa+JBXE0Jekghj6klQQQ1+SCmLoS1JBDH1JKoihL0kFMfQlqSCGviQVxNCXpIIY+pJUkKEI/Yh4e0Q8GBFPRMTtEfHiLn03R0S2LE8sZr2SNKpqD/2IeANwGfCnwMnAXcDNEfE7XTabAo5oWp476DolaSmoPfSBfwdcnZk7MvN7wIXA/wMu6LJNZuaepmXvolQqSSOu1tCPiIOBFwG3zLRl5oHG6/VdNj00In4UEQ9FxOci4oVd9rEiIsZnFmBlv+qXpFFT95n+s4ExoPVMfS+wusM291H9FfBa4M1U/4bbIuKoDv23AZNNy8MLrFmSRlbdod+zzNyVmddm5p2ZeSvwr4CfAG/tsMl2YKJp6fTLQZKWvINq3v9PgSeBVS3tq4A9c3mDzNwfEd8Fjumwfh+wb+Z1RMyvUklaAmo908/MXwN3AGfOtEXEssbrXXN5j4gYA04AHh1EjZK0lNR9pg/VdM1PRcS3gW8C7wQOAXYARMS1wCOZua3x+gPAN4D7gWcC76GasvmJxS5ckkZN7aGfmZ+NiMOBP6O6eHsn8KqmaZhHAweaNnkWcHWj78+p/lI4tTHdU5LURWRm3TUsqsa0zcnJyUnGx8frLkeSFmxqaoqJiQmAicyc6tZ35GbvSJLmz9CXpIIY+pJUEENfkgpi6EtSQQx9SSqIoS9JBTH0Jakghr4kFcTQl6SCGPqSVBBDX5IKYuhLUkEMfUkqiKEvSQUx9CWpIIa+JBXE0Jekghj6klQQQ1+SCmLoS1JBDH1JKoihL0kFMfQlqSCGviQVxNCXpIIY+pJUEENfkgpi6EtSQQx9SSqIoS9JBTH0Jakghr4kFcTQl6SCGPqSVJChCP2IeHtEPBgRT0TE7RHx4ln6nxsR32/0vycizh5UbZnJjh072LhxI+vXr2fjxo3s2LGDzBzULiVpYA6qu4CIeANwGXAhcDvwTuDmiDg2Mx9r0/9U4DPANuDzwJuA6yPi5My8t5+1ZSabN2/m05/+NAcOHPhN+0033cTOnTu55ppriIh+7lKSBirqPmONiNuBb2XmlsbrZcBDwH/NzL9o0/+zwCGZeU5T2zeAOzPzwjnsbxyYnJycZHx8vGvfHTt28Ja3vOUpgT9jbGyMj3/841xwwQWz7VKSBmpqaoqJiQmAicyc6ta31uGdiDgYeBFwy0xbZh5ovF7fYbP1zf0bbu7UPyJWRMT4zAKsnGt91113XdvAB3jyySe57rrr5vpWkjQU6h7TfzYwBuxtad8LrO6wzeoe+28DJpuWh+da3E9+8pOu6x9//PG5vpUkDYW6Q38xbAcmmpaj5rrh4Ycf3nX9YYcdtqDCJGmx1R36PwWeBFa1tK8C9nTYZk8v/TNzX2ZOzSzA9FyL27RpE2NjY23XjY2NsWnTprm+lSQNhVpDPzN/DdwBnDnT1riQeyawq8Nmu5r7N5zVpf+8bd68mfPOO+9pwT82NsZ5553H+eef3+9dStJA1T5lk2q65qci4tvAN6mmbB4C7ACIiGuBRzJzW6P/R4FbI+LdwBeANwKnAH/c78IigmuuuYbTTz+d6667jscff5zDDjuMTZs2cf755ztdU9LIqX3KJkBEbAHeQ3Ux9k5ga2be3li3E3gwMzc39T8XuARYC/wAeG9m3jTHfc15yqYkjYJepmwORegvJkNf0lIzMvP0JUmLy9CXpIIY+pJUEENfkgpi6EtSQQx9SSrIMHw5qxZTU11nNUnSyOglz0qcp38kPdxpU5JGyFGZ+Ui3DiWGfgDPoYcbrzWspPplcdQ8tq2btdfD2utRau0rgR/nLKFe3PBO44B0/U3YTtN9dqZn+8bbsLH2elh7PQqufU79vZArSQUx9CWpIIb+3O0D/rTxc9RYez2svR7W3kVxF3IlqWSe6UtSQQx9SSqIoS9JBTH0Jakghn6TiHh7RDwYEU9ExO0R8eJZ+p8bEd9v9L8nIs5erFrb1DLn2iNic0Rky/LEYtbbVMtpEXFjRPy4Ucfr5rDNGRHxnYjYFxH3R8TmwVfato6eam/U3XrcMyJWL1LJM3Vsi4hvRcR0RDwWEddHxLFz2K72z/t8ah+Wz3tEvC0i7o6IqcayKyJePcs2fT/mhn5DRLwBuIxqutTJwF3AzRHxOx36nwp8BvgkcBJwPXB9RBy/KAU/tZaeam+YAo5oWp476Do7OISq3rfPpXNErAO+APwDcCJwOfCJiHjlgOrrpqfamxzLU4/9Y32uazanA1cCLwXOApYDX46IQzptMESf955rbxiGz/vDwPuAFwGnAF8FPhcRL2zXeWDHPDNdqmmrtwNXNL1eRnW7hvd16P9Z4PMtbd8ArhqB2jcDv6j7mLepK4HXzdLnL4F7W9r+DvjSCNR+RqPfM+s+1i11Hd6o67QufYbm8z6P2ofy896o7WfAv13MY+6ZPhARB1P99r1lpi0zDzRer++w2frm/g03d+k/EPOsHeDQiPhRRDwUER3PNobQUBz3BbozIh6NiK9ExL+ouxhgovHzZ136DOtxn0vtMGSf94gYi4g3Uv21uKtDt4Ecc0O/8mxgDNjb0r4X6DTeurrH/oMyn9rvAy4AXgu8mepzcFtEHDWoIvuo03Efj4hn1FBPLx4FLgT+dWN5CNgZESfXVVBELKMaIvt6Zt7bpeuwfN5/o4fah+bzHhEnRMQvqb5xexWwKTO/16H7QI55cXfZFGTmLprOLiLiNmA38Fbg/XXVtdRl5n1UATTjtoh4HvAu4A/rqYorgeOBl9W0/4WYU+1D9nm/j+pa1ATweuBTEXF6l+DvO8/0Kz8FngRWtbSvAvZ02GZPj/0HZT61P0Vm7ge+CxzT39IGotNxn8rMX9VQz0J9k5qOe0RcAZwDvDwzZ3uw0LB83oGea3+KOj/vmfnrzLw/M+/IzG1UEwHe0aH7QI65oU/1PwRwB3DmTFvjT8cz6Tzetqu5f8NZXfoPxDxrf4qIGANOoBp+GHZDcdz76EQW+bhH5QpgE/CKzHxgDpsNxXGfZ+2t7zFMn/dlwIoO6wZzzOu+ej0sC/AG4Angj4DjgP8G/BxY1Vh/LbC9qf+pwH7g3cDzgQ8CvwaOH4HaPwD8S+B3qaZ4fgb4FfCCGmo/lCr4TqSahfGuxn8f3Vi/Hbi2qf864P8CH2oc9z8B/gl45QjU/k6qceVjqIYlLqf6K+3MRa77r4FfUE1/XN20PKOpz1B+3udZ+1B83hufh9OAtVS/dLYDB4CzFvOYL+r/SYZ9AbYAP6K6yHI78JKmdTuBa1r6n0s1RrcPuBc4exRqBz7S1HcP1bz3k2qq+4xGYLYu1zTWXwPsbLPNdxv1/x9g8yjUDrwXuL8ROI9Tfdfg5TXU3a7mbD6Ow/p5n0/tw/J5p5pv/2CjjseoZuactdjH3FsrS1JBHNOXpIIY+pJUEENfkgpi6EtSQQx9SSqIoS9JBTH0Jakghr4kFcTQlwZsWB7vKIGhLw3UkD3eUfI2DNJCRMThwD3AX2XmnzfaTqW6j8qrqW709ZrMPL5pm7+jemTiqxa/YpXOM31pATLzJ1RPZfpgRJwSESuB/071zOK/Z3gfM6hC+eQsaYEy86aIuBr4W+DbVLd+3tZY3fXxjjmaD37RCPNMX+qPf091EnUucF5m7qu5HqktQ1/qj+cBz6H6/9Tapval9nhHjTiHd6QFioiDgU8Dn6V64MUnIuKEzHyM6tF2Z7dsMsqPd9SIc/aOtEAR8WHg9cDvAb8EbgUmM/OcxpTNe4Ergb8BXgH8FdWMnptrKlkFM/SlBYiIM4CvUD328H822tYCdwHvy8yPNfp8BHgB8DDwnzPzmsWvVjL0JakoXsiVpIIY+pJUEENfkgpi6EtSQQx9SSqIoS9JBTH0Jakghr4kFcTQl6SCGPqSVBBDX5IKYuhLUkH+Px4lRdkEW/MfAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "import matplotlib.pyplot as plt\n", + "axes = plt.figure(figsize=(4, 4), dpi=100).gca()\n", + "historyGD = onp.asarray(history)\n", + "updatesGD = onp.asarray(updates) # for later\n", + "axes.scatter(historyGD[:,0], historyGD[:,1], lw=0.5, color='blue')\n", + "axes.scatter([0], [0], lw=0.25, color='black') # target at 0,0\n", + "axes.set_xlabel('x0')\n", + "axes.set_ylabel('x1')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Newton\n", + "\n", + "For Newton's method, let's consider the combined function, i.e., \n", + "$L(\\mathbf{x}) = x_0^2 + x_1^4$\n", + "...\n", + "\n", + "Newton's method update step is given by\n", + "$\n", + "\\Delta \\mathbf{x} = \n", + "- \\eta \\left( \\frac{\\partial^2 L }{ \\partial \\mathbf{x}^2 } \\right)^{-1}\n", + " \\frac{\\partial L }{ \\partial \\mathbf{x} }\n", + "=\n", + "- \\eta \\ H_L^{-1} \\ ( J_{L} J_{\\mathbf{z}} )^T\n", + "$\n", + "\n", + "Hence we need to compute the Hessian and its inverse, easy in JAX, and very cheap here, but not so simple in practice...\n" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Newton iter 0: [2. 2.6666667]\n", + "Newton iter 1: [1.3333333 2.3703704]\n", + "Newton iter 2: [0.88888884 2.1069958 ]\n", + "Newton iter 3: [0.59259254 1.8728852 ]\n", + "Newton iter 4: [0.39506167 1.6647868 ]\n", + "Newton iter 5: [0.26337445 1.4798105 ]\n", + "Newton iter 6: [0.17558296 1.315387 ]\n", + "Newton iter 7: [0.1170553 1.1692328]\n", + "Newton iter 8: [0.07803687 1.0393181 ]\n", + "Newton iter 9: [0.05202458 0.92383826]\n" + ] + } + ], + "source": [ + "x = np.asarray([3.,3.])\n", + "eta = 1./3.\n", + "history = [x]; updates = []\n", + "\n", + "for i in range(10):\n", + " G = jax.grad(fun)(x)\n", + " H = jax.jacobian(jax.jacobian(fun))(x)\n", + " #H = jax.jacfwd(jax.jacrev(fun_Lz))(x) # more efficient? \n", + " Hinv = np.linalg.inv(H)\n", + " \n", + " x += -eta * np.matmul( Hinv , G)\n", + " history.append(x); updates.append( np.matmul( Hinv , G) )\n", + " print( \"Newton iter %d: \"%i + format(x) )\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let's compare, Newton actually does better, note - uses different step size..." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Text(0, 0.5, 'x1')" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAFtCAYAAADrr7rKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAg5ElEQVR4nO3df7TcdX3n8ef7XkKWDfdeLdKAAUwKBxcLDQqnNawHUJa2oglkt6wesW1gPRbXnKjr2jUL7bZdUvrjVKVrWqu2Qeoe6z8LJoqFWhp6KhEFJZU1QrGgEElQwHuvXYkX7nv/+M7V4Wbm3js3M/P9zv0+H+d8zzCf+Xxn3hkmr3zn+/3M5xOZiSSpHobKLkCS1D+GviTViKEvSTVi6EtSjRj6klQjhr4k1YihL0k1YuhLUo0cVXYB/RYRAbwYmCy7FknqohHg2znPL25rF/oUgf9Y2UVIUg+cBOyfq0MdQ38S4NFHH2V0dLTsWiTpiE1MTHDyySfDAs5g1DH0ARgdHTX0JdWOF3IlqUYMfUmqEUNfkmrE0JekGjH0JalGDH1JqpCc7u1qhqWGfkS8LSL+MSImGtueiHjtPPtcHhFfj4hnIuKrEXFJv+qVpF6YfHqS3X+0hUe3r+HAn53Mo9vXsPuPtjD5dPcnDogy18iNiPXAc8A/AQH8KvAe4OWZ+X9b9D8P+HtgK/Bp4E3AfwNekZn3L/A1R4Hx8fFxx+lLKt3k05McuGkdP/UT+xgenv5R+7PPDfHwU2dwwq/sYeSFI3M+x8TEBGNjYwBjmTkxV99Sj/Qzc1dm3pqZ/5SZD2bmNcD3gVe22eUdwF9n5h9m5r7M/A3gy8DmftUsSd10719cc1jgAxw1PM2an9jHvTuu7errVeacfkQMR8QbgRXAnjbd1gGfm9V2W6O93fMuj4jRmY1iUiJJqoRT/9WuwwJ/xlHD05x69M6uvl7poR8RZ0XE94FDwIeAjZn5tTbdTwAOzmo72GhvZysw3rQ52ZqkSsjp5KihqTn7HDU01dWLu6WHPvAAcDbwc8CfAh+LiJd18fmvB8aatpO6+NyStGgxFDw7vWzOPs/mMmIouvaapYd+Zv4wMx/KzHszcyuwl+LcfSsHgJWz2lY22ts9/6HMnJjZcB59SRXyjWfW8+xzraP42eeG+MahDV19vdJDv4UhYHmbx/YAF81qu5j21wAkqdLOuWobDz91xmHB/+xzQzz85Bmcc+V1XX29UqdWjojrgc8C36K4wPom4ELgFxqP3wTsb3wDALgBuDMi3g18BngjcC7w1v5WLkndMfLCEfiVPfzDjms59eidHDU0xbPTy/jGDzdwzpXXzTtcs1Nlj9P/c4oj9xMpLrL+I/D7mfk3jcd3A49k5qamfS4HrgNWU4zv//XMvLWD13ScvqTKyuns+Bx+J+P0Sw39Mhj6kpaagflxliSpvwx9SaoRQ1+SasTQl6QaMfQlqUYMfUmqEUNfkmrE0JekGjH0JalGDH1JqhFDX5JqxNCXpBox9CWpRgx9SaoRQ1+SasTQl6QaMfQlqUYMfUmqEUNfkmrE0JfUWs3Wz66Lo8ouQFKFTE3C3mtg/y6YnoKhZbBqPazdBstGyq5OXWDoSypMTcLt62B8HzD94/YHt8PBO+Dn9xj8S4CndyQV9l5zeOBDcX98H+y9toyq1GWGvqTC/l0cHvgzpmH/zn5Wox4x9CUVF22np+buMz3lxd0lwNCXBBHFRdu5DC0r+mmgGfqSCqvW0z4ShmDVhn5Wox4x9CUV1m6DsTM4PBaGiva115VRlbrM0JdUWDZSDMs8fTOsWA3HrCpuT9/scM0lJLJmF2YiYhQYHx8fZ3R0tOxypOrK9Bz+gJiYmGBsbAxgLDMn5urrkb6k1gz8JcnQl6QaMfQlqUYMfUmqEUNfkmqk1NCPiK0R8aWImIyIJyLiloh46Tz7bIqInLU906+aJWmQlX2kfwGwHXglcDGwDLg9IlbMs98EcGLT9pJeFilJS0Wp8+ln5i8234+ITcATwDnA38+9ax7oYWmStCSVfaQ/21jj9ql5+h0bEd+MiEcj4lMR8dPtOkbE8ogYndkAf1YoqbYqE/oRMQR8APh8Zt4/R9cHgKuAS4E3U/wZ7oqIk9r03wqMN22PdatmSRo0lZmGISL+FHgt8KrMXHAwR8QyYB/wicz8jRaPLweWNzWNAI85DYOkpaKTaRgqsUZuRHwQeD1wfieBD5CZUxHxFeC0No8fAg41vdaRlCpJA63sIZvRCPyNwGsy8+FFPMcwcBbweLfrk6Slpuwj/e3AmyjOz09GxAmN9vHM/AFARNwE7M/MrY37vwl8AXgIeAHwHoohmx/tb+mSNHjKDv23NW53z2q/Erix8d+n8PzVml8IfAQ4AXgauBc4LzO/1rMqJWmJqMyF3H5xPn1JS43z6UvdVLMDIy1tZZ/ekappahL2XgP7d8H0FAwtKxYOX7vNZQM10Ax9abapSbh9HYzv43mXkx7cDgfvcL1YDTRP70iz7b3m8MCH4v74Pth7bRlVSV1h6Euz7d/F4YE/Yxr27+xnNVJXGfpSs8ziHP5cpqe8uKuBZehLzSKKi7ZzGVpW9JMGkKEvzbZqPe3/agzBqg39rEbqKkNfmm3tNhg7g8P/egwV7WuvK6MqqSsMfWm2ZSPFsMzTN8OK1XDMquL29M0O19TAcxoGaT6ZnsNXpTkNg9RNBr6WEENfkmrE0JekGjH0JalGDH1JqhFDX5JqxNCXpBox9CWpRgx9SaoRQ1+SasTQl6QaMfQlqUYMfUmqEUNf5anZDK9SFRxVdgGqmalJ2HtNsfj49FSx9OCq9cXCJc5TL/Wcoa/+mZqE29fB+D5g+sftD26Hg3e4QInUB57eUf/svebwwIfi/vg+2HttGVVJtWLoq3/27+LwwJ8xDft39rMaqZYMffVHZnEOfy7TU17clXrM0Fd/RBQXbecytMylCaUeM/TVP6vW0/4jNwSrNvSzGqmWDH31z9ptMHYGh3/shor2tdeVUZVUK4a++mfZSDEs8/TNsGI1HLOquD19s8M1pT6JrNmFs4gYBcbHx8cZHR0tu5x6y/QcvtQFExMTjI2NAYxl5sRcfUs90o+IrRHxpYiYjIgnIuKWiHjpAva7PCK+HhHPRMRXI+KSftSrLjPwpb4r+/TOBcB24JXAxcAy4PaIWNFuh4g4D/gE8OfAy4FbgFsi4syeVytJA65Sp3ci4njgCeCCzPz7Nn0+CazIzNc3tX0BuC8zr17Aa3h6R9KSMjCnd1oYa9w+NUefdcDnZrXd1mg/TEQsj4jRmQ3waqGk2qpM6EfEEPAB4POZef8cXU8ADs5qO9hob2UrMN60PXZklUrS4KpM6FOc2z8TeGOXn/d6im8QM9tJXX5+SRoYlZhaOSI+CLweOD8z5zsSPwCsnNW2stF+mMw8BBxqeq0jqFSSuqeMUctlD9mMRuBvBF6TmQ8vYLc9wEWz2i5utEtSpU1OwpYtsGYNnHxycbtlS9HeD6WO3omIPwHeBFwKPND00Hhm/qDR5yZgf2Zubdw/D7gTeC/wGYrTQf8deMU81wJmXtPRO5JKMTkJ69bBvn0w3TTL+NAQnHEG7NkDI4sYajJIo3feRnGefTfweNP2hqY+pwAnztzJzLso/qF4K7AX+CXgsoUEviSV6ZprDg98KO7v2wfX9mEdoUqN0+8Hj/QllWXNGnjkkfaPr14NDy/kJPcsg3Skryqp2QGA1E+ZMDXPOkJTfVhHqBKjd1Siqcli7dr9u4qVq4aWFfPer93mrJdSF0XAsnnWEVrWh3WEPNKvs6lJuH0dPLgd/uUR+MH+4vbB7UX7VJ+GE0g1sX59cdG2laEh2NCHdYQM/Trbew2M7+Pwxcqni/a9fbiqJNXItm3FKJ3ZwT8zeue6PqwjZOjX2f5dHB74M6Zh/85+ViMteSMjxbDMzZuLi7arVhW3mzcvfrhmpzynX1eZxTn8uUxPudCJ1GUjI3DDDcVWu1/kqkQRxUXbuQz14aqSVGNl/PUy9Ots1XrafwSGYFUfripJ6itDv87WboOxMzj8YzBUtK/tw1UlSX1l6NfZshH4+T1w+mZYsRqOWVXcnr65aHecvrTkOA2DfsyLttJAchoGLY6BLy15hr4k1YihL0k1YuhLUo0Y+pJUI4a+JNWIoS9JNWLoS1KNGPqDrmY/rpN0ZJxaeRC5xKGkRTL0B83MEoezV7x6cDscvMM5c6QKqtIMJ57eGTQucSgNhMlJ2LIF1qyBk08ubrdsKdrL5IRrg+ZTa4rFy9tZsRoufbhf1UhqYXIS1q2Dfftguun4bGYt3G4vjeiEa0tVJ0scSirNNdccHvhQ3N+3D64t8Qu5oT9IXOJQGgi7dh0e+DOmp2Hnzv7W08zQHzQucShVWiZMzfOFfKrEL+SG/qBxiUOp0iJg2TxfyJeV+IXc0B80LnEoVd769cVF21aGhmBDiV/IHb0z6Ko0AFgS4Ogd9ZKBL1XOyEgR7Js3w+rVsGpVcbt5c/cDv1Me6UtSj/X6C7lH+pJUIVX6Qm7oS1KNdC30I+KMiPjnDvc5PyJ2RcS3IyIj4rJ5+l/Y6Dd7O+GIipekmujmkf7RwEs63GcFsBd4e4f7vRQ4sWl7osP9JamWFjy1ckS8b54ux3f64pn5WeCzjefvZNcnMvN7nb7eQHAIpqQe6mQ+/XcA9wHtrgwfe8TVLNx9EbEcuB/4rcz8fB9fu/tcFEVaMqp+3NZJ6D8EvD8zP97qwYg4G7i3G0XN4XHgauAeYDnwFmB3RPxcZn65TV3LG31nVCtFXRRFGniTk8XMmrt2FfPqLFtW/Cp327Zyx+S30kno3wOcA7QMfSCBnv77lpkPAA80Nd0VEacC7wJ+uc1uW4H/0cu6jshCFkU594YyKpO0AO1+fbt9O9xxR/k/xpqtkwu57wY+0O7BzNybmWUMAf0icNocj18PjDVtJ/WjqAXbv4vDA3/GNOwvcQ5WSfOq8tz5rSw4pDPzQGZ+MyJe3a5PRPxad8rqyNkUp31aysxDmTkxswElL1bWxEVRpIFX5bnzW1nMkflfR8QfRsSPJg+NiBdFxC7g9zp5oog4NiLOblwPAFjTuH9K4/HrI+Kmpv7vjIhLI+K0iDgzIj4AvAbYvog/R/lcFEUaaFWfO7+VxYT+q4GNwJci4mUR8TqKUTSjFEfdnTgX+EpjA3hf479/p3H/ROCUpv5HA38EfBW4E1gL/LvM/NvO/xgV4aIo0sCq+tz5rXQc+pl5F0W43w98GbgZeD9wYWZ+s8Pn2p2Z0WLb1Hh8U2Ze2NT/DzLztMw8JjOPy8xXZ+bfdfpnqBQXRZEGWpXnzm9lsRdeT6c4Sn8MeJbiF7L/ultF1YqLokgDbdu2Yo782cE/M3f+dRU7but4auWIeC/w28CHgfdQjJz5S4rTO2/OzD3dLrKbKj+1ctV/2SHpMJOTxSidnTt/PE5/w4Yi8PsxXLOTqZUXE/qPA1c1plCYaVsG/C6wJTOXt925Aiof+pIGWhnHbZ2Efic/zppxVmZ+t7khM6eA90TEpxfxfJK0ZFT9i/piLuR+d47H7jyyciRpMFRpGGYnXERFkhZochK2bIE1a+Dkk4vbLVuK9kHhGrn95oVaaSC1m2NnZpROmXPsuEZu1UxNwj1b4FNr4JaTi9t7thTtkgbCoM2x045H+r3WburkmR9fORZfGghr1sAjj7R/fPVqePjhflXzfB7pV8lCpk6WVGmDOMdOO4Z+rzl1sjTwBnGOnXYM/V5y6mRpyRi0OXbaMfR7yamTpSUhc/Dm2GnH0O81p06WBtLsMfk/8zPwqlfBW99aXLRdtaq43by5eksizsXRO73m6B1p4CxkTP6xx1bnS7qjd6rEqZOlgbOQMflVCfxOeaTfb/4iV6q8Ko/Jb8Uj/Soz8KVKW0pj8lsx9CWpyVIak9+Kod8rg3oYINVc5tIZk9/KYhZRUTtTk8W0C/t3FT+6GlpWDNlcu80LtlKFTU4WF2937SpO3QwPwwteAE8//fzjt0Ebk9+Kod8t7YZmPrgdDt7hSB2pouYanvnCFxZDM597rv/r3vaKp3e6xYnVpIE01/DM730PLrsMHn20GK1zww2DHfhg6HePE6tJA2nXrsMDf8b0NOzcObgXbVsx9LvBidWkgbTUh2e2Yuh3gxOrSQNrKQ/PbMXQ7xYnVpMGwuyJ1J58sn3fQR+e2YrTMHSLE6tJlddupE4rVVjwfKGchqEMTqwmVV67kTpQnMIZGRncKZMXyiP9XnFiNalyFjKR2j//8+D91fVIvwoG7VMjLXHT0wsbqbPUGfrdVLNvTVLVNV+0PeUUOHBg7v5LbaROK07DcKScb0eqpE4u2sLSHKnTiqF/JJxvR6qsuS7azrYUJlJbKE/vHAnn25EqKXPu6RUAjjpq6Y/UaaXU0I+I8yNiV0R8OyIyIi5bwD4XRsSXI+JQRDwUEZt6X2kbzrcjVcbsH109+ujc/VeuhG99a+lMpLZQZZ/eWQHsBf4C+D/zdY6INcBngA8BVwAXAR+NiMcz87ZeFnqYTubbWepXhqSSdXr+HoqLtu0WSlnKSg39zPws8FmAWFgwXg08nJnvbtzfFxGvAt4F9Df0nW9HqoxOzt9DfS7atjJo/86tAz43q+22RntLEbE8IkZnNqB7X+Kcb0eqhJ07Owv8uly0bWXQQv8E4OCstoPAaEQc02afrcB40/ZY16pZu62YV+ewt7Ex387amn6qpD5oPoc/3/n74WF48Yvrd9G2lbLP6ffD9cD7mu6P0K3gn5lvZ++1xUXbH43T31AEvsM1pZ7o9Bz+yScP5vQKvTBooX8AWDmrbSUwkZk/aLVDZh4CDs3cX+C1g4U76lg494Zi86Kt1BedjsHfsMG/mjMGLfT3AJfMaru40d4//gpXKtV8Y/Bn1P38fSulhn5EHAuc1tS0JiLOBp7KzG9FxPXAqsz8lcbjHwI2R8QfUAzzfA3wH4HX9a1of4UrlWohSxwODRWndC69tAj8up6/b6XsI/1zgb9ruj9z7v1jwCbgROCUmQcz8+GIeB3wfuAdFOfm39LXMfoL+RXuuTf0rRypbiLmX+LwlFOKH13pcKWO3snM3ZkZLbZNjcc3ZeaFLfZ5eWYuz8xTM/PGvhbtr3Cl0q1f3/6HVXUeg78QgzZks1yd/ApXUs9s21acq58d/J7Dn5+h3wl/hStVwshIMdZ+8+Zi7H0dJ05brLLP6Q+eVeuLi7YtT/H4K1ypX0ZGionSbnC0dEc80u+Uv8KVKsfAXzhDv1Mzv8I9fTOsWA3HrCpuT9/scE1JlRdZs4uOjUnXxsfHxxkdHT3yJ/R7paSSTUxMMDY2BjCWmRNz9fVI/0gZ+JIGiKEvSTVi6EtSjRj6klQjhr4k1YihL0k1YuhLUo0Y+pJUI4a+JNWIoS9JNWLoS1KNGPqSVCOGviTViKEvSTVi6EtSjRj6klQjhr4k1YihL0k1YuhLUo0Y+pJUI4a+JNWIoS9JNWLoS1KNGPqSVCOGviTViKEvSTVi6EtSjRj6klQjhr4k1UglQj8i3h4Rj0TEMxFxd0T87Bx9N0VEztqe6We9kjSoSg/9iHgD8D7gt4FXAHuB2yLiJ+fYbQI4sWl7Sa/rlKSloPTQB/4L8JHM3JGZXwOuBv4fcNUc+2RmHmjaDvalUkkacKWGfkQcDZwDfG6mLTOnG/fXzbHrsRHxzYh4NCI+FRE/PcdrLI+I0ZkNGOlW/ZI0aMo+0n8RMAzMPlI/CJzQZp8HKL4FXAq8meLPcFdEnNSm/1ZgvGl77AhrlqSBVXbodywz92TmTZl5X2beCfx74DvAr7XZ5XpgrGlr94+DJC15R5X8+t8FngNWzmpfCRxYyBNk5lREfAU4rc3jh4BDM/cjYnGVStISUOqRfmb+ELgXuGimLSKGGvf3LOQ5ImIYOAt4vBc1StJSUvaRPhTDNT8WEfcAXwTeCawAdgBExE3A/szc2rj/m8AXgIeAFwDvoRiy+dF+Fy5Jg6b00M/MT0bE8cDvUFy8vQ/4xaZhmKcA0027vBD4SKPv0xTfFM5rDPeUJM0hMrPsGvqqMWxzfHx8nNHR0bLLkaQjNjExwdjYGMBYZk7M1XfgRu9IkhbP0JekGjH0JalGDH1JqhFDX5JqxNCXpBox9CWpRgx9SaoRQ1+SasTQl6QaMfQlqUYMfUmqEUNfkmrE0JekGjH0JalGDH1JqhFDX5JqxNCXpBox9CWpRgx9SaoRQ1+SasTQl6QaMfQlqUYMfUmqEUNfkmrE0JekGjH0JalGDH1JqhFDX5JqxNCXpBox9CWpRgx9SaoRQ1+SasTQl6QaMfQlqUYqEfoR8faIeCQinomIuyPiZ+fpf3lEfL3R/6sRcUmvastMduzYwYYNG1i3bh0bNmxgx44dZGavXlKSeuaosguIiDcA7wOuBu4G3gncFhEvzcwnWvQ/D/gEsBX4NPAm4JaIeEVm3t/N2jKTTZs28fGPf5zp6ekftd96663s3r2bG2+8kYjo5ktKUk9F2UesEXE38KXM3Ny4PwQ8CvyvzPy9Fv0/CazIzNc3tX0BuC8zr17A640C4+Pj44yOjs7Zd8eOHbzlLW95XuDPGB4e5sMf/jBXXXXVfC8pST01MTHB2NgYwFhmTszVt9TTOxFxNHAO8LmZtsycbtxf12a3dc39G25r1z8ilkfE6MwGjCy0vptvvrll4AM899xz3HzzzQt9KkmqhLLP6b8IGAYOzmo/CJzQZp8TOuy/FRhv2h5baHHf+c535nz8ySefXOhTSVIllB36/XA9MNa0nbTQHY8//vg5Hz/uuOOOqDBJ6reyQ/+7wHPAylntK4EDbfY50En/zDyUmRMzGzC50OI2btzI8PBwy8eGh4fZuHHjQp9Kkiqh1NDPzB8C9wIXzbQ1LuReBOxps9ue5v4NF8/Rf9E2bdrEFVdccVjwDw8Pc8UVV3DllVd2+yUlqadKH7JJMVzzYxFxD/BFiiGbK4AdABFxE7A/M7c2+t8A3BkR7wY+A7wROBd4a7cLiwhuvPFGLrjgAm6++WaefPJJjjvuODZu3MiVV17pcE1JA6f0IZsAEbEZeA/Fxdj7gC2ZeXfjsd3AI5m5qan/5cB1wGrgn4Bfz8xbF/haCx6yKUmDoJMhm5UI/X4y9CUtNQMzTl+S1F+GviTViKEvSTVi6EtSjRj6klQjhr4k1UgVfpxViomJOUc1SdLA6CTP6jhOfxUdzLQpSQPkpMzcP1eHOoZ+AC+mg4nXGkYo/rE4aRH7ls3ay2Ht5ahr7SPAt3OeUK/d6Z3GGzLnv4StNM2zMznfL96qxtrLYe3lqHHtC+rvhVxJqhFDX5JqxNBfuEPAbzduB421l8Pay2Htc6jdhVxJqjOP9CWpRgx9SaoRQ1+SasTQl6QaMfSbRMTbI+KRiHgmIu6OiJ+dp//lEfH1Rv+vRsQl/aq1RS0Lrj0iNkVEztqe6We9TbWcHxG7IuLbjTouW8A+F0bElyPiUEQ8FBGbel9pyzo6qr1R9+z3PSPihD6VPFPH1oj4UkRMRsQTEXFLRLx0AfuV/nlfTO1V+bxHxNsi4h8jYqKx7YmI186zT9ffc0O/ISLeALyPYrjUK4C9wG0R8ZNt+p8HfAL4c+DlwC3ALRFxZl8Kfn4tHdXeMAGc2LS9pNd1trGCot63L6RzRKwBPgP8HXA28AHgoxHxCz2qby4d1d7kpTz/vX+iy3XN5wJgO/BK4GJgGXB7RKxot0OFPu8d195Qhc/7Y8B7gXOAc4E7gE9FxE+36tyz9zwz3Yphq3cDH2y6P0QxXcN72/T/JPDpWW1fAD40ALVvAr5X9nveoq4ELpunz+8D989q+yvgrweg9gsb/V5Q9ns9q67jG3WdP0efynzeF1F7JT/vjdqeAv5TP99zj/SBiDia4l/fz820ZeZ04/66Nruta+7fcNsc/XtikbUDHBsR34yIRyOi7dFGBVXifT9C90XE4xHxNxHxb8suBhhr3D41R5+qvu8LqR0q9nmPiOGIeCPFt8U9bbr15D039AsvAoaBg7PaDwLtzree0GH/XllM7Q8AVwGXAm+m+BzcFREn9arILmr3vo9GxDEl1NOJx4Grgf/Q2B4FdkfEK8oqKCKGKE6RfT4z75+ja1U+7z/SQe2V+bxHxFkR8X2KX9x+CNiYmV9r070n73ntZtkUZOYemo4uIuIuYB/wa8BvlFXXUpeZD1AE0Iy7IuJU4F3AL5dTFduBM4FXlfT6R2JBtVfs8/4AxbWoMeCXgI9FxAVzBH/XeaRf+C7wHLByVvtK4ECbfQ502L9XFlP782TmFPAV4LTultYT7d73icz8QQn1HKkvUtL7HhEfBF4PvDoz51tYqCqfd6Dj2p+nzM97Zv4wMx/KzHszcyvFQIB3tOnek/fc0Kf4HwHcC1w009b46ngR7c+37Wnu33DxHP17YpG1P09EDANnUZx+qLpKvO9ddDZ9ft+j8EFgI/CazHx4AbtV4n1fZO2zn6NKn/chYHmbx3rznpd99boqG/AG4BngV4EzgD8DngZWNh6/Cbi+qf95wBTwbuDfAL8F/BA4cwBq/03g54Gfohji+QngB8DLSqj9WIrgO5tiFMa7Gv99SuPx64GbmvqvAf4F+IPG+/6fgWeBXxiA2t9JcV75NIrTEh+g+JZ2UZ/r/hPgexTDH09o2o5p6lPJz/sia6/E573xeTgfWE3xj871wDRwcT/f877+Jan6BmwGvklxkeVu4OeaHtsN3Dir/+UU5+gOAfcDlwxC7cD7m/oeoBj3/vKS6r6wEZiztxsbj98I7G6xz1ca9X8D2DQItQO/DjzUCJwnKX5r8OoS6m5Vcza/j1X9vC+m9qp83inG2z/SqOMJipE5F/f7PXdqZUmqEc/pS1KNGPqSVCOGviTViKEvSTVi6EtSjRj6klQjhr4k1YihL0k1YuhLPVaV5R0lMPSlnqrY8o6S0zBIRyIijge+CvxxZv5uo+08inlUXksx0dfrMvPMpn3+imLJxF/sf8WqO4/0pSOQmd+hWJXptyLi3IgYAf6SYs3iv6W6ywyqplw5SzpCmXlrRHwE+N/APRRTP29tPDzn8o45mAu/aIB5pC91x3+lOIi6HLgiMw+VXI/UkqEvdcepwIsp/k6tbmpfass7asB5ekc6QhFxNPBx4JMUC158NCLOyswnKJa2u2TWLoO8vKMGnKN3pCMUEX8I/BKwFvg+cCcwnpmvbwzZvB/YDvwF8BrgjylG9NxWUsmqMUNfOgIRcSHwNxTLHv5Do201sBd4b2b+aaPP+4GXAY8B/zMzb+x/tZKhL0m14oVcSaoRQ1+SasTQl6QaMfQlqUYMfUmqEUNfkmrE0JekGjH0JalGDH1JqhFDX5JqxNCXpBox9CWpRv4/lCpJxjSNE4EAAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "axes = plt.figure(figsize=(4, 4), dpi=100).gca()\n", + "historyNt = onp.asarray(history)\n", + "updatesNt = onp.asarray(updates) \n", + "axes.scatter(historyGD[:,0], historyGD[:,1], lw=0.5, color='blue')\n", + "axes.scatter(historyNt[:,0], historyNt[:,1], lw=0.5, color='orange')\n", + "axes.scatter([0], [0], lw=0.25, color='black') # target at 0,0\n", + "axes.set_xlabel('x0')\n", + "axes.set_ylabel('x1')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Physical Gradients\n", + "\n", + "Now we also use the _inverse physics_, i.e. the inverse of z:\n", + "$\\mathbf{z}^{-1}(\\mathbf{x}) = [x_0 \\ x_1^{1/2}]^T$\n", + "\n", + "Let’s call intermediate physical result $\\mathbf{y}$, i.e. \n", + "$\\mathbf{z}(\\mathbf{x}) = \\mathbf{y}$\n", + "\n", + "With an update step:\n", + "$\n", + "\\Delta \\mathbf{x} = \n", + "\\mathbf{z}^{-1} \\left( \\mathbf{y} - \\eta_L \n", + " \\left( \\frac{\\partial^2 L }{ \\partial \\mathbf{z}^2 } \\right)^{-1}\n", + " \\frac{\\partial L }{ \\partial \\mathbf{z} }\n", + "\\right) - \\mathbf{x}\n", + "$\n", + "\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "PG iter 0: [2.1 2.5099802]\n", + "PG iter 1: [1.4699999 2.1000001]\n", + "PG iter 2: [1.0289999 1.7569861]\n", + "PG iter 3: [0.72029996 1.47 ]\n", + "PG iter 4: [0.50421 1.2298902]\n", + "PG iter 5: [0.352947 1.029 ]\n", + "PG iter 6: [0.24706289 0.86092323]\n", + "PG iter 7: [0.17294402 0.7203 ]\n", + "PG iter 8: [0.12106082 0.60264623]\n", + "PG iter 9: [0.08474258 0.50421 ]\n" + ] + } + ], + "source": [ + "x = np.asarray([3.,3.])\n", + "eta = 0.3\n", + "history = [x]; updates = []\n", + "\n", + "def fun_z_inv_analytic(y):\n", + " return np.array( [y[0], np.power(y[1],0.5)] )\n", + "\n", + "for i in range(10):\n", + " y = fun_z(x)\n", + " \n", + " # Newton step for L(y)\n", + " GL = jax.grad(fun_L)(y)\n", + " HL = jax.jacobian(jax.jacobian(fun_L))(y)\n", + " HLinv = np.linalg.inv(HL)\n", + " y += -eta * np.matmul( HLinv , GL)\n", + "\n", + " # \"inverse physics\" step via z-inverse\n", + " x = fun_z_inv_analytic(y)\n", + " history.append(x)\n", + " updates.append( history[-2] - history[-1] )\n", + " print( \"PG iter %d: \"%i + format(x) )\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "All together now" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Text(0, 0.5, 'x1')" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAFtCAYAAADrr7rKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAjEElEQVR4nO3df5RcZZ3n8fe3Ok0Om3Q3I2ICncQEOcEg2CgcNVkOIAz+QBKIK6tHnJ3AuopLn6jrOks2jMcZiZlxziKZpdXxxyRmcBj+IqQVB2QY2F0To6AJohFkTSRpSECQ7taF0KS++8dThZXqquqq7qq6z637eZ1zT1G3nlv1pKj+1K17n/t9zN0REZFsyCXdARERaR+FvohIhij0RUQyRKEvIpIhCn0RkQxR6IuIZIhCX0QkQxT6IiIZMivpDrSbmRlwCjCedF9ERJqoB3jSp7jiNnOhTwj8g0l3QkSkBRYAI7UaZDH0xwEOHDhAb29v0n0REZmxsbExFi5cCHUcwchi6APQ29ur0BeRzNGJXBGRDFHoi4hkiEJfRCRDFPoiIhmi0BcRyRCFvohIRPxovqXPn2jom9nHzOxhMxsrLDvN7N1TbHOlmf3CzF40s5+a2aXt6q+ISCuMP/EkB88b4OUTZ5F/dTcvnziLg+cNMP7Ek01/raT39A8C1wPnAOcC9wF3mtkbKjU2sxXAbcA3gDcB24BtZnZmW3orItJk4088yXFvPZX+HQ8z67mjdD2fZ9ZzR+nf+TDHvfXUpge/xTYxupk9B3za3b9R4bHbgTnuflnJuh8Au9392jqfvxcYHR0d1cVZIpK4g+cN0L/jYaxCFHsORlYMsOB/7675HGNjY/T19QH0uftYrbZJ7+m/wsy6zOwDwBxgZ5Vmy4F7y9bdXVhf7Xlnm1lvcSEUJRIRicL8vT+rGPgAlof5P/9ZU18v8dA3s7PM7HfAEeArwGp3/3mV5vOBw2XrDhfWV7MOGC1ZVGxNRKLgR/NYvvbRFsvnm3pyN/HQBx4FzgbeCnwZ+KaZndHE598I9JUsC5r43CIi02ZdOTxnNdt4zrCu5kV14qHv7i+5++Pu/pC7rwP2AB+v0vwQMK9s3bzC+mrPf8Tdx4oLqqMvIhE5tOwNeJUk9hwcOqO541QSD/0KcsDsKo/tBC4uW3cJ1c8BiIhEre8fv8tLr5k9Kfg9By+9ZjZ937qrqa+X9Dj9jWZ2vpktLhzb3whcCHyr8PjWwrqiTcC7zOxTZvZ6M/ssYajnLe3uu4hIM/QsOoWXdv2KkRUDvPyqWRw9IcfLr5rFyIoBXtr1K3oWndLU10u6nv5rgK3AyYSTrA8D73T37xUeXwS8cgbD3XeY2QeBG4HPA78ErnD3R9raaxGRJupZdAo9hWGZfjSPdeVadvIxunH6raZx+iLSaVI5Tl9ERFpPoS8ikiEKfRGRDFHoi4hkiEJfRCRDFPoiIhmi0BcRyRCFvohIhij0RUQyRKEvIpIhCn0RkQxR6IuIZIhCX0QkQxT6IiIZotAXEckQhb6ISIYo9EVEMkShLyKSIQp9EZEMUeiLSGUZmz87K2Yl3QERicjEOOxZDyPDkJ+AXDf0r4SBDdDdk3TvpAkU+iISTIzDPcthdC+Q/8P6x4bg8H3wjp0K/g6gwzsiEuxZPznwIdwf3Qt7bkiiV9JkCn0RCUaGmRz4RXkY2d7O3kiLKPRFJJy0zU/UbpOf0MndDqDQFxEwCydta8l1h3aSagp9EQn6V1I9EnLQv6qdvZEWUeiLSDCwAfqWMTkWcmH9wI1J9EqaTKEvIkF3TxiWuXQQ5iyG4/vD7dJBDdfsIOYZOzFjZr3A6OjoKL29vUl3RyRe7jqGnxJjY2P09fUB9Ln7WK222tMXkcoU+B1JoS8ikiEKfRGRDFHoi8jMZOy8YNop9EWkcePjsHYtLFkCCxeG27Vrw3qJWqKhb2brzOxHZjZuZk+b2TYzO32KbdaYmZctL7arzyKZNz4Oy5fD0BDs3w8jI+F2aCisV/BHLek9/QuAIeBtwCVAN3CPmc2ZYrsx4OSS5bWt7KSIlFi/HvbuhXxZcbZ8Pqy/QdU4Y5Zo6Lv7u9x9i7v/zN33AGuARcA5U2/qh0qWwy3vrIgEw8OTA78on4ftqsYZs6T39Mv1FW6fm6LdXDP7tZkdMLM7zewN1Rqa2Wwz6y0ugC4rFJkud5iYohrnhKpxxiya0DezHHAz8H13f6RG00eBa4DLgQ8R/g07zGxBlfbrgNGS5WCz+iySOWbQPUU1zm5V44xZNKFPOLZ/JvCBWo3cfae7b3X33e7+APBe4Bngo1U22Uj4BVFcqn05iEg9Vq6EXJXoyOVglapxxiyK0DezW4DLgLe7e0N74u4+AfwEOK3K40fcfay4ABpaIDITGzbAsmWTgz+XC+tvVDXOmCU9ZNMKgb8auMjd903jObqAs4Cnmt0/Eamgpwd27oTBQVi8GPr7w+3gYFjfo9NmMUu0yqaZfQn4IOH4/KMlD426+wuFNluBEXdfV7j/GeAHwOPACcCngSuAc9z953W8pqpsijSTqnEmrpEqm7Pa06WqPla4vb9s/dXAlsJ/L+LY2Zr/CPgaMB/4LfAQsKKewBeRFlDgp4rq6YuIpJzq6Ys0U8Z2jKSzJX14RyROE+OwZz2MDEN+AnLdYeLwgQ2aNlBSTaEvUm5iHO5ZDqN7OeZ00mNDcPg+zRcrqabDOyLl9qyfHPgQ7o/uhT0qKCbppdAXKTcyzOTAL8rDiAqKSXop9EVKuYdj+LXkVVBM0kuhL1LKLJy0rSWngmKSXgp9kXL9K6n+p5GDfhUUk/RS6IuUG9gAfcuY/OeRC+sHVFBM0kuhL1KuuycMy1w6CHMWw/H94XbpoIZrSuqpDIPIVFRQTCKnMgwizaTAlw6i0BcRyRCFvohIhij0RUQyRKEvIpIhCn0RkQxR6IuIZIhCX0QkQxT6IiIZotAXEckQhb6ISIYo9EVEMkShLyKSIQp9SU7GKryKxGBW0h2QjJkYhz3rw+Tj+Ykw9WD/yjBxSafUqVcpZomY9vSlfSbG4Z7l8NgQ/H4/vDASbh8bCusnxpPu4fSNj8PatbBkCSxcGG7Xrg3rRSKiSVSkfR5cGwKefIUHc2FmqnM3tbtXMzc+DsuXw969kC/5t+VysGwZ7NwJPR3yK0aipElUJE4jw1QOfML6ke3t7E3zrF8/OfAh3N+7F264IZl+iVSg0Jf2cA/H8GvJT6Tz5O7w8OTAL8rnYXtKv8ykIyn0pT3MwknbWnLd6TsB6g4TU3yZTaT0y0w6kkJf2qd/JdU/cjnoX9XO3jSHGXRP8WXWncIvM+lYCn1pn4EN0LeMyR+7XFg/cGMSvZq5lSvDSdtKcjlYlcIvM+lYGr0j7TUxDntuCCdtXxmnvyoEflrH6Wv0jiSskdE7Cn1JTiddxDQ+HkbpbN8ejuF3d4c9/BtvVOBLy6Um9M1sHfBe4PXAC8AO4L+5+6NTbHcl8DlgMfDLwjZ31fmaCn1prU76MpNUSNM4/QuAIeBtwCVAN3CPmc2ptoGZrQBuA74BvAnYBmwzszNb3luReijwJWJRHd4xs5OAp4EL3P1/VWlzOzDH3S8rWfcDYLe7X1vHa2hPX0Q6Spr29Mv1FW6fq9FmOXBv2bq7C+snMbPZZtZbXAAdYBWRzIom9M0sB9wMfN/dH6nRdD5wuGzd4cL6StYBoyXLwZn1VEQkvaIJfcKx/TOBDzT5eTcSfkEUlwVNfn4RkdSIop6+md0CXAac7+5T7YkfAuaVrZtXWD+Jux8BjpS81gx6KiLSPEkM9Ep0T9+CW4DVwEXuvq+OzXYCF5etu6SwXkQkaklPvZD0OP0vAR8ELgdKx+aPuvsLhTZbgRF3X1e4vwJ4ALge+A7hcNB/B948xbmA4mtq9I6IJKJVF2+nafTOxwjH2e8HnipZ3l/SZhFwcvGOu+8gfFF8BNgDvA+4op7AFxFJUgxTL0Q1Tr8dtKcvIklZsgT276/++OLFsK+eg9xl0rSnLzHJ2A6ASDvFMvVCFKN3JEET47BnfZjK8JWqlytDGeS0Vr0UiVAsUy9oTz/LJsbhnuVhsvLf74cXRsLtY0Nh/USbhhOIZEQMUy8o9LNsz3oY3cvkycrzYf0eTegt0kwbNoRROuXBXxy9c2Mb5hFS6GfZyDCTA78oHyY6EZGm6ekJwzIHB8NJ2/7+cDs42L65dnRMP6vcwzH8WvITqg0v0mQ9PbBpU1gyd0WuJMgsnLStJacJvUVaKYk/L4V+lvWvpPpHIBfmrhWRjqLQz7KBDdC3jMkfg1xYP9CGs0oi0lYK/Szr7oF37ISlgzBnMRzfH26XDob1Gqcv0nFUhkH+IA0nbdPQR5E2UxkGmZ5YwzTpWrQiHUR7+hK3VtWiFekg2tOXzhFDLVqRDqLQl7gND08O/KJ8HrbrqmGRRij0JV6x1KIV6SAKfYlXLLVoRTqIQl/iFkMtWpEOotCXuMVQi1akgyj0067Tj2fHUItWpINonH4aZXmKQ12RKzJJI+P0VU8/bYpTHJbPePXYEBy+r/Nr5ijwJYVi2lfR4Z200RSHIqkQa/UQHd5JmzuXhMnLq5mzGC7f167eiEgF7a4eojIMnaqRKQ5FJDExVw9R6KeJpjgUSYWYq4co9NNGUxyKRC326iEK/bTRFIciUYu9eohCP200xaFI9GKuHqLRO2kX0wBgEQE0ekdaSYEvEp2Yq4doT19EpMVa/YNce/oiIhGJ6Qe5Ql9EJEOaFvpmtszMftXgNueb2bCZPWlmbmZXTNH+wkK78mX+jDovIpIRzdzTPw54bYPbzAH2ANc1uN3pwMkly9MNbi8ikkl1l1Y2s5umaHJSoy/u7t8Fvlt4/kY2fdrdn2/09VIhTUMw09RXEQEaq6f/cWA3UO3M8NwZ96Z+u81sNvAI8Fl3/34bX7v50jQpyvh4qCY1PByuJe/uDleibNigWaxEiH9fqO4hm2b2KPA5d7+1yuNnAw+5e9e0OmLmwGp331ajzenAhcCDwGzgw8CfAG919x9X2WZ2oW1RD3AwmiGb1SZFKZZViOkq23ZfcSKSEknvC7VqyOaDwDk1Hnegpd9v7v6ou/+duz/k7jvc/RpgB/DJGputA0ZLloOt7GPD0jQpSsz1YkUSUtwXGhqC/fthZCTcDg2F9UlPmlKukdD/FHBztQfdfY+7JzEE9IfAaTUe3wj0lSwL2tGpuo0MMznwi/IwkmAN1nIx14sVSUja9oXqDml3P+Tuvzazt1drY2YfbU63GnI28FS1B939iLuPFRcgnu/dNE2KEnu9WJGEpG1faDp75v9sZn9jZq8UDzWzV5vZMPBXjTyRmc01s7ML5wMAlhTuLyo8vtHMtpa0/4SZXW5mp5nZmWZ2M3ARMDSNf0fy0jQpSuz1YkUSkMZ9oemE/tuB1cCPzOwMM3sPYRRNL2GvuxHnAj8pLAA3Ff77Lwv3TwYWlbQ/DvgfwE+BB4AB4I/d/V8a/2dEIk2TosRcL1YkAWncF5pWwTUzmwt8BXgfIbH+HPiCp6B6W3QF1zR6RyTV1q4NJ20rHeLJ5UJlzU2bWtuHdhRcW0rYSz8IvEy4QvbfTPO5si1Nk6LEXC9WJCEbNoR9nvIfwcV9oRsjm8yu4T19M7se+Avgq8CnCSNn/oFweOdD7r6z2Z1spuj29MvFfmVHqTT1VaSFxsfDKJ3t2/8wTn/VqhD4sY3Tn07oPwVcUyihUFzXDXweWOvus6tuHIHoQ19EUi2JfaFGQr+RMgxFZ7n7b0pXuPsE8Gkz+/Y0nk9EpGPE/uO34WP65YFf9tgDM+uOiEg6xD9spTJNoiIiUqfx8TBaZ8kSWLgw3K5dG1+phVo0R2676eSnSCrFPGJZc+TGZmIcHlwLdy6BbQvD7YNrw3oRSYW01dipRnv6rZami69EpKolS0L1zGoWL4Z9+9rVm2NpTz8maSqdLCIVpbHGTjUK/VZLU+lkEakojTV2qlHot1KaSieLSE2dUm9Qod9KaSqdLCJVuaevxk41Cv1Wi6l0sn5RiNStfEz+G98I550HH/lIuusNavROqyU9eifpGZtFUqieMflz58bzI72lBdfSLpGLsybGwyidke3hGH6uO+zhD9zY+sCP9WoSkYjFUCO/EQr9GjJ1RW7aPrkikYh5TH4lGqcfs3b+HkzbjM0iEeikMfmVKPQ7Vad/ckVapJPG5Fei0G+VpMO00z+5Ii3i3jlj8itR6DdTbIXVOvmTK9JE5cMz77wTTjhh8j5R2sbkV6ITuc2S9NDMSjR6R2RKtf5MTjghDM08erT98942QidykxBjYbWenhDsg4PpvppEpIVqlUx+/nm44go4cCCM1tm0Kf1/NtrTb5Y7l8Dv91d/fM5iuDzhMV6awEVkkrQNz6xEe/rtlpbCagp8kWNkcZCbQr8ZVFhNJLWyNshNod8sMRVWE5GqykfqPPts9badOMhNx/SbJcbROyJyjGojdSpJ0yA3HdNPQndPCPalg+Gk7fH94XbpoAJfJBLVRupAOITT09P5g9y0p98qGikjEp16Rur86lfp+9PVnn4M0vapEelw+Xx9I3U6nUK/mTL2q0kkdqUnbRctgkOHarfvtJE6lcxKugOpNzEersYdGS6ZIGUlDGxo3XF8HToSmVIjJ22hM0fqVKI9/Zkojth5bChcjfvCSLh9bCisb2ahtfJxZkuWhPvjCRVzE4lcrZO25TqhkFq9dCJ3Jh5cGwJ+Ur0dgFwYuXNuE2amUuE0kYa4w6mn1j5pO2sWzJsXdyG1eqXmRK6ZnW9mw2b2pJm5mV1RxzYXmtmPzeyImT1uZmta39MqRoapHPiE9SNNmpmqVkWovXvhhgSKuYlEpvzH8IEDtdvPmwdPPNE5hdTqlfThnTnAHuC6ehqb2RLgO8C/AmcDNwNfN7N3tqh/1bWz3o6mPRSpqfhjeGgo7N2PjIRyyLV0d1efbqKTJXoi192/C3wXwOo7MXktsM/dP1W4v9fMzgM+Cdzdkk5W0656O41UhNLJXcmoRo7fQ3ZO2laStu+55cC9ZevuLqyvyMxmm1lvcQGa9yOuHfV2NO2hyJS2b28s8LNy0raStIX+fOBw2brDQK+ZHV9lm3XAaMlysGm9GdgQ6upMehsL9XYGmvSp0rSHIpOUHsOf6vh9Vxecckpnl1eoVxbG6W8Ebiq530Ozgr9Yb2fPDeGk7Svj9FeFwG/WOP0NG+C++6qP3snqLotkVqNj8BcuTGd5hVZIW+gfAuaVrZsHjLn7C5U2cPcjwJHi/TrPHdRv1twwLPPcTa07rl6c9vCGG8Lv2ImJzhhnJjJNjY7BX7VKgV8UzTh9M3Ngtbtvq9Hmr4FL3f2sknX/CLzK3d9V5+vMfJx+ElfhltJJW8m4qQqnFWXlUpZGxuknuqdvZnOB00pWLTGzs4Hn3P0JM9sI9Lv7fyg8/hVg0My+APw9cBHw74H3tK3T1ermPzYEh+9rTxllBb5kWD0D2nK5cEjn8sv1Y7hc0od3ziWMuS8qHnv/JrAGOBlYVHzQ3feZ2XuALwIfJxyb/7C7t2+45p71FSZKIdwf3RuO7zfjKlwRqaieAW2LFsU/mXlSEh294+73u7tVWNYUHl/j7hdW2OZN7j7b3V/n7lva2ul2XYUrIlVpQNv0pW3IZrLaeRWuiFS1YUM4Vl8e/BrQNjWFfiNafRWuvixE6lIc0DY4GMbed/oUh82k0G9Us6/CVclkkWnp6QmF0vbtCxdnZa1w2nRFM2SzXWY8ZLPa6J3iVbiNjN5RyWQRaYLUlFZOpeJVuEsHYc5iOL4/3C4dbHy4pkomi0ibaU9/pmZyodRUV5gsXqxxZyIyJe3pt9NMTtrWWzJZRKRJFPpJUclkEUmAQj9JusJERNpMoZ8kXWEiIm2m0E+SrjARkTbT6J2YqGSyiEyDRu+klQJfRFpMoS8ikiEK/XbL2OE0EYmLQr8dVFRNRCKhE7mtpqJqItJiOpEbExVVE5GIKPRbbXh4cuAX5fOwXdMrikj7KPRbSUXVRCQyCv1WUlE1EYmMQr/VVFRNRCKi0G81FVUTkYgo9FtNRdVEJCIap99uKqomIk2mcfqxqPSFqsAXkQQp9JtNJRdEJGI6vNNMKrkgIgnQ4Z2kqOSCiEROod9MKrkgIpFT6DeLSi6ISAoo9JtFJRdEJAUU+s2kkgsiEjmFfjOp5IKIRE6h30xz56rkgohELYrQN7PrzGy/mb1oZrvM7C012q4xMy9bXmxnf49RfjHWG98YTtY+/DAcOAD79sGmTQp8EYnCrKQ7YGbvB24CrgV2AZ8A7jaz09396SqbjQGnl9xPZkhMtYuxhobgvvu0dy8i0YlhT/+/AF9z983u/nNC+P8/4Joa27i7HypZDrelp+V0MZaIpEyioW9mxwHnAPcW17l7vnB/eY1N55rZr83sgJndaWZvqPEas82st7gAzdv11sVYIpIySe/pvxroAsr31A8D86ts8yjhV8DlwIcI/4YdZragSvt1wGjJcnCGfQ50MZaIpFDSod8wd9/p7lvdfbe7PwC8F3gG+GiVTTYCfSVLtS+HxuhiLBFJoaRD/zfAUWBe2fp5wKF6nsDdJ4CfAKdVefyIu48VF6B5NY51MZaIpEyioe/uLwEPARcX15lZrnB/Zz3PYWZdwFnAU63oY026GEtEUibpPX0IwzX/k5n9qZktA74MzAE2A5jZVjPbWGxsZp8xs3eY2alm9mbgVuC1wNfb3nPNfysiKZP4OH13v93MTgL+knDydjfwrpJhmIuA0iEyfwR8rdD2t4RfCisKwz3br6cnXHy1aZPmvxWR6GnmLBGRlNPMWSIiUpFCX0QkQxT6IiIZotAXEckQhb6ISIYo9EVEMkShLyKSIQp9EZEMUeiLiGSIQl9EJEMU+iIiGaLQFxHJEIW+iEiGKPRFRDJEoS8ikiEKfRGRDFHoi4hkiEJfRCRDFPoiIhmi0BcRyRCFvohIhij0RUQyRKEvIpIhCn0RkQxR6IuIZIhCX0QkQxT6IiIZotAXEckQhb6ISIYo9EVEMkShLyKSIQp9EZEMUeiLiGSIQl9EJEMU+iIiGRJF6JvZdWa238xeNLNdZvaWKdpfaWa/KLT/qZld2qq+uTubN29m1apVLF++nFWrVrF582bcvVUvKSLSMrOS7oCZvR+4CbgW2AV8ArjbzE5396crtF8B3AasA74NfBDYZmZvdvdHmtk3d2fNmjXceuut5PP5V9bfdddd3H///WzZsgUza+ZLioi0lCW9x2pmu4Afuftg4X4OOAD8T3f/qwrtbwfmuPtlJet+AOx292vreL1eYHR0dJTe3t6abTdv3syHP/zhYwK/qKuri69+9atcc801U72kiEhLjY2N0dfXB9Dn7mO12iZ6eMfMjgPOAe4trnP3fOH+8iqbLS9tX3B3tfZmNtvMeosL0FNv/+64446KgQ9w9OhR7rjjjnqfSkQkCkkf03810AUcLlt/GJhfZZv5DbZfB4yWLAfr7dwzzzxT8/Fnn3223qcSEYlC0qHfDhuBvpJlQb0bnnTSSTUfP/HEE2fUMRGRdks69H8DHAXmla2fBxyqss2hRtq7+xF3HysuwHi9nVu9ejVdXV0VH+vq6mL16tX1PpWISBQSDX13fwl4CLi4uK5wIvdiYGeVzXaWti+4pEb7aVuzZg1XXXXVpODv6uriqquu4uqrr272S4qItFTiQzYJwzW/aWYPAj8kDNmcA2wGMLOtwIi7ryu03wQ8YGafAr4DfAA4F/hIsztmZmzZsoULLriAO+64g2effZYTTzyR1atXc/XVV2u4poikTuJDNgHMbBD4NOFk7G5grbvvKjx2P7Df3deUtL8SuBFYDPwS+DN3v6vO16p7yKaISBo0MmQzitBvJ4W+iHSa1IzTFxGR9lLoi4hkiEJfRCRDFPoiIhmi0BcRyRCFvohIhsRwcVYixsZqjmoSEUmNRvIsi+P0+2mg0qaISIoscPeRWg2yGPoGnEIDhdcKeghfFgumsW3S1PdkqO/JyGrfe4AnfYpQz9zhncIbUvObsJKSOjvjU13xFhv1PRnqezIy3Pe62utErohIhij0RUQyRKFfvyPAXxRu00Z9T4b6ngz1vYbMncgVEcky7emLiGSIQl9EJEMU+iIiGaLQFxHJEIV+CTO7zsz2m9mLZrbLzN4yRfsrzewXhfY/NbNL29XXCn2pu+9mtsbMvGx5sZ39LenL+WY2bGZPFvpxRR3bXGhmPzazI2b2uJmtaX1PK/ajob4X+l3+vruZzW9Tl4v9WGdmPzKzcTN72sy2mdnpdWyX+Od9On2P5fNuZh8zs4fNbKyw7DSzd0+xTdPfc4V+gZm9H7iJMFzqzcAe4G4ze02V9iuA24BvAG8CtgHbzOzMtnT42L401PeCMeDkkuW1re5nFXMI/b2unsZmtgT4DvCvwNnAzcDXzeydLepfLQ31vcTpHPveP93kfk3lAmAIeBtwCdAN3GNmc6ptENHnveG+F8TweT8IXA+cA5wL3AfcaWZvqNS4Ze+5u2sJw1Z3AbeU3M8RyjVcX6X97cC3y9b9APhKCvq+Bng+6fe8Qr8cuGKKNn8NPFK27p+Af05B3y8stDsh6fe6rF8nFfp1fo020Xzep9H3KD/vhb49B/zHdr7n2tMHzOw4wrfvvcV17p4v3F9eZbPlpe0L7q7RviWm2XeAuWb2azM7YGZV9zYiFMX7PkO7zewpM/uemf3bpDsD9BVun6vRJtb3vZ6+Q2SfdzPrMrMPEH4t7qzSrCXvuUI/eDXQBRwuW38YqHa8dX6D7VtlOn1/FLgGuBz4EOFzsMPMFrSqk01U7X3vNbPjE+hPI54CrgX+XWE5ANxvZm9OqkNmliMcIvu+uz9So2ksn/dXNND3aD7vZnaWmf2OcMXtV4DV7v7zKs1b8p5nrsqmgLvvpGTvwsx2AHuBjwJ/nlS/Op27P0oIoKIdZvY64JPAnyTTK4aAM4HzEnr9mair75F93h8lnIvqA94HfNPMLqgR/E2nPf3gN8BRYF7Z+nnAoSrbHGqwfatMp+/HcPcJ4CfAac3tWktUe9/H3P2FBPozUz8koffdzG4BLgPe7u5TTSwUy+cdaLjvx0jy8+7uL7n74+7+kLuvIwwE+HiV5i15zxX6hP8RwEPAxcV1hZ+OF1P9eNvO0vYFl9Ro3xLT7PsxzKwLOItw+CF2UbzvTXQ2bX7fLbgFWA1c5O776tgsivd9mn0vf46YPu85YHaVx1rznid99jqWBXg/8CLwp8Ay4O+A3wLzCo9vBTaWtF8BTACfAl4PfBZ4CTgzBX3/DPAO4FTCEM/bgBeAMxLo+1xC8J1NGIXxycJ/Lyo8vhHYWtJ+CfB74AuF9/0/Ay8D70xB3z9BOK58GuGwxM2EX2kXt7nfXwKeJwx/nF+yHF/SJsrP+zT7HsXnvfB5OB9YTPjS2QjkgUva+Z639Y8k9gUYBH5NOMmyC3hryWP3A1vK2l9JOEZ3BHgEuDQNfQe+WNL2EGHc+5sS6veFhcAsX7YUHt8C3F9hm58U+v9/gTVp6DvwZ8DjhcB5lnCtwdsT6HelPnvp+xjr5306fY/l804Yb7+/0I+nCSNzLmn3e67SyiIiGaJj+iIiGaLQFxHJEIW+iEiGKPRFRDJEoS8ikiEKfRGRDFHoi4hkiEJfRCRDFPoiLRbL9I4ioNAXaanIpncUURkGkZkws5OAnwJ/6+6fL6xbQaij8m5Coa/3uPuZJdv8E2HKxHe1v8eSddrTF5kBd3+GMCvTZ83sXDPrAf6BMGfxvxDvNIOSUZo5S2SG3P0uM/sa8C3gQULp53WFh2tO7+jpnPhFUkx7+iLN8V8JO1FXAle5+5GE+yNSkUJfpDleB5xC+JtaXLK+06Z3lJTT4R2RGTKz44BbgdsJE1583czOcvenCVPbXVq2SZqnd5SU0+gdkRkys78B3gcMAL8DHgBG3f2ywpDNR4Ah4O+Bi4C/JYzouTuhLkuGKfRFZsDMLgS+R5j28P8U1i0G9gDXu/uXC22+CJwBHAQ+5+5b2t9bEYW+iEim6ESuiEiGKPRFRDJEoS8ikiEKfRGRDFHoi4hkiEJfRCRDFPoiIhmi0BcRyRCFvohIhij0RUQyRKEvIpIhCn0RkQz5/1GPYD8t02iXAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "axes = plt.figure(figsize=(4, 4), dpi=100).gca()\n", + "historyPG = onp.asarray(history)\n", + "updatesPG = onp.asarray(updates) \n", + "axes.scatter(historyGD[:,0], historyGD[:,1], lw=0.5, color='blue')\n", + "axes.scatter(historyNt[:,0], historyNt[:,1], lw=0.5, color='orange')\n", + "axes.scatter(historyPG[:,0], historyPG[:,1], lw=0.5, color='red')\n", + "axes.scatter([0], [0], lw=0.25, color='black') # target at 0,0\n", + "axes.set_xlabel('x0')\n", + "axes.set_ylabel('x1')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "So, PG does best...\n", + "Also most _diagonal_, i.e., direct path to target\n", + "\n", + "Difference also shows in first update steps: compare...\n", + "GD mostly \"down\" along x1, Newton mostly along x0, PG is nicely diagonal." + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Diagonal lengths (larger is better): GD 1.053930, Nt 1.264911, PG 1.356443 \n" + ] + } + ], + "source": [ + "def mag(x):\n", + " return np.sqrt(np.sum(np.square(x)))\n", + "\n", + "def one_len(x):\n", + " return np.dot( x/mag(x), np.array([1,1])) \n", + "\n", + "print(\"Diagonal lengths (larger is better): GD %f, Nt %f, PG %f \" % \n", + " (one_len(updatesGD[0]) , one_len(updatesNt[0]) , one_len(updatesPG[0])) )\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Approximate Inversions\n", + "\n", + "We can use BFGS if inverse of z is not readily available\n", + "\n", + "Just a bit ugly for calling `fmin_l_bfgs_b()` from scipy:" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "BFGS optimization test run, find x such that y=[2,2]:\n" + ] + }, + { + "data": { + "text/plain": [ + "array([2.00000003, 1.41421353])" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "#import numpy as np\n", + "import jax.numpy as np\n", + "\n", + "def fun_z_inv_opt(target_y, x_ini):\n", + " # a bit ugly, we switch to pure scipy here inside each iteration for BFGS\n", + " import numpy as np\n", + " from scipy.optimize import fmin_l_bfgs_b\n", + " target_y = onp.array(target_y)\n", + " x_ini = onp.array(x_ini)\n", + "\n", + " def fun_z_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", + " ret = onp.sum( onp.square(y-target_y) )\n", + " return ret\n", + " \n", + " ret = fmin_l_bfgs_b(lambda x: fun_z_opt(x,target_y), x_ini, approx_grad=True )\n", + " #print( ret ) # full BFGS details\n", + " return ret[0]\n", + "\n", + "print(\"BFGS optimization test run, find x such that y=[2,2]:\")\n", + "fun_z_inv_opt([2,2], [3,3])\n" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "PG iter 0: [2.09999967 2.50998022]\n", + "PG iter 1: [1.46999859 2.10000011]\n", + "PG iter 2: [1.02899871 1.75698602]\n", + "PG iter 3: [0.72029824 1.4699998 ]\n", + "PG iter 4: [0.50420733 1.22988982]\n", + "PG iter 5: [0.35294448 1.02899957]\n", + "PG iter 6: [0.24705997 0.86092355]\n", + "PG iter 7: [0.17294205 0.72030026]\n", + "PG iter 8: [0.12106103 0.60264817]\n", + "PG iter 9: [0.08474171 0.50421247]\n" + ] + } + ], + "source": [ + "x = np.asarray([3.,3.])\n", + "eta = 0.3\n", + "history = [x]; updates = []\n", + "\n", + "for i in range(10): \n", + " # same as before, Newton step for L(y)\n", + " y = fun_z(x)\n", + " GL = jax.grad(fun_L)(y)\n", + " y += -eta * np.matmul( np.linalg.inv( jax.jacobian(jax.jacobian(fun_L))(y) ) , GL)\n", + "\n", + " # optimize for inverse physics, assuming we dont have access to an inverse for fun_z\n", + " x = fun_z_inv_opt(y,x)\n", + " history.append(x)\n", + " updates.append( history[-2] - history[-1] )\n", + " print( \"PG iter %d: \"%i + format(x) )\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Nice! It works, just like PG. Not much point plotting this, it's basiclly the PG version, but let's measure the difference..." + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Diff analytic PG versus approximated: 0.000022\n" + ] + } + ], + "source": [ + "historyPGa = onp.asarray(history)\n", + "updatesPGa = onp.asarray(updates) \n", + "\n", + "print(\"Diff analytic PG versus approximated: %f\" % (np.sum(np.abs(historyPGa-historyPG))) )\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Conclusions \n", + "That concludes simple example...\n", + "\n", + "Main takeaways:\n", + "* PGs work neatly\n", + "* much faster convergence\n", + " \n", + "Coming up, complex example." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Next Steps\n", + "\n", + "Try\n", + "- X\n", + "- Y" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.8.5" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/physgrad-discuss.md b/physgrad-discuss.md new file mode 100644 index 0000000..06cf6b7 --- /dev/null +++ b/physgrad-discuss.md @@ -0,0 +1,14 @@ +Discussion +======================= + +The training via physical gradients, +... **TODO** ... + +## Summary + +✅ Pro: +- ... + +❌ Con: +- ... + diff --git a/physgrad.md b/physgrad.md new file mode 100644 index 0000000..86a62d1 --- /dev/null +++ b/physgrad.md @@ -0,0 +1,17 @@ +Physical Gradients +======================= + +It's getting better and better... + +```{figure} resources/placeholder.png +--- +height: 220px +name: pg-training +--- +TODO, visual overview of PG training +``` + +## Derivation + +... **TODO** ... +