2023-03-28 01:04:30 +02:00
{
"cells": [
{
"cell_type": "markdown",
"id": "9785abf6",
"metadata": {},
"source": [
"We have:\n",
"$$\n",
"\\nabla f = (I + \\beta A) \\vec s - \\vec y\n",
" = \\vec s - \\vec y + \\beta A \\vec s\n",
"$$ \n",
"We can substitute in the partial derivatives:\n",
"$$ \n",
"\\left(\\begin{matrix}\n",
"s_0 - y_0 + \\beta (s_0 - s_1) \\\\\n",
"s_1 - y_1 + \\beta (2s_1 - s_0 - s_2) \\\\\n",
"\\vdots \\\\\n",
"s_N - y_N + \\beta (s_N - s_{N - 1})\n",
"\\end{matrix}\\right) \n",
"= \\vec s - \\vec y + \\beta A \\vec s\n",
"$$\n",
"\n",
"Adding $\\vec y - \\vec s$ to both sides and then dividing by $\\beta$ yields:\n",
"$$ \n",
"\\left(\\begin{matrix}\n",
"s_0 - s_1 \\\\\n",
"2s_1 - s_0 - s_2 \\\\\n",
"\\vdots \\\\\n",
"s_N - s_{N - 1}\n",
"\\end{matrix}\\right) \n",
"= A \\vec s\n",
"$$\n",
"\n",
"Solving for $A$ is trivial and yields:\n",
"$$ \n",
"A = \\left(\\begin{matrix}\n",
"1 & -1 & & & & & \\\\\n",
"-1 & 2 & -1 & & & & \\\\\n",
" & -1 & 2 & -1 & & & \\\\\n",
" & & \\ddots & \\ddots & \\ddots & & \\\\\n",
" & & & -1 & 2 & -1 & \\\\\n",
" & & & & -1 & 2 & -1\\\\\n",
" & & & & & -1 & 1\n",
"\\end{matrix}\\right) \n",
"$$"
]
},
{
"cell_type": "code",
2023-03-29 03:05:50 +02:00
"execution_count": 1,
2023-03-28 01:04:30 +02:00
"id": "e7d69a97",
"metadata": {},
2023-03-28 02:33:56 +02:00
"outputs": [],
2023-03-28 01:04:30 +02:00
"source": [
"import numpy as np\n",
2023-03-29 22:12:50 +02:00
"import scipy.sparse as sparse\n",
2023-03-29 03:05:50 +02:00
"import modules.tridiagonal as td\n",
2023-03-29 22:12:50 +02:00
"import matplotlib.pyplot as plot\n",
"import matplotlib.pyplot as plot\n",
"import modules.tridiagonal as td\n",
"import modules.image_tools as image_tools\n",
"from PIL import Image # PIL is an image processing module\n",
2023-03-28 01:04:30 +02:00
"\n",
"def sample_with_noise(points, vectorized_g):\n",
" \"\"\"\n",
" Sample a function at a given set of points,\n",
" introducing noise in the process.\n",
"\n",
" Returns a tuple of the real values of the function\n",
" and the values with introduced noise.\n",
"\n",
" For performance reasons, this function takes in\n",
" a vectorized version of `g`. If you want to pass\n",
" a non vectorized function use `np.vectorize`. Eg:\n",
" ```py\n",
" sample_with_noise(points, np.vectorize(g))\n",
" ```\n",
" \"\"\"\n",
" gs = vectorized_g(points)\n",
" ys = gs + 0.05 * np.random.randn(len(gs))\n",
"\n",
" return (gs, ys)\n",
"\n",
2023-03-28 02:00:58 +02:00
"def p5(x):\n",
2023-03-28 01:04:30 +02:00
" \"\"\"\n",
" Vectorized version of the function\n",
" described in the pdf.\n",
" \"\"\"\n",
2023-03-28 02:00:58 +02:00
" n = 5\n",
2023-03-28 01:04:30 +02:00
" result = 0\n",
"\n",
" for k in range(n + 1):\n",
" result += x**k\n",
"\n",
" return result/(n + 1)\n",
"\n",
2023-03-29 03:05:50 +02:00
"def create_A(n):\n",
" \"\"\"\n",
" Creates the (n+1) x (n + 1) A matrix.\n",
" (see the latex above for it's definition)\n",
" \"\"\"\n",
" a = 2*np.ones(n + 1)\n",
" a[0] = 1\n",
" a[-1] = 1\n",
"\n",
" c = -np.ones(n)\n",
" e = -np.ones(n)\n",
"\n",
" return td.create(a, c, e)\n",
"\n",
2023-03-28 01:04:30 +02:00
"def smooth(data, beta):\n",
" \"\"\"\n",
" Smooth noisy data (y) by means of solving \n",
" a minimization problem.\n",
"\n",
" Args:\n",
" data (numpy.array): noisy data to be smoothed (y)\n",
" beta (float): parameter >= 0 that balances fit and smoothing\n",
"\n",
" Returns:\n",
" numpy array of smoothed data (s)\n",
" \"\"\"\n",
" # We need to solve for s in:\n",
" # (I + βA)s - y = 0\n",
" # <=> (I + βA)s = y\n",
" #\n",
" # Because A and I are tridiagonal, I + βA\n",
" # is tridiagonal as well (linear combinations\n",
" # of tridiagonal matrices are themselves\n",
" # tridiagonal)\n",
" #\n",
" # This means we can use our existing implementation\n",
" # of tridiagonal linear equation system solving.\n",
" # \n",
" # We start by constructing \n",
" # iba := I + βA\n",
2023-03-29 03:05:50 +02:00
" n = len(data) - 1\n",
2023-03-28 01:04:30 +02:00
"\n",
" # I + βA\n",
2023-03-29 03:05:50 +02:00
" iba = td.add(\n",
" td.identity(n + 1),\n",
" td.scale(beta, create_A(n))\n",
2023-03-28 01:04:30 +02:00
" )\n",
"\n",
2023-03-29 03:05:50 +02:00
" data_smoothed = td.solve(*iba, data)\n",
2023-03-28 01:04:30 +02:00
"\n",
" return data_smoothed"
]
},
{
"cell_type": "code",
2023-03-29 03:05:50 +02:00
"execution_count": 2,
2023-03-28 01:04:30 +02:00
"id": "2b9f2b79",
"metadata": {},
"outputs": [
{
"data": {
2023-03-29 22:12:50 +02:00
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAHgCAYAAABZ+0ykAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAACLmUlEQVR4nOzdd3gUVdvA4d/uprcN6QklCb2E3gkQEOkiiAqoL10/UVEUlWIDfJVir6CilFcQESkqHelVutJrIJRASID0unu+P2IWlhTSN9k893XtBXv2zMwzs4F5ctpolFIKIYQQQggrobV0AEIIIYQQxUmSGyGEEEJYFUluhBBCCGFVJLkRQgghhFWR5EYIIYQQVkWSGyGEEEJYFUluhBBCCGFVJLkRQgghhFWR5EYIIYQQVkWSGyGKaN68eWg0GtPLwcEBPz8/OnfuzLRp04iKiir0vo8fP87kyZO5cOFC8QWchy1btqDRaNiyZUupHO9+OnXqRKdOnQq0TVk7B0vK6VqsXr2ayZMnWywmIUqDJDdCFJO5c+eye/duNmzYwNdff02TJk2YMWMG9erV488//yzUPo8fP86UKVNKLbkpa2bOnMnMmTMLtE2zZs3YvXs3zZo1K6GoyrfVq1czZcoUS4chRImysXQAQliLkJAQWrRoYXr/6KOP8sorr9C+fXv69+/PmTNn8PX1tWCE5U/9+vULvI2bmxtt2rQpgWgKJzk5GQcHBzQajaVDEaLCkJYbIUpQtWrV+Pjjj4mPj+fbb781le/fv59BgwYRFBSEo6MjQUFBPPHEE1y8eNFUZ968eTz++OMAdO7c2dTtNW/ePAA2bNhA3759qVKlCg4ODtSsWZNnn32W6OjofMV28uRJevTogZOTE15eXowaNYr4+Pgc6/7555906dIFNzc3nJycCA0NZePGjWZ1Jk+ejEaj4dixYzzxxBPo9Xp8fX0ZMWIEsbGxZnVTUlKYOHEiwcHB2NnZUblyZV544QVu375tVi+nbqlZs2bRuHFjXFxccHV1pW7durzxxhumz3Pqihk2bBguLi6cPXuWXr164eLiQtWqVXn11VdJTU012//ly5d57LHHcHV1xd3dnaeeeop9+/aZXfvcZHVRrl+/nhEjRuDt7Y2Tk5PpGIsXL6Zt27Y4Ozvj4uJC9+7dOXTokNk+zp8/z6BBgwgICMDe3h5fX1+6dOnC4cOHTXU0Gk2OXUtBQUEMGzYs1/iGDRvG119/bdpH1iurZXDJkiW0bt0avV6Pk5MT1atXZ8SIEXmesxBlkbTcCFHCevXqhU6nY9u2baayCxcuUKdOHQYNGoSHhweRkZHMmjWLli1bcvz4cby8vOjduzdTp07ljTfe4OuvvzZ1s9SoUQOAc+fO0bZtW55++mn0ej0XLlzgk08+oX379hw5cgRbW9tcY7p+/TphYWHY2toyc+ZMfH19WbhwIaNHj85Wd8GCBQwZMoS+ffsyf/58bG1t+fbbb+nevTvr1q2jS5cuZvUfffRRBg4cyMiRIzly5AgTJ04EYM6cOQAopejXrx8bN25k4sSJdOjQgX/++YdJkyaxe/dudu/ejb29fY5x//zzzzz//PO8+OKLfPTRR2i1Ws6ePcvx48fv+z2kp6fz8MMPM3LkSF599VW2bdvGf//7X/R6Pe+88w4AiYmJdO7cmZs3bzJjxgxq1qzJ2rVrGThw4H33f7cRI0bQu3dvfvzxRxITE7G1tWXq1Km89dZbDB8+nLfeeou0tDQ+/PBDOnTowN69e02tVL169cJgMPDBBx9QrVo1oqOj2bVrV7bErzDefvttEhMT+fXXX9m9e7ep3N/fn927dzNw4EAGDhzI5MmTcXBw4OLFi2zatKnIxxWi1CkhRJHMnTtXAWrfvn251vH19VX16tXL9fOMjAyVkJCgnJ2d1eeff24qX7JkiQLU5s2b84zBaDSq9PR0dfHiRQWo3377Lc/648ePVxqNRh0+fNisvGvXrmbHS0xMVB4eHqpPnz5m9QwGg2rcuLFq1aqVqWzSpEkKUB988IFZ3eeff145ODgoo9GolFJq7dq1OdZbvHixAtR3331nKgsLC1NhYWGm96NHj1bu7u55ntvmzZuzXbOhQ4cqQP3yyy9mdXv16qXq1Kljev/1118rQK1Zs8as3rPPPqsANXfu3DyPnfWzMGTIELPyiIgIZWNjo1588UWz8vj4eOXn56cGDBiglFIqOjpaAeqzzz7L8ziAmjRpUrbywMBANXToUNP7nK7FCy+8oHL6r/+jjz5SgLp9+3aexxaiPJBuKSFKgVLK7H1CQgLjx4+nZs2a2NjYYGNjg4uLC4mJiZw4cSJf+4yKimLUqFFUrVoVGxsbbG1tCQwMBLjvPjZv3kyDBg1o3LixWfmTTz5p9n7Xrl3cvHmToUOHkpGRYXoZjUZ69OjBvn37SExMNNvm4YcfNnvfqFEjUlJSTLPGsloC7u0+efzxx3F2ds7W3XW3Vq1acfv2bZ544gl+++23fHfBQWY3TJ8+fbLFdndX4NatW3F1daVHjx5m9Z544ol8HwcyW6/utm7dOjIyMhgyZIjZdXRwcCAsLMzUhebh4UGNGjX48MMP+eSTTzh06BBGo7FAxy6sli1bAjBgwAB++eUXrly5UirHFaIkSHIjRAlLTEwkJiaGgIAAU9mTTz7JV199xdNPP826devYu3cv+/btw9vbm+Tk5Pvu02g00q1bN5YtW8a4cePYuHEje/fuZc+ePQD33UdMTAx+fn7Zyu8tu379OgCPPfYYtra2Zq8ZM2aglOLmzZtm23h6epq9z+piyoopJiYGGxsbvL29zeppNBr8/PyIiYnJNe7BgwczZ84cLl68yKOPPoqPjw+tW7dmw4YNeZ4vgJOTEw4ODtliS0lJMb2PiYnJcdB3QQeC+/v7m73Puo4tW7bMdh0XL15sStI0Gg0bN26ke/fufPDBBzRr1gxvb29eeumlXMdDFZeOHTuyYsUKUxJWpUoVQkJCWLRoUYkeV4iSIGNuhChhq1atwmAwmAbGxsbGsnLlSiZNmsSECRNM9VJTU7MlCrk5evQof//9N/PmzWPo0KGm8rNnz+Zre09PT65du5at/N4yLy8vAL788stcZyAV9Mbv6elJRkYGN27cMEtwlFJcu3bN1IKQm+HDhzN8+HASExPZtm0bkyZN4qGHHuL06dOmlqvC8vT0ZO/evdnKc7pWebl3ZlTWdfz111/vG2NgYCA//PADAKdPn+aXX35h8uTJpKWl8c033wCZSdm9A6GBPBPD/Ojbty99+/YlNTWVPXv2MG3aNJ588kmCgoJo27ZtkfYtRGmSlhshSlBERASvvfYaer2eZ599Fsi88Smlsg2a/f777zEYDGZl97Z6ZMm6ed67j7tnZOWlc+fOHDt2jL///tus/KeffjJ7Hxoairu7O8ePH6dFixY5vuzs7PJ1zCxZA5AXLFhgVr506VISExOzDVDOjbOzMz179uTNN98kLS2NY8eOFSiOnISFhREfH8+aNWvMyn/++eci7bd79+7Y2Nhw7ty5XK9jTmrXrs1bb71Fw4YNOXjwoKk8KCiIf/75x6zupk2bSEhIuG8suf1M3VsnLCyMGTNmAGSb0SVEWSctN0IUk6NHj5rGUkRFRbF9+3bmzp2LTqdj+fLlplYKNzc3OnbsyIcffoiXlxdBQUFs3bqVH374AXd3d7N9hoSEAPDdd9/h6uqKg4MDwcHB1K1blxo1ajBhwgSUUnh4ePDHH3/kq3sG4OWXX2bOnDn07t2b9957zzRb6uTJk2b1XFxc+PLLLxk6dCg3b97ksccew8fHhxs3bvD3339z48YNZs2aVaDr1LVrV7p378748eOJi4sjNDTUNFuqadOmDB48ONdtn3nmGRwdHQkNDcXf359r164xbdo09Hr9fVt88mPo0KF8+umn/Oc//+G9996jZs2arFmzhnXr1gGg1Rbu98GgoCDeffdd3nzzTc6fP0+PHj2oVKkS169fZ+/evTg7OzNlyhT++ecfRo8ezeOPP06tWrWws7Nj06ZN/PPPP2atfIMHD+btt9/mnXfeISwsjOP
2023-03-28 01:04:30 +02:00
"text/plain": [
"<Figure size 640x480 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Reusable function to plot the data\n",
2023-03-28 02:35:48 +02:00
"def plot_denoise_results(beta, N, ax):\n",
2023-03-28 01:04:30 +02:00
" # Create the data\n",
" x = np.linspace(-1, 1, num=N+1)\n",
2023-03-28 02:00:58 +02:00
" g, y = sample_with_noise(x, p5)\n",
2023-03-29 03:05:50 +02:00
" s = smooth(y, beta)\n",
2023-03-28 01:04:30 +02:00
"\n",
2023-03-28 03:23:12 +02:00
" ax.set_title(f\"β={beta} and {N=}\")\n",
2023-03-28 01:04:30 +02:00
" ax.set_xlabel(\"x\")\n",
" ax.set_ylabel(\"y\")\n",
" ax.plot(x, y, \".\", label=\"raw data points\")\n",
" ax.plot(x, g, \"-\", label=\"original function\")\n",
2023-03-28 02:29:13 +02:00
" ax.plot(x, s, \"-\", label=\"denoised values\")\n",
2023-03-28 01:04:30 +02:00
" ax.legend()\n",
"\n",
"\n",
"N = 100\n",
2023-03-28 02:35:48 +02:00
"beta = 10\n",
2023-03-28 01:04:30 +02:00
"\n",
"fig, ax = plot.subplots()\n",
2023-03-29 22:12:50 +02:00
"fig.suptitle(\"Data denoising results\")\n",
2023-03-28 01:04:30 +02:00
"\n",
"# Plot the data\n",
2023-03-28 02:35:48 +02:00
"plot_denoise_results(beta, N, ax)"
2023-03-28 01:04:30 +02:00
]
},
{
"cell_type": "markdown",
"id": "f2e5c860",
"metadata": {},
"source": [
2023-03-29 22:12:50 +02:00
"We notice that finding a point where $\\nabla f(s) = 0$ is the same as\n",
"finding a point where $\\frac 1 \\beta \\nabla f(s) = \\nabla \\frac 1 \\beta f(s) = 0$.\n",
"We can use this to our advantage, by defining $h(s) = \\frac 1 \\beta f(s)$ and\n",
"finding the points where $\\nabla h(s) = 0$ instead. When $\\beta \\to \\infty$, we have:\n",
2023-03-28 01:04:30 +02:00
"$$\n",
"\\lim_{\\beta \\to \\inf} h _\\beta (s) =\n",
"\\lim_{\\beta \\to \\inf} \n",
" \\frac 1 {2\\beta} \\sum_{k=0}^n (y_k - s_k)^2\n",
" +\\frac 1 2 \\sum_{k=1}^n (s_k - s_{k - 1})^2\n",
" = \\frac 1 2 \\sum_{k=1}^n (s_k - s_{k - 1})^2\n",
"$$\n",
"\n",
2023-03-29 22:12:50 +02:00
"It now becomes trivial to verify that $\\nabla f(s) = As$.\n",
"Looking at the value of $A$ we notice that any vector with all values\n",
"equal to some constant are solutions for this equation.\n",
"This means the result will be a constant function\n",
"(which might differ based on the linear equation solver we are using)."
2023-03-28 01:04:30 +02:00
]
},
{
"cell_type": "code",
2023-03-29 03:05:50 +02:00
"execution_count": 3,
2023-03-28 01:04:30 +02:00
"id": "d7f7dbdf",
"metadata": {
2023-03-28 03:23:12 +02:00
"scrolled": false
2023-03-28 01:04:30 +02:00
},
"outputs": [
{
"data": {
2023-03-29 22:12:50 +02:00
"image/png": "iVBORw0KGgoAAAANSUhEUgAABQsAAAPLCAYAAAD4+99NAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOzdd1yV5f/H8dfhsIeggIoTcpt775GlqZk50rJfZmrfzJZZrqwcmTu1oQ0rsbQyczQ0R5qamoYzV6aGYooiDkBAxjn37w/iJAIKCByE9/PxOI/gvq/7uj/n7ng+3Nd9DZNhGAYiIiIiIiIiIiJS5DnYOwAREREREREREREpGNRYKCIiIiIiIiIiIoAaC0VERERERERERORfaiwUERERERERERERQI2FIiIiIiIiIiIi8i81FoqIiIiIiIiIiAigxkIRERERERERERH5lxoLRUREREREREREBFBjoYiIiIiIiIiIiPxLjYUiIiIiIiIiIiICqLFQRERERERERERE/qXGQhEREUknODgYk8lke7m6ulK6dGnat2/PlClTiIiIyHHd27dvZ/z48Vy5ciX3Ar6J8ePHYzKZ8uVct5J6XU+ePJmnxxRmGV2P/P5M5Za4uDgmTJhA9erVcXV1xdfXl0cffZTLly/bOzQREREpwtRYKCIiIplasGABv/32G+vXr2fu3LnUq1ePadOmUaNGDX7++ecc1bl9+3YmTJhwxzXs5IauXbvy22+/ERAQkKfHFDV34mfKMAweffRR5syZw9NPP81PP/3EG2+8wZIlS3j77bftHZ6IiIgUYY72DkBEREQKrlq1atGoUSPb77169eKll16iVatW9OzZk2PHjlGqVCk7Rnhn8ff3x9/fP8+PyWtxcXG4u7vbO4w72ubNm/n+++9ZunQpvXv3BrD13I2NjbVzdCIiIlKUqWehiIiIZEuFChV4++23iYmJ4aOPPgLg+PHjPPnkk1SpUgV3d3fKli1Lt27dOHDgQJpjx48fz4gRIwAICgqyDXPetGlTluu4mVWrVlGvXj1cXFwICgpi5syZGZY7duwY/fr1o2TJkri4uFCjRg3mzp2bLlaTycShQ4d49NFH8fb2plSpUgwcOJCoqKh0dW7dupUOHTrg5eWFu7s7LVq0YNWqVWnKZDSE9sKFC/zvf/+jfPnyuLi44O/vT8uWLW09NzM6Jruxfffdd9SpUwcXFxfuuusu3nnnnSwPz04tt2fPHnr37k3x4sWpVKlSlq9jVt7jgAEDCAwMzPTct4ovs8/Urc6bmax+FnNaP8DSpUspXrw4PXr0sG3bsmUL58+f55577rnl8almz57NypUrb1omOTk53TaLxZLlc4iIiEjRop6FIiIikm1dunTBbDazZcsWAM6ePYuvry9Tp07F39+fS5cusXDhQpo2bcrevXupVq0aAIMHD+bSpUu89957LF++3Da0tmbNmuzbty9LdWRmw4YNdO/enebNm/P1119jsViYPn0658+fT1Pu8OHDtGjRwtboWbp0adauXcsLL7xAZGQk48aNS1O+V69e9O3bl0GDBnHgwAHGjBkDwGeffWYrs3nzZu677z7q1KnDp59+iouLC/PmzaNbt2589dVX9O3bN9O4H3/8cfbs2cNbb71F1apVuXLlCnv27OHixYu3/P+QldjWrFlDz549adOmDUuWLCE5OZmZM2emuy630rNnTx555BGGDBlCbGxstq7j7bzHW7nZZ6pPnz45Om9WP8+38762b99O06ZNbedbu3YtI0eOpEOHDjzwwANZfv+7du1i1KhRfPPNNzz00ENp9u3Zs4d+/fpx7Ngx7r33XhYvXkxMTAy9e/dm3759NGjQgC+//JIqVapk+XwiIiJSBBgiIiIiN1iwYIEBGCEhIZmWKVWqlFGjRo0M9yUnJxuJiYlGlSpVjJdeeinNvhkzZhiAERoaetMYblZHRpo2bWqUKVPGiI+Pt22Ljo42SpQoYVz/J0+nTp2McuXKGVFRUWmOf+655wxXV1fj0qVLhmEYxrhx4wzAmD59eppyQ4cONVxdXQ2r1Wrb1qxZM6NkyZJGTExMmvhr1apllCtXzlY29bpe/949PT2NYcOGZfq+MjomO7E1btzYKF++vJGQkGDbFhMTY/j6+hpZ+VMw9VxvvPFGmu1ZvY5ZeY9PPPGEUbFixUzPfb2Mrkdmn6lbnTerMvss5rT++Ph4w9HR0Rg3bpwxceJEAzAAo0KFCsbp06ezHVu/fv0MJycnY8WKFWn21alTx1izZo1x8eJFo2/fvkanTp2Mxo0bG6NGjTIuXrxoLFmyxGjZsmW24xcREZHCTcOQRUREJEcMw7D9nJyczOTJk6lZsybOzs44Ojri7OzMsWPHOHLkSJbqu506YmNjCQkJoWfPnri6utq2e3l50a1bN9vv165dY8OGDfTo0QN3d3eSk5Ntry5dunDt2jV27NiRpu4HH3wwze916tTh2rVrthWhY2Nj2blzJ71798bT09NWzmw28/jjj/PPP/9w9OjRTGNv0qQJwcHBTJo0iR07dpCUlHTri5WN2Hbt2sVDDz2Es7OzrZynp2ea65IVvXr1sv2c3et4O+/xduT0vFn9LOa0/j179pCcnEyTJk147LHHWLt2LRMmTCAmJoY2bdpw9epVACIjI9OsSp7Ry9HRkS+//JKkpCT69Olj6zH6zz//UK5cOTp16kSJEiVYuHAhx48ft/WYLFGiBH369AGwnU9EREQENAxZREREciA2NpaLFy9Su3ZtAIYPH87cuXMZNWoUbdu2pXjx4jg4ODB48GDi4+OzVOft1HH58mWsViulS5dOt+/6bRcvXiQ5OZn33nuP9957L8O6IiMj0/zu6+ub5ncXFxcAW0yXL1/GMIwMVysuU6aM7byZWbJkCZMmTeKTTz7h9ddfx9PTkx49ejB9+vQM309OYstoEZrsLkxz/fvL7nW8nfd4O3J63qx+FnNa/++//w6kNDb6+flx11130bFjR6pWrcqjjz7Kjh07uPfee/Hy8mL+/Pm3fJ9r1qxh2bJldO/e3faZMAwDB4f/+gW4uLgQEBDAxYsXsVqtafZd3/AvIiIiosZCERERybZVq1ZhsVho164dAIsWLaJ///5Mnjw5TbnIyEh8fHyyVOft1FG8eHFMJhPnzp1Lt+/6bcWLF7f1+Hv22WczrCsoKChL8V5fp4ODA+Hh4en2nT17FgA/P79Mj/fz82POnDnMmTOHsLAwvv/+e0aPHk1ERARr1qzJViwZxWYymTKcnzCja3Uz1y80kt3reKv36OrqSkJCQro6bmy4za6cXtusfhZzWv/vv//OXXfdlennIrWR2cXFhcGDB9/0Pa5atYoff/yR3r1789VXX+HomPLnfbly5Th16hTr1q2jefPmTJ48GScnJ+Lj43n++eeZPn06P//8MxaLBS8vr5ueQ0RERIoWNRaKiIhItoSFhfHKK6/g7e3N008/DaQ0JKX2aku1atUqzpw5Q+XKldNsv7H3W6rs1HEjDw8PmjRpwvLly5kxY4ZtKHJMTAw//PCDrZy7uzvt27dn79691KlTJ83Q3Jzy8PCgadOmLF++nJkzZ+Lm5gaA1Wpl0aJFlCtXjqpVq2aprgoVKvDcc8+xYcMGtm3bliuxNWrUiJUrVzJz5kzb+7169So//vhjjuu9neuY0XsMDAwkIiKC8+fP23o8JiYmsnbt2izVmdln6lbnzUxOPovZqf/3339P1xPVMAw++eQTatWqRc2aNW96/PVmzJhhW0gntaEw9T0EBwfz2GOP8eeff9KsWTOWL19OYmIivXr1wtPTk2rVqvHtt99m+VwiIiJSNKixUERERDJ18OBB21x0ERER/PrrryxYsACz2cyKFSvw9/cH4IEHHiA4OJjq1atTp04ddu/ezYwZMyhXrly6OlOHLr/zzjs88cQTODk5Ua1atWzVkZE333yT+++/n/vuu4+XX34Zi8XCtGnT8PDw4NKlS7Zy77zzDq1ataJ169Y888wzBAYGEhMTw/Hjx/nhhx/YuHFjtq/TlClTuO+++2jfvj2vvPIKzs7OzJs3j4MHD/L
2023-03-28 01:04:30 +02:00
"text/plain": [
"<Figure size 1280x960 with 9 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fig, axs = plot.subplots(\n",
" 3,\n",
" 3,\n",
" layout=\"constrained\",\n",
" figsize=(12.8, 9.6)\n",
")\n",
"\n",
"fig.suptitle(r\"Data denoising results as $\\beta \\to \\infty$\")\n",
"\n",
"N = 100\n",
"\n",
"for i in range(len(axs.flat)):\n",
" ax = axs.flat[i]\n",
2023-03-28 02:33:56 +02:00
" plot_denoise_results(4**i, N, ax)"
2023-03-28 01:04:30 +02:00
]
},
2023-03-28 03:23:12 +02:00
{
"cell_type": "markdown",
"id": "3c9f50cc",
"metadata": {},
"source": [
"When $\\beta \\to 0$ we have:\n",
"$$\n",
"\\lim_{\\beta \\to 0} \\nabla f(s) =\n",
"\\lim_{\\beta \\to 0} (I + \\betaA)s - y = s - y\n",
"$$\n",
"\n",
2023-03-29 22:12:50 +02:00
"This means when solving for $\\nabla f(s) = 0$, we are solving\n",
"the equation $s - y = 0$, which has the solution $s = y$.\n",
"The resulting function will then be one equal to the data containing the noise."
2023-03-28 03:23:12 +02:00
]
},
2023-03-28 01:04:30 +02:00
{
"cell_type": "code",
2023-03-29 03:05:50 +02:00
"execution_count": 4,
2023-03-28 01:04:30 +02:00
"id": "dfc485d5",
2023-03-29 22:12:50 +02:00
"metadata": {
"lines_to_next_cell": 1
},
2023-03-28 03:23:12 +02:00
"outputs": [
{
"data": {
2023-03-29 22:12:50 +02:00
"image/png": "iVBORw0KGgoAAAANSUhEUgAABQsAAAPLCAYAAAD4+99NAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOzdd3gUZdfH8e+m9wAJJKGk0IN0AlKkiYCIqICC5RFR9FVREZAiViwPRVGxtwcEC8UCNpAiUgREutJUSkIogRBKEtKzO+8fS1aWFBIg2U3y+1xXLsjsPTNnxnUPc/YuJsMwDERERERERERERKTSc3F0ACIiIiIiIiIiIuIcVCwUERERERERERERQMVCEREREREREREROUfFQhEREREREREREQFULBQREREREREREZFzVCwUERERERERERERQMVCEREREREREREROUfFQhEREREREREREQFULBQREREREREREZFzVCwUERERERERERERQMVCEREREREREREROUfFQhERESnQrFmzMJlMth8vLy9CQ0Pp3r07kydPJjEx8ZKPvX79eiZOnMiZM2euXMBFmDhxIiaTqUzOdTF59zUuLq5U96nICrofZf2eulLS09N54YUXaNy4MV5eXgQFBXHHHXdw+vTpK3L8s2fPMnLkSGrWrImXlxctW7Zk3rx5V+TYIiIiUjGpWCgiIiJF+uSTT/jtt99Yvnw57777Li1btmTq1KlER0fz888/X9Ix169fzwsvvFDuCjtXQt++ffntt98ICwsr1X0qm/L4njIMgzvuuIPp06fz4IMP8tNPP/Hcc88xf/58XnvttStyjgEDBjB79myef/55fvrpJ9q2bcsdd9zBnDlzrsjxRUREpOJxc3QAIiIi4tyaNm1KTEyM7feBAwcyatQorrnmGgYMGMDevXsJCQlxYITlS/Xq1alevXqp71Pa0tPT8fHxcXQY5drq1av5/vvv+eqrr7j11lsBbD1309LSLvv4ixcvZvny5cyZM4c77rjDdvyDBw8yduxYBg8ejKur62WfR0RERCoW9SwUERGREgsPD+e1114jNTWVDz/8EIB9+/Zx77330qBBA3x8fKhVqxb9+vVjx44ddvtOnDiRsWPHAhAVFWUb5rxq1apiH6MoixYtomXLlnh6ehIVFcW0adMKbLd3717uvPNOatSogaenJ9HR0bz77rv5YjWZTOzatYs77riDwMBAQkJCuO+++0hOTs53zLVr19KjRw/8/f3x8fGhY8eOLFq0yK5NQUNoT5w4wf/93/9Rp04dPD09qV69Op06dbL13Cxon5LG9t1339G8eXM8PT2pW7cub775ZrGHZ+e127p1K7feeitVq1alXr16xb6PxbnGoUOHEhkZWei5LxZfYe+pi523MMV9L17q8QG++uorqlatSv/+/W3b1qxZw/Hjx7n22msvun+eN954g2+//Tbf9oULF+Ln58dtt91mt/3ee+/l6NGj/P7778U+h4iIiFQe6lkoIiIil+SGG27A1dWVNWvWAHD06FGCgoKYMmUK1atX59SpU8yePZurr76abdu20ahRIwDuv/9+Tp06xdtvv82CBQtsQ2ubNGnC9u3bi3WMwqxYsYKbb76ZDh06MG/ePMxmM6+88grHjx+3a7d79246duxoK3qGhoaydOlSRowYQVJSEs8//7xd+4EDBzJ48GCGDRvGjh07mDBhAgAzZ860tVm9ejU9e/akefPmzJgxA09PT9577z369evH3LlzGTx4cKFx33333WzdupX//ve/NGzYkDNnzrB161ZOnjx50f8OxYltyZIlDBgwgC5dujB//nxyc3OZNm1avvtyMQMGDOD222/noYceIi0trUT38XKu8WKKek8NGjToks5b3Pfz5VzX+vXrufrqq23nW7p0KePGjaNHjx7ceOONxb7+zZs3M378eL788ktuueUW2/adO3cSHR2Nm5v9P/mbN29ue71jx47FPo+IiIhUEoaIiIhIAT755BMDMDZt2lRom5CQECM6OrrA13Jzc43s7GyjQYMGxqhRo+xee/XVVw3AiI2NLTKGoo5RkKuvvtqoWbOmkZGRYduWkpJiVKtWzTj/nz29e/c2ateubSQnJ9vt/+ijjxpeXl7GqVOnDMMwjOeff94AjFdeecWu3fDhww0vLy/DYrHYtrVv396oUaOGkZqaahd/06ZNjdq1a9va5t3X86/dz8/PGDlyZKHXVdA+JYmtbdu2Rp06dYysrCzbttTUVCMoKMgozj8H88713HPP2W0v7n0szjXec889RkRERKHnPl9B96Ow99TFzltchb0XL/X4GRkZhpubm/H8888bL774ogEYgBEeHm4cOnSoxLHdeeedhru7u7Fw4ULb9gYNGhi9e/fO1/7o0aMGYEyaNKnEcYuIiEjFp2HIIiIicskMw7D9PTc3l0mTJtGkSRM8PDxwc3PDw8ODvXv3smfPnmId73KOkZaWxqZNmxgwYABeXl627f7+/vTr18/2e2ZmJitWrKB///74+PiQm5tr+7nhhhvIzMxkw4YNdse+6aab7H5v3rw5mZmZthWh09LS+P3337n11lvx8/OztXN1deXuu+/m8OHD/P3334XG3q5dO2bNmsXLL7/Mhg0byMnJufjNKkFsmzdv5pZbbsHDw8PWzs/Pz+6+FMfAgQNtfy/pfbyca7wcl3re4r4XL/X4W7duJTc3l3bt2nHXXXexdOlSXnjhBVJTU+nSpQtnz54FICkpyW5V8oJ+3NzcmDNnDjk5OQwaNMiux2hRQ7idZYVwERERcS4ahiwiIiKXJC0tjZMnT9KsWTMARo8ezbvvvsv48ePp2rUrVatWxcXFhfvvv5+MjIxiHfNyjnH69GksFguhoaH5Xjt/28mTJ8nNzeXtt9/m7bffLvBYSUlJdr8HBQXZ/e7p6Qlgi+n06dMYhlHgasU1a9a0nbcw8+fP5+WXX+Z///sfzz77LH5+fvTv359XXnmlwOu5lNgKWoSmpAvTnH99Jb2Pl3ONl+NSz1vc9+KlHn/jxo2AtdgYHBxM3bp16dWrFw0bNuSOO+5gw4YNXHfddfj7+/Pxxx9f9DqXLFnCN998w80332x7TwQFBRX4vjt16hQA1apVu+hxRUREpPJRsVBEREQuyaJFizCbzXTr1g2Azz//nCFDhjBp0iS7dklJSVSpUqVYx7ycY1StWhWTycSxY8fyvXb+tqpVq9p6/D3yyCMFHisqKqpY8Z5/TBcXFxISEvK9dvToUQCCg4ML3T84OJjp06czffp04uPj+f7773nyySdJTExkyZIlJYqloNhMJlOB8xMWdK+Kcn5PtJLex4tdo5eXF1lZWfmOcWHhtqQu9d4W9714qcffuHEjdevWLfR9kVdk9vT05P777y/yGhctWsSPP/7Irbfeyty5c21zFDZr1oy5c+eSm5trN29h3iItTZs2LfK4IiIiUjlpGLKIiIiUWHx8PGPGjCEwMJAHH3wQsBaS8nq15Vm0aBFHjhzJt/+Fvd/ylOQYF/L19aVdu3YsWLCAzMxM2/bU1FR++OEH2+8+Pj50796dbdu20bx5c2JiYvL9XNhbrzjnvvrqq1mwYIHdNVksFj7//HNq165Nw4YNi3Ws8PBwHn30UXr27MnWrVtLFEdhscXExPDtt9+SnZ1t23727Fl+/PHHSz7u5dzHgq4xMjKSxMREu6JmdnY2S5cuLVY8hb2nLnbewlzKe7Ekx9+4cWO+nqiGYfC///2Ppk2b0qRJkyL3P9+rr75qW0jn/KJg//79OXv2LN98841d+9mzZ1OzZk3b4ioiIiIi51PPQhERESnSzp07bXPRJSYm8uuvv/LJJ5/g6urKwoULqV69OgA33ngjs2bNonHjxjRv3pwtW7bw6quvUrt27XzHzBu6/Oabb3LPPffg7u5Oo0aNSnSMgrz00ktcf/319OzZkyeeeAKz2czUqVPx9fW1Db3MO+8111xD586defjhh4mMjCQ1NZV9+/bxww8/8Msvv5T4Pk2ePJmePXvSvXt3xowZg4eHB++99x47d+5k7ty5hc4Pl5ycTPfu3bnzzjt
2023-03-28 03:23:12 +02:00
"text/plain": [
"<Figure size 1280x960 with 6 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fig, axs = plot.subplots(\n",
" 2,\n",
" 3,\n",
" layout=\"constrained\",\n",
" figsize=(12.8, 9.6)\n",
")\n",
"\n",
"fig.suptitle(r\"Data denoising results as $\\beta \\to 0$\")\n",
"\n",
"N = 100\n",
"\n",
"for i in range(len(axs.flat)):\n",
" ax = axs.flat[i]\n",
" plot_denoise_results(2**-i, N, ax)"
]
},
2023-03-29 22:12:50 +02:00
{
"cell_type": "markdown",
"id": "7cd79ba5",
"metadata": {
"lines_to_next_cell": 0
},
"source": [
"We know from the last lecture that picking any\n",
"$\\alpha$ which satisfies $\\alpha < \\frac 2 {\\lambda_H}$\n",
"leads to a converging sequence, where $\\lambda_H$ is\n",
"the largest eigenvalue of the hessian.\n",
"Let $\\lambda_A$ be the largest eigenvalue of $A$.\n",
"Our hessian is $I + \\beta A$, therefore we have $Ax = \\lambda_A x$ \n",
"(for some vector $x$, from the definition of the eigenvalue for $\\lambda_A$).\n",
"We can multiply by $\\beta$ and add $x$ to both sides to get\n",
"$(I + \\beta A)I = (1 + \\beta\\lambda_A) x$,\n",
"which means $1 + \\beta\\lambda_A$ is indeed an eigenvalue of the hessian.\n",
"\n",
"Is it the biggest one though? \n",
"Let's assume a bigger eigenvalue $\\lambda_M$ exists (for the hessian).\n",
"This means that $(I + \\beta A)y = \\lambda_M y$ (for some vector $y$).\n",
"By distributivity we then also have $y + \\beta A y = \\lambda_M y$,\n",
"which makes it easy to then notice that $\\frac {\\lambda_M-1} \\beta$\n",
"is an eigenvalue for $A$.\n",
"Earlier we have assumed that $\\lambda_M > \\lambda_H$,\n",
"therefore $\\frac {\\lambda_M - 1} \\beta > \\frac {\\lambda_H - 1} \\beta = \\lambda_A$,\n",
"which would imply $\\lambda_A$ is not the biggest eigenvalue for $A$,\n",
"hence a contradiction.\n",
"This let's us conclude that $\\lambda_H$ is the biggest eigenvalue for the hessian,\n",
"and therefore our sequence will converge given\n",
"$\\alpha < \\frac 2 {\\lambda _H} = \\frac 2 {1 + \\beta \\lambda_A}$."
]
},
2023-03-29 03:05:50 +02:00
{
"cell_type": "code",
2023-03-29 22:12:50 +02:00
"execution_count": 5,
2023-03-29 03:05:50 +02:00
"id": "cef74c06",
"metadata": {
2023-03-29 22:12:50 +02:00
"lines_to_next_cell": 1,
2023-03-29 03:05:50 +02:00
"scrolled": true
},
"outputs": [
{
"data": {
2023-03-29 22:12:50 +02:00
"image/png": "iVBORw0KGgoAAAANSUhEUgAABQsAAAPLCAYAAAD4+99NAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOzdd1yV1R/A8c9lb1BABFHAgZl7b3HkTHOVM7fmNrWcufdIc5RmmdJwVY7SzJErdzhzmwSSiiKk7H3P7w/i/rwyVeCCft+v1329uOc5z3m+z3PvhS/nPuccjVJKIYQQQgghhBBCCCGEeOUZGToAIYQQQgghhBBCCCFE/iCdhUIIIYQQQgghhBBCCEA6C4UQQgghhBBCCCGEEP+RzkIhhBBCCCGEEEIIIQQgnYVCCCGEEEIIIYQQQoj/SGehEEIIIYQQQgghhBACkM5CIYQQQgghhBBCCCHEf6SzUAghhBBCCCGEEEIIAUhnoRBCCCGEEEIIIYQQ4j/SWSiEEELH19cXjUbDmTNndGW7d+9mxowZhgsqG3F4enrSt2/fPI3HkDQajd61SH3dAgMD8zyWVatW4evr+8LtPH1OhjRjxgw0Gs0z75efzuFV0rhxYxo3bmzoMDK0ZcsWypcvj6WlJRqNhgsXLuTq8W7cuMFbb71FkSJFsLW1pUqVKqxbty7H2j927BgDBw6kevXqmJubZ/m7Z+XKlbz22muYm5vj5eXFzJkzSUxMzLF4ckLfvn3x9PTMsl7jxo3RaDS0atUqzbbAwEA0Gg0ff/zxC8cTGRnJ+PHjadGiBc7Ozln+bjl37hxvvPEGNjY2ODg40KlTJ/7+++906xaE10MIIYThSWehEEKITO3evZuZM2caOoxM49i+fTtTp07N44jyjzfffJOTJ0/i6uqa58fOqc7C/GTgwIGcPHnymfc7efIkAwcOzIWIRGZWrVrFqlWrDB1Guh4+fEivXr0oVaoUe/bs4eTJk3h7e+fa8eLi4mjTpg1+fn4sWrSIn3/+mYYNGzJgwAB++eWXHDnGgQMH+O233yhRogT16tXLtO7cuXN5//336dSpE3v37mXYsGHMmzeP4cOH50gshrJ3714OHjyYa+2HhYXxxRdfEB8fT4cOHTKte/36dRo3bkxCQgLff/8969at4+bNmzRs2JCHDx/q1X1ZXw8hhBA5z8TQAQghhHg1xcTEYGVllSNtVa1aNUfaKaicnZ1xdnbOsl5OXvOXmbu7O+7u7s+8X506dXIhGpGR1Pfz66+/buhQMnTz5k0SExN599138fHxyZE2M/sc//HHH/z999989dVXurutmzRpwoYNG9i/fz9vvvnmCx9/6tSpTJ8+HYCPP/6Yw4cPp1svLCyMOXPmMGjQIObNmwek3JmXmJjIlClTGD16dL5+7TLi7e1NUlIS48ePx8/P77nuQs6Kh4cHjx49QqPREBoaytq1azOsO23aNMzNzdm1axd2dnYAVK9enTJlyvDxxx+zcOFC4OV9PYQQQuQOubNQCCFEhvr27ctnn30GpAyxTH2kDjlTSrFq1SqqVKmCpaUlhQoV4u23304z/Klx48ZUqFCB33//nXr16mFlZUX//v2BlCF6LVq0wNXVFUtLS8qVK8fEiROJjo7OdhxPDkN++PAhZmZm6d5peP36dTQaDStWrNCV3b9/n8GDB+Pu7o6ZmZluWFZSUlKm16ZDhw54eHig1WrTbKtduzbVqlXTPf/hhx+oXbs29vb2WFlZUbJkSd35ZyYiIoJBgwbh6OiIjY0NrVq14ubNm2nqpTcMObNrHhERwYcffoiXlxdmZmYUK1aM0aNH611zAK1Wy8qVK3Wvr4ODA3Xq1OHnn38GUq77lStXOHLkiO41yWooX3bPCeCvv/6iR48eFClSBHNzc8qVK6d7H6Q6fPgwGo2GTZs28dFHH+Hm5oadnR1vvPEGN27cSNPmunXrqFy5MhYWFhQuXJiOHTty7do1vTrpDUM+ePAgjRs3xtHREUtLS0qUKEHnzp2JiYnR1cloePihQ4cYOnQoTk5OODo60qlTJ+7du6fXfnx8PB988AFFixbFysqKRo0acfbs2WwPsY+Pj2fWrFmUK1cOCwsLHB0dadKkCSdOnNDViYuLY9KkSXqv+/Dhw3n8+LFeW56enrRt25Zdu3ZRtWpV3edy165duvMqV64c1tbW1KpVS2/aAkj5vNrY2HDlyhWaNWuGtbU1zs7OjBgxQu96AXz22Wc0atSIIkWKYG1tTcWKFVm0aFGaYZGZvZ/TG4a8evVqKleujI2NDba2trz22mtMnjxZr87ly5dp3749hQoVwsLCgipVqvD111/r1XnW99fT16FBgwYAdO3aFY1Goxfnzz//TN26dbGyssLW1pbmzZunuaM19b147tw53n77bQoVKkSpUqUyPGZQUBCAXqdPZGQkkZGRmJqaZhpvdhkZZe/fhz179hAXF0e/fv30yvv164dSih07dmS6/8OHDxk2bBivv/46NjY2FClShKZNm3L06FG9ek8O/126dCleXl7Y2NhQt25dTp06laZdX19fypYtq/ud8s0332TrfFKZmpoyd+5czp49y5YtW55p3+xK/X2alaSkJHbt2kXnzp11HYWQ0tnYpEkTtm/frit70ddDCCHEq0XuLBRCCJGhqVOnEh0dzY8//qj3T2zqcNfBgwfj6+vLqFGjWLhwIf/++y+zZs2iXr16XLx4ERcXF90+wcHBvPvuu4wfP5558+bp/uH866+/aNOmDaNHj8ba2prr16+zcOFC/vjjD90wr6zieJKzszNt27bl66+/ZubMmXr/2K5fvx4zMzN69uwJpHQU1qpVCyMjI6ZNm0apUqU4efIkc+bMITAwkPXr12d4bfr370/79u05ePAgb7zxhq78+vXr/PHHH7oOyZMnT9K1a1e6du3KjBkzsLCw4Pbt21kOYVNK0aFDB06cOMG0adOoWbMmx48fp3Xr1pnu96T0rnlMTAw+Pj7cuXOHyZMnU6lSJa5cucK0adO4dOkSv/32m+6f1L59+/Ldd98xYMAAZs2ahZmZGefOndN1Sm7fvp23334be3t73TBQc3PzHDmnq1evUq9ePUqUKMGSJUsoWrQoe/fuZdSoUYSGhurubEo1efJk6tevz9q1a4mIiGDChAm0a9eOa9euYWxsDMD8+fOZPHky3bt3Z/78+YSFhTFjxgzq1q2Ln58fZcqUSTfuwMBA3nzzTRo2bMi6detwcHDg7t277Nmzh4SEhCzv1hw4cCBvvvkmGzdu5J9//mHcuHG8++67eu+Bfv36sWXLFsaPH0/Tpk25evUqHTt2JCIiItO2IaXDoHXr1hw9epTRo0fTtGlTkpKSOHXqFEFBQdSrV0937Q8cOMCkSZNo2LAhf/75J9OnT+fkyZOcPHlS77W7ePEikyZN4qOPPsLe3p6ZM2fSqVMnJk2axIEDB5g3bx4ajYYJEybQtm1bAgICsLS01O2fmJhImzZtGDx4MBMnTuTEiRPMmTOH27dvs3PnTl09f39/evTooevAvHjxInPnzuX69etp5tnL6HfI0zZv3sywYcMYOXIkH3/8MUZGRty6dYurV6/q6ty4cYN69epRpEgRVqxYgaOjI9999x19+/blwYMHjB8/Xq/N7Ly/njZ16lRq1arF8OHDmTdvHk2aNNF16GzcuJGePXvSokULNm3aRHx8PIsWLaJx48YcOHBA18mYqlOnTnTr1o0hQ4ak6dR/+r0AYGJiQkxMDDdu3GDy5MlotVq6dOmiV1er1ab7ZcfTNBpNhueYmcuXLwNQsWJFvXJXV1ecnJx02zPy77//AjB9+nSKFi1KVFQU27dv112jpzuIP/vsM1577TWWLVsGpFz/Nm3aEBAQgL29PZDSUdivXz/at2/PkiVLCA8PZ8aMGcTHx2e7ExRSOn8//vhjpkyZQufOnTPtiM3qi6dUxsbGz3yXor+/P7GxsVSqVCnNtkqVKrF//37i4uKwsLB44ddDCCHEK0YJIYQQ/1m/fr0ClJ+fn65s+PDhKr0/FydPnlSAWrJkiV75P//8oywtLdX48eN1ZT4
2023-03-29 03:05:50 +02:00
"text/plain": [
"<Figure size 1280x960 with 6 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"def smooth_iterative(data, beta, num_iterations):\n",
" \"\"\"\n",
" Smooth noisy data (y) by means of solving a minimization problem\n",
" iteratively with the gradient descent method.\n",
"\n",
" Args:\n",
" data (numpy.array): of the noisy data to be smoothed (y)\n",
" beta (float): parameter >= 0 that balances fit and smoothing\n",
" num_iterations (int): number of gradient descent iterations (>0)\n",
"\n",
" Returns:\n",
" numpy array of smoothed data (s)\n",
" \"\"\"\n",
"\n",
" n = len(data) - 1\n",
"\n",
" A = create_A(n)\n",
"\n",
" # Our initial eigenvalue guess\n",
" initial_guess = np.zeros(n + 1)\n",
" initial_guess[0] = 1\n",
"\n",
" max_eigenvalue = td.largest_eigenvalue(\n",
" *A, \n",
" initial_guess,\n",
" 30\n",
" )\n",
"\n",
" alpha = 1.5 / (1 + beta * max_eigenvalue)\n",
" x = np.zeros(n + 1)\n",
"\n",
" for _ in range(num_iterations):\n",
" gradient = td.multiply_vector(\n",
" *td.add(\n",
" td.identity(n + 1), \n",
" td.scale(beta, A)\n",
" ),\n",
" x\n",
" ) - data\n",
" # print(f\"{alpha=}\")\n",
" # print(f\"{gradient=}\")\n",
" # print(f\"{x=}\")\n",
" x = x - alpha * gradient\n",
"\n",
" return x\n",
"\n",
"fig, axs = plot.subplots(\n",
" 3,\n",
" 2,\n",
" # 1,\n",
" # 1,\n",
" layout=\"constrained\",\n",
" figsize=(12.8, 9.6)\n",
")\n",
"\n",
"N = 100\n",
"beta = 10\n",
"\n",
"fig.suptitle(f\"Iterative vs direct denoising comparison for β={beta} and {N=}\")\n",
"\n",
"for i in range(len(axs.flat)):\n",
"# for _ in range(1):\n",
" ax = axs.flat[i]\n",
" iterations = 3**(i + 1)\n",
" # ax = axs\n",
" # iterations = 3 ** 5\n",
"\n",
" # Create the data\n",
" x = np.linspace(-1, 1, num=N+1)\n",
" g, y = sample_with_noise(x, p5)\n",
" sd = smooth(y, beta) # Indirect\n",
" si = smooth_iterative(y, beta, iterations) # Direct\n",
"\n",
" # Plot the data\n",
" ax.set_title(f\"{iterations=}\")\n",
" ax.set_xlabel(\"x\")\n",
" ax.set_ylabel(\"y\")\n",
" ax.plot(x, y, \".\", label=\"raw data points\")\n",
" ax.plot(x, g, \"-\", label=\"original function\")\n",
" ax.plot(x, sd, \"-\", label=\"directly denoised values\")\n",
" ax.plot(x, si, \"-\", label=\"iteratively denoised values\")\n",
" ax.legend()"
]
},
2023-03-28 03:23:12 +02:00
{
"cell_type": "code",
2023-03-29 22:12:50 +02:00
"execution_count": 6,
"id": "8079a686",
"metadata": {
"lines_to_next_cell": 0
},
"outputs": [],
"source": [
"def hessian(beta, nx, ny):\n",
" \"\"\"\n",
" Compute the Hessian of the image denoising problem, given by\n",
" `H = I + beta C^T C`.\n",
"\n",
" Args:\n",
" beta (float): parameter >= 0 that balances fit and smoothing\n",
" nx (int): number of pixels in x-direction\n",
" ny (int): number of pixels in y-direction\n",
"\n",
" Returns:\n",
" H: sparse hessian matrix (scipy.sparse.csr_array)\n",
" \"\"\"\n",
" size = nx * ny\n",
"\n",
" b_upper_diagonal = np.ones(nx)\n",
" b_upper_diagonal[-1] = 0\n",
"\n",
" b_main_diagonal = -np.ones(nx)\n",
" b_main_diagonal[-1] = 0\n",
"\n",
" # We want to reuse the tridiagonal -> csr conversion code,\n",
" # so we first define C as a tridiagonal matrix.\n",
" tridiagonal_C = td.create(\n",
" np.resize(b_main_diagonal, size),\n",
" np.resize(b_upper_diagonal, size - 1),\n",
" np.zeros(size - 1)\n",
" )\n",
"\n",
" C = td.to_csr(*tridiagonal_C)\n",
" H = sparse.identity(size, format=\"csr\") + beta * (C.T @ C)\n",
"\n",
" return H"
]
},
{
"cell_type": "code",
"execution_count": 13,
"id": "1775b494",
"metadata": {
"lines_to_next_cell": 0
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAV4AAAGxCAYAAAAqOd5NAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOydd3xb1fn/39rDsrxnbCfO3ntvCISGQAiQEKCsAoU2BZqGTb+MsEdYZY8CKWW1FCi0rAwSKIEsMp3txHtPSba29Psjv3O4ku3ETpzJfb9eftmWrq7uvdL9nOc85xmacDgcRkVFRUXlmKE93gegoqKi8ktDFV4VFRWVY4wqvCoqKirHGFV4VVRUVI4xqvCqqKioHGNU4VVRUVE5xqjCq6KionKMUYVXRUVF5RijCq+KiorKMUYVXpU2+fHHH5k7dy4ZGRkYjUbS09OZM2cOP/zwQ4f2c99996HRaA7rGFauXIlGo2HlypWH9fr2MnXqVKZOndpp26moHAxVeFVa5bnnnmPChAmUlJTw+OOPs2zZMhYvXkxpaSkTJ07k+eefb/e+rr322g6LtWD48OH88MMPDB8+/LBe39m8+OKLvPjii8f7MFROcjRqrQaVaL7//nsmT57M2Wefzccff4xer5fPBQIBzj//fD7//HO+/fZbJkyY0OZ+mpubsVqtx+KQjxhhxR5ty1pFBVSLV6UVHnnkETQaDS+99FKE6ALo9XpefPFFNBoNjz76qHxcuBN++ukn5syZQ0JCAj169Ih4TonX6+Xmm28mPT0dq9XK5MmT2bBhA926deOqq66S27Xmarjqqquw2Wzs3buXs88+G5vNRnZ2NjfffDNerzfifRYtWsSYMWNITEzEbrczfPhw/vrXv3K49ka0q6GgoACNRsMTTzzBY489Rrdu3bBYLEydOpXdu3fj9/u54447yMzMJC4ujvPPP5+qqqqIfX7wwQdMnz6djIwMLBYL/fr144477qCpqanF+7/22mv07t0bk8lE//79effdd7nqqqvo1q1bxHY+n48HH3yQvn37YjKZSElJ4Te/+Q3V1dWHdd4qnYv+0Juo/JIIBoN88803jBw5kqysrFa3yc7OZsSIEaxYsYJgMIhOp5PPXXDBBVx88cX87ne/a1U4BL/5zW/44IMPuO222zj99NPZvn07559/Pg6Ho13H6ff7mTVrFtdccw0333wz3377LQ888ABxcXHcc889cruCggKuv/56cnJygAN+6xtvvJHS0tKI7Y6UF154gcGDB/PCCy/Q0NDAzTffzLnnnsuYMWMwGAy88cYbFBYWcsstt3Dttdfy6aefytfu2bOHs88+mwULFhATE8POnTt57LHHWLt2LStWrJDbvfrqq1x//fVceOGFPP300zQ2NrJo0aIWg00oFOK8887ju+++47bbbmP8+PEUFhZy7733MnXqVNavX4/FYum0c1c5DMIqKgoqKirCQPjiiy8+6Hbz5s0LA+HKyspwOBwO33vvvWEgfM8997TYVjwnyMvLCwPh22+/PWK79957LwyEr7zySvnYN998EwbC33zzjXzsyiuvDAPhf/zjHxGvP/vss8N9+vRp85iDwWDY7/eH77///nBSUlI4FArJ56ZMmRKeMmXKQc+5te32798fBsJDhgwJB4NB+fgzzzwTBsKzZs2KeP2CBQvCQLixsbHV/YdCobDf7w+vWrUqDIQ3b94sjz09PT08ZsyYiO0LCwvDBoMh3LVrV/mYuI7/+te/IrZdt25dGAi/+OKLhzxPlaOL6mpQOSzC/3+qHu1CuPDCCw/52lWrVgFw0UUXRTw+Z86cFq6NttBoNJx77rkRjw0ePJjCwsKIx1asWMEZZ5xBXFwcOp0Og8HAPffcQ21tbYsp/5Fw9tlno9X+fDv169cPgJkzZ0ZsJx4vKiqSj+3bt49LL72U9PR0eYxTpkwBYMeOHQDs2rWLioqKFtcsJyenhZ/9P//5D/Hx8Zx77rkEAgH5M3ToUNLT01U/9gmA6mpQiSA5ORmr1cr+/fsPul1BQQFWq5XExMSIxzMyMg75HrW1tQCkpaVFPK7X60lKSmrXcVqtVsxmc8RjJpMJj8cj/1+7di3Tp09n6tSpvPbaa2RlZWE0Gvnkk0946KGHcLvd7Xqv9hB9HYxG40EfF8fpcrmYNGkSZrOZBx98kN69e2O1WikuLuaCCy6Qx9jWNROPKT+vyspKGhoa5HtFU1NTczinqNKJqMKrEoFOp+O0007jyy+/pKSkpFU/b0lJCRs2bGDGjBkR/l1oaQG3hhDXyspKunTpIh8PBAJSYDqD999/H4PBwH/+858Ikf7kk0867T2OlBUrVlBWVsbKlSullQvQ0NAQsZ3ymkVTUVER8X9ycjJJSUl8+eWXrb5nbGzsER61ypGiuhpUWnDnnXcSDoeZP38+wWAw4rlgMMjvf/97wuEwd95552Htf/LkycCB1XwlH374IYFA4PAOuhU0Gg16vT5icHC73bz99tud9h5HihioTCZTxOOvvPJKxP99+vQhPT2df/zjHxGPFxUVsXr16ojHzjnnHGprawkGg4wcObLFT58+fY7Cmah0BNXiVWnBhAkTeOaZZ1iwYAETJ07khhtuICcnh6KiIl544QXWrFnDM888w/jx4w9r/wMGDOCSSy7hySefRKfTcfrpp5OXl8eTTz5JXFxchK/0SJg5cyZPPfUUl156Kddddx21tbUsXry4hcgdT8aPH09CQgK/+93vuPfeezEYDLzzzjts3rw5YjutVsuiRYu4/vrrmTNnDldffTUNDQ0sWrSIjIyMiGt28cUX884773D22Wfzxz/+kdGjR2MwGCgpKeGbb77hvPPO4/zzzz/Wp6qiQBVelVa58cYbGTVqFE8++SQ333wztbW1JCYmMnHiRP73v/8xbty4I9r/m2++SUZGBn/96195+umnGTp0KP/4xz/41a9+RXx8fKecw+mnn84bb7zBY489xrnnnkuXLl347W9/S2pqKtdcc02nvMeRkpSUxH//+19uvvlmLrvsMmJiYjjvvPP44IMPWmTrXXfddWg0Gh5//HHOP/98unXrxh133MG///3viMU6nU7Hp59+yrPPPsvbb7/NI488gl6vJysriylTpjBo0KBjfZoqUaiZayonDKtXr2bChAm88847XHrppcf7cE4KGhoa6N27N7Nnz+bVV1893oej0k5U4VU5LixdupQffviBESNGYLFY2Lx5M48++ihxcXFs2bKlRcSCyoFFtIceeojTTjuNpKQkCgsLefrpp9m5cyfr169nwIABx/sQVdqJ6mpQOS7Y7Xa+/vprnnnmGZxOJ8nJycyYMYNHHnlEFd02MJlMFBQUMH/+fOrq6rBarYwdO5aXX35ZFd2TDNXiVVFRUTnGHNdwshdffJHc3FzMZjMjRozgu+++O56Ho6KionJMOG7C+8EHH7BgwQL+/Oc/s3HjRiZNmsSMGTMiVmdVVFRUTkWOm6thzJgxDB8+nJdeekk+1q9fP2bPns0jjzxyPA5JRUVF5ZhwXBbXfD4fGzZs4I477oh4fPr06S2ycOBA7VZl6btQKERdXR1JSUmH3VJGRUVFpbMJh8M4nU4yMzMPmgh0XIS3pqaGYDDYouBHWlpai7xzOFCYe9GiRcfq8FRUVFSOiOLi4jbrWcNxDieLtlbD4XCrFuydd97JwoUL5f+NjY3k5ORQXFyM3W4/6sepoqKi0h4cDgfZ2dmHLER0XIQ3OTkZnU7XwrqtqqpqteydyWRqNb/ebrerwquionLCcSgX6HGJajAajYwYMYKlS5dGPL506dLDLryioqKicrJw3FwNCxcu5PLLL2fkyJGMGzeOV199laKiIn73u98dr0NSUVFROSYcN+GdN28etbW13H///ZSXlzNw4EA+//xzunbterwOSUVFReWYcFKmDDscDuLi4mhsbFR9vCoqKicM7dUmtUiOyimB0n4Qf4soGbHQIR7vrELrKiqHi/oNVDnpUQotHFhRFqIbCoUIhUIttlFROZ6owqty0qMU1dbazodCIflb/K2icjxRXQ0qpxyhUEg2zXS73RQ
"text/plain": [
"<Figure size 640x480 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAV4AAAGxCAYAAAAqOd5NAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOz9Z1yc55n3jX/pM8MMdShD7713AUIICQn1Ltty7yV2ipPY60120/zEiZ3YcRx3x7Js2bIkq1i9oY5A9N57h4EBBgYGGJj/Cz9zPuv17ib/vXMnt/fm9/nwgmGYel3HdZ7H8StmRqPRyBKWsIQlLOHvBvN/9AtYwhKWsIT/27BUeJewhCUs4e+MpcK7hCUsYQl/ZywV3iUsYQlL+DtjqfAuYQlLWMLfGUuFdwlLWMIS/s5YKrxLWMISlvB3xlLhXcISlrCEvzOWCu8SlrCEJfydsVR4l/DfxkcffYSZmRkSiYSurq5v/D07O5uoqKj/1mNnZ2eTnZ39v/gK/3p0dnZiZmbGRx999Hd7ziX83wvLf/QLWMK3H7Ozs/z0pz/lk08++Zs95ltvvfU3e6y/BiqVisLCQgIDA/+uz7uE/zuxtOJdwv8y8vLy+Oyzz6iqqvqbPWZERAQRERF/s8f7S7CxsSEtLQ0XF5e/23Mu4f9eLBXeJfwv47nnnsPZ2Znnn3/+L95Xr9fzwgsv4O/vj7W1NZ6ennznO99hfHz8a/f7j1oNb7/9NrGxscjlchQKBWFhYfzzP/8z8FWrwNLSkpdeeukbz3n9+nXMzMw4fPjwf/q6/qNWw89//nPMzMyorq5m165d2Nvb4+TkxLPPPovBYKCpqYm8vDwUCgV+fn68/PLL33ivP/zhD4mLixP/u2zZMr788stvPP/4+DgPP/wwTk5OyOVyNmzYQHt7O2ZmZvz85z//2n1bWlrYs2cPrq6u2NjYEB4ezptvvvmfvrcl/J+HpcK7hP9lKBQKfvrTn3L+/HkuX778n97PaDSydetWfve733Hvvfdy+vRpnn32Wfbt20dOTg6zs7P/6f9+/vnnPPXUU6xYsYJjx45x/PhxfvCDH6DT6QDw8/Nj8+bNvPPOOywsLHztf//0pz/h4eHBtm3b/lvvb/fu3cTGxnLkyBEeffRRXnvtNX7wgx+wdetWNmzYwLFjx8jJyeH555/n6NGj4v9mZ2fRaDT86Ec/4vjx4xw4cIDMzEy2b9/Oxx9/LO63uLjIpk2b+Oyzz3j++ec5duwYqamp5OXlfeO11NfXk5ycTG1tLb///e85deoUGzZs4Lvf/S6/+MUv/lvvbwn/ABiXsIT/Jvbu3WsEjCUlJcbZ2VljQECAMSkpybi4uGg0Go3GFStWGCMjI8X9z507ZwSML7/88tce5+DBg0bA+N5774nbVqxYYVyxYoX4/emnnzY6ODj8l6/nypUrRsB47NgxcVtfX5/R0tLS+Itf/OK//N+Ojg4jYNy7d6+47Wc/+5kRMP7+97//2n3j4uKMgPHo0aPitvn5eaOLi4tx+/bt/+lzGAwG4/z8vPHhhx82xsfHi9tPnz5tBIxvv/321+7/0ksvGQHjz372M3Hb2rVrjV5eXsaJiYmv3ffpp582SiQSo0aj+S/f5xL+z8DSincJfxNYW1vz4osvUlpayqFDh/7D+5hWww888MDXbt+1axe2trbk5+f/p4+fkpLC+Pg4d911F19++SUjIyPfuE92djaxsbFf23a/8847mJmZ8dhjj/033tVX2Lhx49d+Dw8Px8zMjHXr1onbLC0tCQoK+ga74/Dhw2RkZCCXy7G0tMTKyoo///nPNDQ0iPtcu3YN+Gpl/W9x1113fe13vV5Pfn4+27ZtQyaTYTAYxM/69evR6/UUFRX9t9/nEv5+WCq8S/ib4c477yQhIYGf/OQnzM/Pf+Pvo6OjWFpafmOAZWZmhru7O6Ojo//pY9977718+OGHdHV1sWPHDlxdXUlNTeXixYtfu993v/td8vPzaWpqYn5+nvfff5+dO3fi7u7+335fTk5OX/vd2toamUyGRCL5xu16vV78fvToUXbv3o2npyf79++nsLCQkpISHnrooa/dz/S5/PvncXNz+9rvo6OjGAwG3njjDaysrL72s379eoD/8IK0hP/zsEQnW8LfDGZmZvz2t78lNzeX99577xt/d3Z2xmAwoFarv1Z8jUYjg4ODJCcn/5eP/+CDD/Lggw+i0+m4fv06P/vZz9i4cSPNzc34+voCsGfPHp5//nnefPNN0tLSGBwc5Dvf+c7f9o3+ldi/fz/+/v4cPHgQMzMzcfu/72WbPheNRvO14js4OPi1+zk6OmJhYcG99977n74nf3//v+E7WML/LiyteJfwN8Xq1avJzc3ll7/8JVNTU1/726pVq4CvCtK/xZEjR9DpdOLvfwm2trasW7eOn/zkJ8zNzVFXVyf+JpFIeOyxx9i3bx+vvvoqcXFxZGRk/C++q/8ezMzMsLa2/lrRHRwc/AarYcWKFQAcPHjwa7d//vnnX/tdJpOxcuVKKioqiImJISkp6Rs/zs7O/5vezRL+llha8S7hb47f/va3JCYmMjw8TGRkpLg9NzeXtWvX8vzzz6PVasnIyKC6upqf/exnxMfHc++99/6nj/noo48ilUrJyMhApVIxODjISy+9hL29/TdWyk899RQvv/wyZWVlfPDBB//b3udfwsaNGzl69ChPPfUUO3fupKenh1/96leoVCpaWlrE/fLy8sjIyOCHP/whWq2WxMRECgsLBfPB3Pz/Wx+9/vrrZGZmsnz5cp588kn8/PyYnJyktbWVkydP/peskiX8H4R/9HRvCd9e/FtWw7/Hnj17jMDXWA1Go9E4MzNjfP75542+vr5GKysro0qlMj755JPGsbGxr93v37Ma9u3bZ1y5cqXRzc3NaG1tbfTw8DDu3r3bWF1d/R++tuzsbKOTk5Nxenr6r3ov/xWrQa1Wf+2+999/v9HW1vYbj/HvWRxGo9H4m9/8xujn52e0sbExhoeHG99//33xuP8WGo3G+OCDDxodHByMMpnMmJubaywqKjICxtdff/0br/Whhx4yenp6Gq2srIwuLi7G9PR044svvvhXvdcl/ONhZjQupQwv4X8WhoeH8fX15ZlnnvmGqOHbhM8++4y7776bgoIC0tPT/9EvZwl/Qyy1GpbwPwa9vb20t7fzyiuvYG5uzve+971/9Ev6q3HgwAH6+vqIjo7G3NycoqIiXnnlFbKyspaK7v9ALBXeJfyPwQcffMAvf/lL/Pz8+PTTT/H09PxHv6S/GgqFgs8//5wXX3wRnU6HSqXigQce4MUXX/xHv7Ql/G/AUqthCUtYwhL+zviH0sneeust/P39kUgkJCYmcuPGjX/ky1nCEpawhL8L/mGF9+DBg3z/+9/nJz/5CRUVFSxfvpx169bR3d39j3pJS1jCEpbwd8E/rNWQmppKQkICb7/9trgtPDycrVu3/ofWfktYwhKW8D8F/5Dh2tzcHGVlZfzTP/3T125fs2YNt27d+sb9Z2dnvyazXFxcRKPR4Ozs/DVV0BKWsIQl/CNhNBqZnJzEw8Pja8KXf49/SOEdGRlhYWHhGyYgbm5u39CnA7z00ktLXqNLWMISvjXo6enBy8vrP/37P5RO9u9Xq0aj8T9cwb7wwgs8++yz4veJiQl8fHy4cuUK5ubmtLe3Mzo6SkJCAiUlJWzatAm5XM7g4CANDQ3U1dWhUqkYGhrCwcGBiIgIenp6cHZ2JjExEa1Wi1qtpqOjA3Nzc3Q6HVZWVgQEBGBtbc3MzAwDAwNs3rwZKysrzp07R1tbG+bm5sTHx6PX64mNjWVycpKxsTHUajXT09NYWFgwPT2Nubk5LS0tZGZmMjs7S3BwMG1tbbS2thIcHCzcrszMzPD39+f27dscOnRISGi7u7tRqVTo9XpCQkIoLCxEKpUSEhKCn58fFy9exNvbGysrKxwdHbl8+TKzs7MsLi7i6emJRqPBzMyM6elp6uvrCQ0NZefOnYyPj2NlZYWrqytSqZSTJ08SEBBAd3c3GzZs4Pjx44SFhREZGcmFCxfo7e3
"text/plain": [
"<Figure size 640x480 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Here so we can experiment with different picutres\n",
"image_name = \"sonic\"\n",
"image_tools.create_noisy_image(\"./images/\", image_name + \".jpg\")"
]
},
{
"cell_type": "code",
2023-03-29 22:25:18 +02:00
"execution_count": 28,
2023-03-29 22:12:50 +02:00
"id": "a5dbb0d2",
"metadata": {
"scrolled": true
},
"outputs": [
{
"data": {
2023-03-29 22:25:18 +02:00
"image/png": "iVBORw0KGgoAAAANSUhEUgAABQsAAADTCAYAAAAmlkUmAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOx9eXhkRbn+21u603unsyeTGWaGGXbZV2UYQNlB8CKoCCgiCHjV64KKP8AFRcAF9Xq9XhAQFRG5IKKgsl5kEWQH2ZktM9mT7k66k04v9fsjvDVf15xOMukZmIF6nydPktOnz6lTX1Wd+t56v69cSikFCwsLCwsLCwsLCwsLCwsLCwsLi3c83G91ASwsLCwsLCwsLCwsLCwsLCwsLCy2DFiy0MLCwsLCwsLCwsLCwsLCwsLCwgKAJQstLCwsLCwsLCwsLCwsLCwsLCws3oAlCy0sLCwsLCwsLCwsLCwsLCwsLCwAWLLQwsLCwsLCwsLCwsLCwsLCwsLC4g1YstDCwsLCwsLCwsLCwsLCwsLCwsICgCULLSwsLCwsLCwsLCwsLCwsLCwsLN6AJQstLCwsLCwsLCwsLCwsLCwsLCwsAFiy0MLCwsLCwsLCwsLCwsLCwsLCwuINWLLQwsLCwsLCYk649tpr4XK58M9//vOtLso7HqeffjoWLFjwVhcDALBy5Uq4XC5ce+21G/W9LekZ3mo41cW3v/1t3HrrrW9JeSwsLCwsLCzeWbBkoYWFhYWFhYWFxSZDW1sbHn74YRx11FEb9b3/9//+H2655ZbNVKqtH5YstLCwsLCwsHiz4H2rC2BhYWFhYWFhsTWhVCqhWCzC7/e/1UXZIuH3+7Hvvvtu9PcWLVq0GUozd+RyOQSDwbe6GBYWFhYWFhYWbzqsstDCwsLCwsJik+H0009HOBzGiy++iMMOOwyhUAhtbW249NJLAQCPPPII3v3udyMUCmHJkiW47rrrKr4/MDCAc845BzvssAPC4TCam5tx8MEH44EHHtjgXt3d3fi3f/s3RCIRxONxfOQjH8Fjjz3mGAL7z3/+E8ceeywaGhoQCASw22674Xe/+92Mz8OQ2ssuuwzf+ta3sM0228Dv9+Pee+8FANx2223Yb7/9EAwGEYlE8N73vhcPP/yw/v7zzz8Pl8uFm266SR97/PHH4XK5sOOOO1bc69hjj8Uee+wxY5muvfZaLF26FH6/H9tvvz1++ctfOp43OTmJb33rW9huu+3g9/vR1NSEj33sYxgYGKg4b8GCBTj66KNx5513Yvfdd0d9fT222247/OIXv9jgms899xyOO+44JBIJBAIB7LrrrhvY0CkMeWBgAJ/85Ccxb948XZYDDjgAd911lz7HKfTW5XLhvPPOw/XXX4/tt98ewWAQ73rXu3D77bdvULY//OEP2GWXXeD3+7Fw4UJceeWVuPjii+FyuWaqUhx00EHYaaed8H//93/Yf//9EQwG8fGPfxwAkMlk8IUvfAHbbLMN6urq0NHRgc9+9rPIZrMV17jpppuwzz77IBaLIRgMYuHChfoawPqw/ZUrV1Z877777oPL5cJ9991XtXwulwvZbBbXXXcdXC4XXC4XDjroIABTpCbLFwgE0NDQgD333BM33HDDjM9tYWFhYWFhYeEEqyy0sLCwsLCw2KQoFAo44YQTcPbZZ+OLX/wifvOb3+ArX/kKMpkMbr75Zpx//vno7OzEj3/8Y5x++unYaaedNEk2PDwMALjooovQ2tqKsbEx3HLLLTjooINw9913a4Ikm81i+fLlGB4exne/+10sXrwYd955J0466aQNynPvvffi8MMPxz777IOf/exniMVi+O1vf4uTTjoJuVwOp59++ozP9KMf/QhLlizBFVdcgWg0im233Ra/+c1v8JGPfATve9/7cMMNNyCfz+Oyyy7TZX33u9+NHXfcEW1tbbjrrrtw4oknAgDuuusu1NfX41//+hfWrVuH9vZ2FItF3H///Tj77LOnLce1116Lj33sYzjuuOPwve99D+l0GhdffDHy+Tzc7vVrwOVyGccddxweeOABfOlLX8L++++PVatW4aKLLsJBBx2Ef/7zn6ivr9fnP/300/j85z+PL3/5y2hpacFVV12FM844A4sXL8aBBx4IAHjppZew//77o7m5GT/60Y+QTCbxq1/9Cqeffjr6+vrwpS99qWq5P/rRj+KJJ57AJZdcgiVLliCVSuGJJ57A0NDQjHX/pz/9CY899hi+8Y1vIBwO47LLLsPxxx+Pl156CQsXLgQA3HnnnTjhhBNw4IEH4sYbb0SxWMQVV1yBvr6+Ga9P9PT04JRTTsGXvvQlfPvb34bb7UYul8OyZcvQ3d2Nr371q9hll13w/PPP48ILL8Szzz6Lu+66Cy6XCw8//DBOOukknHTSSbj44osRCASwatUq3HPPPbO+/3R4+OGHcfDBB2P58uX4f//v/wEAotEoAOA//uM/cP311+Nb3/oWdtttN2SzWTz33HOzqlsLCwsLCwsLC0coCwsLCwsLC4s54JprrlEA1GOPPaaPnXbaaQqAuvnmm/WxQqGgmpqaFAD1xBNP6ONDQ0PK4/Go//iP/6h6j2KxqAqFgjrkkEPU8ccfr4//53/+pwKg7rjjjorzzzrrLAVAXXPNNfrYdtttp3bbbTdVKBQqzj366KNVW1ubKpVKVe+/YsUKBUAtWrRITU5O6uOlUkm1t7ernXfeueL7o6Ojqrm5We2///762CmnnKIWLlyo/z/00EPVmWeeqRKJhLruuuuUUko9+OCDCoD661//WrUsvOfuu++uyuWyPr5y5Url8/nU/Pnz9bEbbrhhAzsopdRjjz2mAKif/vSn+tj8+fNVIBBQq1at0sfGx8dVQ0ODOuuss/Sxk08+Wfn9frV69eqKax5xxBEqGAyqVCpVUWfSBuFwWH32s5+t+mxKTbUd+QxKKQVAtbS0qEwmo4/19vYqt9utvvOd7+hje+21l5o3b57K5/P62OjoqEomk2o2091ly5YpAOruu++uOP6d73xHud3uijaulFK///3vFQD15z//WSml1BVXXKEA6DpwAvvLihUrKo7fe++9CoC699579TGnugiFQuq0007b4Lo77bSTev/73z/jM1pYWFhYWFhYzBY2DNnCwsLCwsJik8LlcuHII4/U/3u9XixevBhtbW3Ybbfd9PGGhgY0Nzdj1apVFd//2c9+ht133x2BQABerxc+nw933303XnjhBX3O/fffj0gkgsMPP7ziux/60Icq/n/11Vfx4osv4iMf+QgAoFgs6p8jjzwSPT09eOmll2Z8pmOPPRY+n0///9JLL2HdunX46Ec/WqHoC4fD+MAHPoBHHnkEuVwOAHDIIYfg9ddfx4oVKzAxMYG///3vOPzww7F8+XL87W9/AzClNvT7/Xj3u99dtQy854c//OGK0Nr58+dj//33rzj39ttvRzwexzHHHFPxzLvuuitaW1s3CHnddddd0dXVpf8PBAJYsmRJhW3uueceHHLIIZg3b17Fd08//XTkcrmK8GsTe++9N6699lp861vfwiOPPIJCoVD1XBPLly9HJBLR/7e0tFS0m2w2i3/+8594//vfj7q6On1eOBzGMcccM+v7JBIJHHzwwRXHbr/9duy0007YddddK+rxsMMOqwgd3muvvQAAH/zgB/G73/0Oa9eunfV9a8Xee++NO+64A1/+8pdx3333YXx8/E27t4WFhYWFhcXbE5YstLCwsLCwsNikCAaDCAQCFcfq6urQ0NCwwbl1dXWYmJjQ/3//+9/Hpz71Keyzzz64+eab8cgjj+Cxxx7D4YcfXkGCDA0NoaWlZYPrmccYhvqFL3wBPp+v4uecc84BAAwODs74TG1tbRX/M8TTPA4A7e3tKJfLGBkZAQAceuihAKYIwb///e8oFAo4+OCDceihh+Luu+/Wnx1wwAEVocEmeM/W1tYNPjOP9fX1IZVKoa6uboPn7u3t3eCZk8nkBtf0+/0b1Hm155Xlc8KNN96I0047DVdddRX2228/NDQ04NRTT0Vvb2/V78y2bCMjI1BKzao9TAenZ+vr68M
2023-03-29 22:12:50 +02:00
"text/plain": [
2023-03-29 22:25:18 +02:00
"<Figure size 1280x240 with 12 Axes>"
2023-03-29 22:12:50 +02:00
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"image = np.asarray(Image.open(f\"./images/{image_name}_gray_noisy.jpg\"))\n",
"\n",
"fig, axs = plot.subplots(\n",
2023-03-29 22:25:18 +02:00
" 1,\n",
" 12,\n",
2023-03-29 22:12:50 +02:00
" layout=\"constrained\",\n",
2023-03-29 22:25:18 +02:00
" figsize=(12.8, 2.4)\n",
2023-03-29 22:12:50 +02:00
")\n",
"\n",
"fig.suptitle(\"Image row denoising results\")\n",
"\n",
"ny, nx = image.shape\n",
"flat_image = image.flatten()\n",
"for i in range(len(axs.flat)):\n",
" beta = 3**i\n",
" ax = axs.flat[i]\n",
"\n",
" H = hessian(beta, nx, ny)\n",
"\n",
" smooth = sparse.linalg.spsolve(H, flat_image).reshape(image.shape)\n",
" smooth_image = Image.fromarray(smooth.astype(np.uint8), mode=\"L\")\n",
"\n",
" ax.set_title(f\"β={beta}\")\n",
2023-03-29 22:25:18 +02:00
"\n",
" # Turn off tick labels\n",
" ax.set_yticklabels([])\n",
" ax.set_xticklabels([])\n",
"\n",
2023-03-29 22:12:50 +02:00
" ax.imshow(smooth_image)"
]
},
{
"cell_type": "markdown",
2023-03-29 22:25:18 +02:00
"id": "5edb8858",
2023-03-29 22:12:50 +02:00
"metadata": {},
"source": [
"Finding the equation for $f$ when we want to denoise by column would simply require us to consider an image $p_c = p^T$:\n",
"$$\n",
"f(s) = \\frac 1 2 \\sum_{i = 1}^{N_x} \\sum_{j = 1}^{N_y} (y_{j, i} - s_{j, i})^2\n",
"+ \\beta \\frac 1 2 \\sum_{i = 1}^{N_x} \\sum_{j = 2}^{N_y} (s_{j, i} - s_{j - 1, i})^2\n",
"$$\n",
"\n",
"In this case, our hessian looks almost the same as the one in the original equation, with the difference being $B$ having the ones under the main diagonal instead of on top.\n",
"\n",
"For the actual python implementation, we can simply transpose the image before/after the denoising:"
]
},
{
"cell_type": "code",
2023-03-29 22:25:18 +02:00
"execution_count": 30,
2023-03-29 03:05:50 +02:00
"id": "d48e32df",
2023-03-28 03:23:12 +02:00
"metadata": {},
2023-03-29 22:12:50 +02:00
"outputs": [
{
"data": {
2023-03-29 22:25:18 +02:00
"image/png": "iVBORw0KGgoAAAANSUhEUgAABQsAAADTCAYAAAAmlkUmAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOx9d5hkVZn+Wzl1VVdXV8eZnswMGYYclJkBXSQLLoKKgKKCgCuuARV/gIEgYEBdl1WQMCoisCCgoAsMiGQlZ5iZntA5VFVXdVV1pfP7o3nPfHWnqrunwwQ47/P00923bt17zvlO+t7zBZtSSsHAwMDAwMDAwMDAwMDAwMDAwMDgfQ/7ti6AgYGBgYGBgYGBgYGBgYGBgYGBwfYBQxYaGBgYGBgYGBgYGBgYGBgYGBgYADBkoYGBgYGBgYGBgYGBgYGBgYGBgcG7MGShgYGBgYGBgYGBgYGBgYGBgYGBAQBDFhoYGBgYGBgYGBgYGBgYGBgYGBi8C0MWGhgYGBgYGBgYGBgYGBgYGBgYGAAwZKGBgYGBgYGBgYGBgYGBgYGBgYHBuzBkoYGBgYGBgYGBgYGBgYGBgYGBgQEAQxYaGBgYGBgYGBgYGBgYGBgYGBgYvAtDFhoYGBgYGLyPcdNNN8Fms+Gf//znti7KDoPly5dj+fLl27oY04ozzzwT8+bN29bFAAC0t7fDZrPhpptu2qLvbU912Nao1BaXX3457r777m1SHgMDAwMDA4MdC4YsNDAwMDAwMDAw2G7Q0tKCJ598Esccc8wWfe///b//h7vuumuGSrXjw5CFBgYGBgYGBhOFc1sXwMDAwMDAwMDAwIDweDw46KCDtvh7CxcunIHSTB7pdBp+v39bF8PAwMDAwMDAYIthLAsNDAwMDAwMynDmmWeipqYGb7zxBo488kgEAgG0tLTgyiuvBAA89dRT+MAHPoBAIIDFixfj5ptvLvt+X18fzj33XOy6666oqalBY2MjDj/8cDz22GObvWvjxo3493//dwSDQYTDYXzqU5/Cs88+W9EN9Z///CeOP/54RCIReL1eLF26FH/84x8nVKeRkRF873vfwy677AKv14v6+nqsWLECTzzxhL4nm83iW9/6FubPnw+3241Zs2bhvPPOQzweH/PZjzzyCGw2Gx555JGy65XcaafatnQbX7VqFb74xS8iGo2ivr4eJ510Ejo7OyfUFjfddBOWLFkCj8eDXXbZBbfcckvF+3K5HH7wgx9g5513hsfjQUNDAz7zmc+gr6+v7L558+bh2GOPxQMPPIB99tkHPp8PO++8M37zm99s9sxXXnkFJ5xwAurq6uD1erH33ntvVsdK7dbX14cvfOELaGtr02U59NBD8eCDD5a1rdX11maz4fzzz8fKlSuxyy67wO/3Y6+99sJ99923Wdn+9Kc/Yc8994TH48GCBQtw7bXX4tJLL4XNZhuvSbF8+XLsvvvu+Pvf/45DDjkEfr8fn/3sZwEAQ0ND+NrXvlbWry644AIMDw+XPeP222/HgQceiNraWvj9fixYsEA/A9gk+/b29rLvVet/1nYYHh7GzTffDJvNBpvNpl3p0+m0Lp/X60UkEsF+++2HW2+9ddx6GxgYGBgYGLw3YSwLDQwMDAwMDDZDPp/HSSedhHPOOQdf//rX8fvf/x7f+ta3MDQ0hDvvvBMXXnghZs+ejZ///Oc488wzsfvuu2PfffcFAAwODgIALrnkEjQ3NyOVSuGuu+7C8uXL8dBDD2mSYnh4GCtWrMDg4CB++MMfYtGiRXjggQdwyimnbFaeVatW4SMf+QgOPPBAXHfddaitrcUf/vAHnHLKKUin0zjzzDOr1qVQKOCoo47CY489hgsuuACHH344CoUCnnrqKaxfvx6HHHIIlFL46Ec/ioceegjf+ta38MEPfhAvvfQSLrnkEjz55JN48skn4fF4tnnbEp/73OdwzDHH4Pe//z02bNiAr3/96zjttNPw8MMPj/num266CZ/5zGdwwgkn4Ec/+hESiQQuvfRSjIyMwG7fdIZcKpVwwgkn4LHHHsM3vvENHHLIIVi3bh0uueQSLF++HP/85z/h8/n0/S+++CK++tWv4pvf/Caamppw/fXX46yzzsKiRYtw2GGHAQDefPNNHHLIIWhsbMTPfvYz1NfX47e//S3OPPNM9PT04Bvf+EbVcn/605/Gc889h8suuwyLFy9GPB7Hc889h4GBgXHb+89//jOeffZZfO9730NNTQ2uuuoqnHjiiXjzzTexYMECAMADDzyAk046CYcddhhuu+02FAoFXHPNNejp6Rn3+URXVxdOO+00fOMb38Dll18Ou92OdDqNZcuWYePGjfj2t7+NPffcE6+++iouvvhivPzyy3jwwQdhs9nw5JNP4pRTTsEpp5yCSy+9FF6vF+vWrRtXnhPFk08+icMPPxwrVqzA//t//w8AEAqFAAD/+Z//iZUrV+IHP/gBli5diuHhYbzyyisTalsDAwMDAwOD9yiUgYGBgYGBwfsWN954owKgnn32WX3tjDPOUADUnXfeqa/l83nV0NCgAKjnnntOXx8YGFAOh0P953/+Z9V3FAoFlc/n1RFHHKFOPPFEff2//uu/FAB1//33l91/9tlnKwDqxhtv1Nd23nlntXTpUpXP58vuPfbYY1VLS4sqFotV33/LLbcoAOrXv/511XseeOABBUBdddVVZddvu+02BUD96le/0teWLVumli1bpv9ftWqVAqBWrVpV9t21a9duVo+pti3lde6555a966qrrlIAVFdXV9U6FotF1draqvbZZx9VKpX09fb2duVyudTcuXP1tVtvvXWzciql1LPPPqsAqF/+8pf62ty5c5XX61Xr1q3T1zKZjIpEIurss8/W10499VTl8XjU+vXry5551FFHKb/fr+LxeNV2q6mpURdccEHVuik12rayDkopBUA1NTWpoaEhfa27u1vZ7XZ1xRVX6Gv777+/amtrUyMjI/paMplU9fX1aiLb5WXLlikA6qGHHiq7fsUVVyi73V42vpRS6o477lAA1F/+8hellFLXXHONAqDboBIo+7Vr15Zdr9T/KrVFIBBQZ5xxxmbP3X333dVHP/rRcetoYGBgYGBg8P6BcUM2MDAwMDAw2Aw2mw1HH320/t/pdGLRokVoaWnB0qVL9fVIJILGxkasW7eu7PvXXXcd9tlnH3i9XjidTrhcLjz00EN4/fXX9T2PPvoogsEgPvKRj5R99xOf+ETZ/++88w7eeOMNfOpTnwIwainIn6OPPhpdXV148803q9bl/vvvh9frLXPptIIWXFYLxZNPPhmBQAAPPfRQ1e9uKabatgBw/PHHl/2/5557AkDFe4k333wTnZ2d+OQnP1nmWjt37lwccsghZffed999CIfDOO6448rae++990Zzc/NmLq9777035syZo//3er1YvHhxWXkefvhhHHHEEWhrayv77plnnol0Oo0nn3yyatkPOOAA3HTTTfjBD36Ap556Cvl8vuq9VqxYsQLBYFD/39TUVNauw8PD+Oc//4mPfvSjcLvd+r6amhocd9xxE35PXV0dDj/88LJr9913H3bffXfsvffeZe145JFHlrkO77///gCAj3/84/jjH/+Ijo6OCb93qjjggANw//3345vf/CYeeeQRZDKZrfZuAwMDAwMDg+0Thiw0MDAwMDAw2Ax+vx9er7fsmtvtRiQS2exet9uNbDar///xj3+ML37xizjwwANx55134qmnnsKzzz6Lj3zkI2VExMDAAJqamjZ7nvUaXUG/9rWvweVylf2ce+65AID+/v6qdenr60Nra2uZm60VAwMDcDqdaGhoKLtus9nQ3Nw8rS6ZU2lbor6+vux/ukiPRfSwDs3NzZt9Zr3W09ODeDwOt9u9WZt3d3dv1t7W8rBMVnm3tLRsdl9ra2tZ+SrhtttuwxlnnIHrr78eBx98MCKRCE4//XR0d3dX/c5EyxaLxaCUmlBfHAuV6tbT04OXXnppszYMBoNQSul2POyww3D33XejUCjg9NNPx+zZs7H77rtvlbiBP/vZz3DhhRfi7rvvxooVKxC
2023-03-29 22:12:50 +02:00
"text/plain": [
2023-03-29 22:25:18 +02:00
"<Figure size 1280x240 with 12 Axes>"
2023-03-29 22:12:50 +02:00
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"image = np.asarray(Image.open(f\"./images/{image_name}_gray_noisy.jpg\")).T\n",
"\n",
"fig, axs = plot.subplots(\n",
2023-03-29 22:25:18 +02:00
" 1,\n",
" 12,\n",
2023-03-29 22:12:50 +02:00
" layout=\"constrained\",\n",
2023-03-29 22:25:18 +02:00
" figsize=(12.8, 2.4)\n",
2023-03-29 22:12:50 +02:00
")\n",
"\n",
"fig.suptitle(\"Image column denoising results\")\n",
"\n",
"ny, nx = image.shape\n",
"flat_image = image.flatten()\n",
"for i in range(len(axs.flat)):\n",
" beta = 3**i\n",
" ax = axs.flat[i]\n",
"\n",
" H = hessian(beta, nx, ny)\n",
"\n",
" smooth = sparse.linalg.spsolve(H, flat_image).reshape(image.shape).T\n",
" smooth_image = Image.fromarray(smooth.astype(np.uint8), mode=\"L\")\n",
"\n",
" ax.set_title(f\"β={beta}\")\n",
2023-03-29 22:25:18 +02:00
"\n",
" # Turn off tick labels\n",
" ax.set_yticklabels([])\n",
" ax.set_xticklabels([])\n",
"\n",
2023-03-29 22:12:50 +02:00
" ax.imshow(smooth_image)"
]
2023-03-28 01:04:30 +02:00
}
],
"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",
2023-03-28 02:35:48 +02:00
"version": "3.10.9"
2023-03-28 01:04:30 +02:00
}
},
"nbformat": 4,
"nbformat_minor": 5
}
