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-06-04 15:11:18 +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-30 02:36:40 +02:00
"execution_count": 33,
2023-03-28 01:04:30 +02:00
"id": "2b9f2b79",
"metadata": {},
"outputs": [
{
"data": {
2023-03-30 02:36:40 +02:00
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAHgCAYAAABZ+0ykAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAACTQklEQVR4nOzdd3hT5dvA8W+S7pXuAXQxyt5llAIFka2ioKCogKCviKCIyFARcICoCCpDcQCKIiqg/hDZe5c9ZRbKKJQOumdy3j9CA6GDtrRNx/3xytXmnOecc+ckkrvPVCmKoiCEEEIIUUmozR2AEEIIIURJkuRGCCGEEJWKJDdCCCGEqFQkuRFCCCFEpSLJjRBCCCEqFUluhBBCCFGpSHIjhBBCiEpFkhshhBBCVCqS3AghhBCiUpHkRogHtGjRIlQqlfFhY2ODt7c3nTt3Zvr06URHRxf73CdPnmTKlClcvHix5AIuwJYtW1CpVGzZsqVMrnc/nTp1olOnTkU6pry9BnPK616sXr2aKVOmmC0mIcqCJDdClJCFCxeye/du1q9fz9y5c2nWrBkzZsygfv36bNiwoVjnPHnyJFOnTi2z5Ka8mTdvHvPmzSvSMS1atGD37t20aNGilKKq2FavXs3UqVPNHYYQpcrC3AEIUVk0atSI4OBg4/N+/frxxhtv0L59e/r27cvZs2fx8vIyY4QVT4MGDYp8jJOTE23bti2FaIonLS0NGxsbVCqVuUMRosqQmhshSpGfnx8zZ84kKSmJb775xrh9//79PP300wQEBGBra0tAQADPPPMMly5dMpZZtGgRTz31FACdO3c2NnstWrQIgPXr19OnTx9q1KiBjY0NtWvX5uWXXyYmJqZQsf3333/06NEDOzs73N3dGT58OElJSXmW3bBhA126dMHJyQk7OztCQ0PZuHGjSZkpU6agUqk4ceIEzzzzDFqtFi8vL4YOHUpCQoJJ2fT0dCZOnEhgYCBWVlZUr16dV199lVu3bpmUy6tZav78+TRt2hQHBwccHR2pV68eb7/9tnF/Xk0xQ4YMwcHBgXPnztGrVy8cHBzw9fXlzTffJCMjw+T8V65c4cknn8TR0RFnZ2eeffZZwsPDTe59fnKaKNetW8fQoUPx8PDAzs7OeI1ly5YREhKCvb09Dg4OdO/enUOHDpmc48KFCzz99NNUq1YNa2trvLy86NKlC4cPHzaWUalUeTYtBQQEMGTIkHzjGzJkCHPnzjWeI+eRUzP4+++/06ZNG7RaLXZ2dtSsWZOhQ4cW+JqFKI+k5kaIUtarVy80Gg3btm0zbrt48SJ169bl6aefxtXVlaioKObPn0+rVq04efIk7u7u9O7dm2nTpvH2228zd+5cYzNLrVq1ADh//jwhISG8+OKLaLVaLl68yOeff0779u05duwYlpaW+cZ048YNwsLCsLS0ZN68eXh5efHzzz8zcuTIXGWXLFnCoEGD6NOnD4sXL8bS0pJvvvmG7t27s3btWrp06WJSvl+/fgwYMIBhw4Zx7NgxJk6cCMAPP/wAgKIoPP7442zcuJGJEyfSoUMHjh49yuTJk9m9eze7d+/G2to6z7h//fVXRowYwahRo/jss89Qq9WcO3eOkydP3vd9yMrK4rHHHmPYsGG8+eabbNu2jQ8++ACtVst7770HQEpKCp07dyYuLo4ZM2ZQu3Zt1qxZw4ABA+57/rsNHTqU3r1789NPP5GSkoKlpSXTpk3j3Xff5YUXXuDdd98lMzOTTz/9lA4dOrBv3z5jLVWvXr3Q6XR88skn+Pn5ERMTw65du3IlfsUxadIkUlJS+OOPP9i9e7dxu4+PD7t372bAgAEMGDCAKVOmYGNjw6VLl9i0adMDX1eIMqcIIR7IwoULFUAJDw/Pt4yXl5dSv379fPdnZ2crycnJir29vfLFF18Yt//+++8KoGzevLnAGPR6vZKVlaVcunRJAZS//vqrwPLjx49XVCqVcvjwYZPtXbt2NbleSkqK4urqqjz66KMm5XQ6ndK0aVOldevWxm2TJ09WAOWTTz4xKTtixAjFxsZG0ev1iqIoypo1a/Ist2zZMgVQFixYYNwWFhamhIWFGZ+PHDlScXZ2LvC1bd68Odc9Gzx4sAIov/32m0nZXr16KXXr1jU+nzt3rgIo//77r0m5l19+WQGUhQsXFnjtnM/CoEGDTLZHRkYqFhYWyqhRo0y2JyUlKd7e3kr//v0VRVGUmJgYBVBmz55d4HUAZfLkybm2+/v7K4MHDzY+z+tevPrqq0pe//R/9tlnCqDcunWrwGsLURFIs5QQZUBRFJPnycnJjB8/ntq1a2NhYYGFhQUODg6kpKRw6tSpQp0zOjqa4cOH4+vri4WFBZaWlvj7+wPc9xybN2+mYcOGNG3a1GT7wIEDTZ7v2rWLuLg4Bg8eTHZ2tvGh1+vp0aMH4eHhpKSkmBzz2GOPmTxv0qQJ6enpxlFjOTUB9zafPPXUU9jb2+dq7rpb69atuXXrFs888wx//fVXoZvgwNAM8+ijj+aK7e6mwK1bt+Lo6EiPHj1Myj3zzDOFvg4Yaq/utnbtWrKzsxk0aJDJfbSxsSEsLMzYhObq6kqtWrX49NNP+fzzzzl06BB6vb5I1y6uVq1aAdC/f39+++03rl69WibXFaI0SHIjRClLSUkhNjaWatWqGbcNHDiQOXPm8OKLL7J27Vr27dtHeHg4Hh4epKWl3fecer2ebt26sWLFCsaNG8fGjRvZt28fe/bsAbjvOWJjY/H29s61/d5tN27cAODJJ5/E0tLS5DFjxgwURSEuLs7kGDc3N5PnOU1MOTHFxsZiYWGBh4eHSTmVSoW3tzexsbH5xv3888/zww8/cOnSJfr164enpydt2rRh/fr1Bb5eADs7O2xsbHLFlp6ebnweGxubZ6fvonYE9/HxMXmecx9btWqV6z4uW7bMmKSpVCo2btxI9+7d+eSTT2jRogUeHh689tpr+faHKikdO3bkzz//NCZhNWrUoFGjRixdurRUrytEaZA+N0KUsn/++QedTmfsGJuQkMCqVauYPHkyEyZMMJbLyMjIlSjk5/jx4xw5coRFixYxePBg4/Zz584V6ng3NzeuX7+ea/u929zd3QH46quv8h2BVNQvfjc3N7Kzs7l586ZJgqMoCtevXzfWIOTnhRde4IUXXiAlJYVt27YxefJkHnnkEc6cOWOsuSouNzc39u3bl2t7XveqIPeOjMq5j3/88cd9Y/T39+f7778H4MyZM/z2229MmTKFzMxMvv76a8CQlN3bERooMDEsjD59+tCnTx8yMjLYs2cP06dPZ+DAgQQEBBASEvJA5xaiLEnNjRClKDIykrFjx6LVann55ZcBwxefoii5Os1+99136HQ6k2331nrkyPnyvPccd4/IKkjnzp05ceIER44cMdn+yy+/mDwPDQ3F2dmZkydPEhwcnOfDysqqUNfMkdMBecmSJSbbly9fTkpKSq4Oyvmxt7enZ8+evPPOO2RmZnLixIkixZGXsLAwkpKS+Pfff022//rrrw903u7du2NhYcH58+fzvY95CQoK4t1336Vx48YcPHjQuD0gIICjR4+alN20aRPJycn3jSW/z9S9ZcLCwpgxYwZArhFdQpR3UnMjRAk5fvy4sS9FdHQ027dvZ+HChWg0GlauXGmspXBycqJjx458+umnuLu7ExAQwNatW/n+++9xdnY2OWejRo0AWLBgAY6OjtjY2BAYGEi9evWoVasWEyZMQFEUXF1d+d///leo5hmA0aNH88MPP9C7d28+/PBD42ip//77z6Scg4MDX331FYMHDyYuLo4nn3wST09Pbt68yZEjR7h58ybz588v0n3q2rUr3bt3Z/z48SQmJhIaGmocLdW8eXOef/75fI996aWXsLW1JTQ0FB8fH65fv8706dPRarX3rfEpjMGDBzNr1iyee+45PvzwQ2rXrs2///7L2rVrAVCri/f3YEBAAO+//z7vvPMOFy5coEePHri4uHDjxg327duHvb09U6dO5ejRo4wcOZKnnnqKOnXqYGVlxaZNmzh69KhJLd/zzz/PpEmTeO+
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-30 02:36:40 +02:00
"execution_count": 34,
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-30 02:36:40 +02:00
"image/png": "iVBORw0KGgoAAAANSUhEUgAABQsAAAPLCAYAAAD4+99NAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOzdd3gUVdvA4d/uphcSkkAIJSSUUAxIb6EjXaRKBF4izVdEQIpUpSuhCyhgA8IrqFGpKlKkF8HQ/GhSA6EEQgKkt92d74+YlSU9JNmU576uvUxmz5x5Zlj3yZw5RaUoioIQQgghhBBCCCGEEKLEU5s6ACGEEEIIIYQQQgghROEgjYVCCCGEEEIIIYQQQghAGguFEEIIIYQQQgghhBD/kMZCIYQQQgghhBBCCCEEII2FQgghhBBCCCGEEEKIf0hjoRBCCCGEEEIIIYQQApDGQiGEEEIIIYQQQgghxD+ksVAIIYQQQgghhBBCCAFIY6EQQgghhBBCCCGEEOIf0lgohBBCCCGEEEIIIYQApLFQCCGEEEIIIYQQQgjxD2ksFEIIIUQaAQEBqFQqw8vKyopy5crRrl07/P39CQsLy3Xdx48fZ/bs2Tx9+jTvAs7E7NmzUalUBXKsrKRe11u3buXrPsVZetejoD9TeSUuLo45c+ZQs2ZNrKyscHZ2ZsCAATx58sTUoQkhhBCiBJPGQiGEEEJkaP369fzxxx/s3buXVatWUa9ePRYuXEitWrX4/fffc1Xn8ePHmTNnTpFr2MkL3bt3548//sDNzS1f9ylpiuJnSlEUBgwYwPLly3n77bf57bffmDlzJoGBgSxdutTU4QkhhBCiBDMzdQBCCCGEKLy8vb1p1KiR4fe+ffsyfvx4WrZsSZ8+fbh27Rqurq4mjLBoKVOmDGXKlMn3ffJbXFwcNjY2pg6jSDt06BA7duzgxx9/pF+/fgCGnruxsbEmjk4IIYQQJZn0LBRCCCFEjri7u7N06VKio6P54osvALh+/TpDhw6levXq2NjYUKFCBXr06MH58+eN9p09ezaTJk0CwNPT0zDM+eDBg9muIzO//vor9erVw9LSEk9PT5YsWZJuuWvXrjFw4EDKli2LpaUltWrVYtWqVWliValUXLx4kQEDBuDg4ICrqyvDhg0jMjIyTZ1Hjx6lQ4cO2NvbY2NjQ4sWLfj111+NyqQ3hPbRo0f897//pVKlSlhaWlKmTBl8fHwMPTfT2yensW3fvp26detiaWlJlSpVWLFiRbaHZ6eWO3PmDP369aN06dJUrVo129cxO+c4ZMgQPDw8Mjx2VvFl9JnK6rgZye5nMbf1A/z444+ULl2a3r17G7YdPnyYhw8f0r59+yz3T/XJJ5+wbdu2TMtotdo023Q6XbaPIYQQQoiSRXoWCiGEECLHunXrhkaj4fDhwwDcv38fZ2dnFixYQJkyZXj8+DEbNmygadOmnD17lho1agAwYsQIHj9+zKeffsqWLVsMQ2tr167NuXPnslVHRvbt20fPnj1p3rw533//PTqdjkWLFvHw4UOjcpcuXaJFixaGRs9y5cqxe/duxo4dS3h4OLNmzTIq37dvX3x9fRk+fDjnz59n2rRpAKxbt85Q5tChQ3Ts2JG6deuydu1aLC0tWb16NT169OC7777D19c3w7gHDx7MmTNn+Pjjj/Hy8uLp06ecOXOGiIiILP8dshPbrl276NOnD61btyYwMBCtVsuSJUvSXJes9OnThzfeeIORI0cSGxubo+v4IueYlcw+U/3798/VcbP7eX6R8zp+/DhNmzY1HG/37t1MnjyZDh068Oqrr2b7/E+dOsWUKVP44Ycf6NWrl9F7Z86cYeDAgVy7do1XXnmFTZs2ER0dTb9+/Th37hwNGjTg22+/pXr16tk+nhBCCCFKAEUIIYQQ4jnr169XACUoKCjDMq6urkqtWrXSfU+r1SpJSUlK9erVlfHjxxu9t3jxYgVQgoODM40hszrS07RpU6V8+fJKfHy8YVtUVJTi5OSkPPsnT+fOnZWKFSsqkZGRRvuPHj1asbKyUh4/fqwoiqLMmjVLAZRFixYZlRs1apRiZWWl6PV6w7ZmzZopZcuWVaKjo43i9/b2VipWrGgom3pdnz13Ozs7Zdy4cRmeV3r75CS2xo0bK5UqVVISExMN26KjoxVnZ2clO38Kph5r5syZRtuzex2zc45vvvmmUrly5QyP/az0rkdGn6msjptdGX0Wc1t/fHy8YmZmpsyaNUuZO3euAiiA4u7urty5cyfHsQ0cOFAxNzdXtm7davRe3bp1lV27dikRERGKr6+v0rlzZ6Vx48bKlClTlIiICCUwMFDx8fHJcfxCCCGEKN5kGLIQQgghckVRFMPPWq2W+fPnU7t2bSwsLDAzM8PCwoJr165x+fLlbNX3InXExsYSFBREnz59sLKyMmy3t7enR48eht8TEhLYt28fvXv3xsbGBq1Wa3h169aNhIQETpw4YVT3a6+9ZvR73bp1SUhIMKwIHRsby8mTJ+nXrx92dnaGchqNhsGDB3P37l2uXLmSYexNmjQhICCAjz76iBMnTpCcnJz1xcpBbKdOnaJXr15YWFgYytnZ2Rldl+zo27ev4eecXscXOccXkdvjZvezmNv6z5w5g1arpUmTJgwaNIjdu3czZ84coqOjad26NTExMQCEh4cbrUqe3svMzIxvv/2W5ORk+vfvb+gxevfuXSpWrEjnzp1xcnJiw4YNXL9+3dBj0snJif79+wMYjieEEEIIATIMWQghhBC5EBsbS0REBHXq1AFgwoQJrFq1iilTptCmTRtKly6NWq1mxIgRxMfHZ6vOF6njyZMn6PV6ypUrl+a9Z7dFRESg1Wr59NNP+fTTT9OtKzw83Oh3Z2dno98tLS0BDDE9efIERVHSXa24fPnyhuNmJDAwkI8++oivv/6aGTNmYGdnR+/evVm0aFG655Ob2NJbhCanC9M8e345vY4vco4vIrfHze5nMbf1//nnn0BKY6OLiwtVqlShU6dOeHl5MWDAAE6cOMErr7yCvb09X331VZbnuWvXLjZv3kzPnj0NnwlFUVCr/+0XYGlpiZubGxEREej1eqP3nm34F0IIIYSQxkIhhBBC5Nivv/6KTqejbdu2AGzcuBE/Pz/mz59vVC48PBxHR8ds1fkidZQuXRqVSsWDBw/SvPfsttKlSxt6/L377rvp1uXp6ZmteJ+tU61WExoamua9+/fvA+Di4pLh/i4uLixfvpzly5cTEhLCjh07mDp1KmFhYezatStHsaQXm0qlSnd+wvSuVWaeXWgkp9cxq3O0srIiMTExTR3PN9zmVG6vbXY/i7mt/88//6RKlSoZfi5SG5ktLS0ZMWJEpuf466+/8ssvv9CvXz++++47zMxS/ryvWLEit2/fZs+ePTRv3pz58+djbm5OfHw8Y8aMYdGiRfz+++/odDrs7e0zPYYQQgghShZpLBRCCCFEjoSEhPD+++/j4ODA22+/DaQ0JKX2akv166+/cu/ePapVq2a0/fneb6lyUsfzbG1tadKkCVu2bGHx4sWGocjR0dH8/PPPhnI2Nja0a9eOs2fPUrduXaOhublla2tL06ZN2bJlC0uWLMHa2hoAvV7Pxo0bqVixIl5eXtmqy93dndGjR7Nv3z6OHTuWJ7E1atSIbdu2sWTJEsP5xsTE8Msvv+S63he5jumdo4eHB2FhYTx8+NDQ4zEpKYndu3dnq86MPlNZHTcjufks5qT+P//8M01PVEVR+Prrr/H29qZ27dqZ7v+sxYsXGxbSSW0oTD2HgIAABg0axN9//02zZs3YsmULSUlJ9O3bFzs7O2rUqMFPP/2U7WMJIYQQomSQxkIhhBBCZOjChQuGuejCwsI4cuQI69evR6PRsHXrVsqUKQPAq6++SkBAADVr1qRu3bqcPn2axYsXU7FixTR1pg5dXrFiBW+++Sbm5ubUqFEjR3WkZ968eXTp0oWOHTsyceJEdDodCxcuxNbWlsePHxvKrVixgpYtW9KqVSveeecdPDw8iI6O5vr16/z888/s378/x9fJ39+fjh070q5dO95//30
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",
2023-03-30 02:36:40 +02:00
"\\lim_{\\beta \\to 0} (I + \\beta A)s - y = s - y\n",
2023-03-28 03:23:12 +02:00
"$$\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-30 02:36:40 +02:00
"execution_count": 35,
2023-03-28 01:04:30 +02:00
"id": "dfc485d5",
2023-03-29 22:12:50 +02:00
"metadata": {
2023-03-30 02:36:40 +02:00
"lines_to_next_cell": 1,
"scrolled": false
2023-03-29 22:12:50 +02:00
},
2023-03-28 03:23:12 +02:00
"outputs": [
{
"data": {
2023-03-30 02:36:40 +02:00
"image/png": "iVBORw0KGgoAAAANSUhEUgAABQsAAAPLCAYAAAD4+99NAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOzdd3gU5drH8e+m9wSSkAKk0FFAlNBBQARERAUVLEdA0VdERUCKWLFRbAcL6jkKwhFEVMAGgoh0pIPSVMCEUAIhlHRSduf9Y8nKkgRCSSYkv8917ZXszDMz9wzL3plnnmIxDMNAREREREREREREKj0XswMQERERERERERGR8kGVhSIiIiIiIiIiIgKoslBEREREREREREROU2WhiIiIiIiIiIiIAKosFBERERERERERkdNUWSgiIiIiIiIiIiKAKgtFRERERERERETkNFUWioiIiIiIiIiICKDKQhERERERERERETlNlYUiIiIiIiIiIiICqLJQRERERERERERETlNloYiIiBRp2rRpWCwWx8vLy4vw8HA6derE+PHjSU5Ovuh9r1mzhrFjx3Ly5MnLF/A5jB07FovFUibHOp+C65qQkFCq21RkRV2Psv5MXS5ZWVm89NJLNGjQAC8vL4KDg7nnnns4ceLEZdl/RkYGQ4cOJTIyEi8vL5o2bcoXX3xxWfYtIiIiFZMqC0VEROScPv30U3799VcWL17M5MmTadq0KRMnTqRhw4b8/PPPF7XPNWvW8NJLL11xFTuXQ48ePfj111+JiIgo1W0qmyvxM2UYBvfccw+TJk3ikUce4ccff+SFF15g9uzZvPXWW5flGL1792b69Om8+OKL/PjjjzRv3px77rmHzz///LLsX0RERCoeN7MDEBERkfKtUaNGxMXFOd7fcccdDBs2jHbt2tG7d292795NWFiYiRFeWUJDQwkNDS31bUpbVlYWPj4+ZodxRVu+fDnfffcdX331FXfeeSeAo+VuZmbmJe9/wYIFLF68mM8//5x77rnHsf99+/YxcuRI+vbti6ur6yUfR0RERCoWtSwUERGRCxYVFcVbb71Feno6//nPfwDYs2cPDzzwAHXr1sXHx4fq1avTs2dPtm3b5rTt2LFjGTlyJACxsbGObs7Lli0r8T7OZf78+TRt2hRPT09iY2N58803iyy3e/du7r33XqpVq4anpycNGzZk8uTJhWK1WCzs2LGDe+65h8DAQMLCwnjwwQdJTU0ttM9Vq1bRuXNn/P398fHxoU2bNsyfP9+pTFFdaI8ePcr//d//UbNmTTw9PQkNDaVt27aOlptFbXOhsX377bc0adIET09PatWqxTvvvFPi7tkF5TZv3sydd95JlSpVqF27domvY0nOccCAAcTExBR77PPFV9xn6nzHLU5JP4sXu3+Ar776iipVqtCrVy/HshUrVnDkyBFuuOGG825f4N///jfffPNNoeXz5s3Dz8+Pu+66y2n5Aw88wKFDh1i3bl2JjyEiIiKVh1oWioiIyEW5+eabcXV1ZcWKFQAcOnSI4OBgJkyYQGhoKMePH2f69Om0bNmSLVu2UL9+fQAeeughjh8/znvvvcfcuXMdXWuvuuoqtm7dWqJ9FGfJkiXcdttttG7dmi+++AKr1crrr7/OkSNHnMrt3LmTNm3aOCo9w8PDWbRoEUOGDCElJYUXX3zRqfwdd9xB3759GThwINu2bWPMmDEATJ061VFm+fLldOnShSZNmjBlyhQ8PT354IMP6NmzJ7NmzaJv377Fxn3//fezefNmXnvtNerVq8fJkyfZvHkzx44dO++/Q0liW7hwIb179+b6669n9uzZ5Ofn8+abbxa6LufTu3dv7r77bgYNGkRmZuYFXcdLOcfzOddnqk+fPhd13JJ+ni/lvNasWUPLli0dx1u0aBGjRo2ic+fO3HLLLSU+/40bNzJ69Gi+/PJLbr/9dsfy7du307BhQ9zcnP/kb9KkiWN9mzZtSnwcERERqSQMERERkSJ8+umnBmBs2LCh2DJhYWFGw4YNi1yXn59v5ObmGnXr1jWGDRvmtO6NN94wACM+Pv6cMZxrH0Vp2bKlERkZaWRnZzuWpaWlGVWrVjXO/LOnW7duRo0aNYzU1FSn7R9//HHDy8vLOH78uGEYhvHiiy8agPH66687lRs8eLDh5eVl2Gw2x7JWrVoZ1apVM9LT053ib9SokVGjRg1H2YLreua5+/n5GUOHDi32vIra5kJia968uVGzZk0jJyfHsSw9Pd0IDg42SvLnYMGxXnjhBaflJb2OJTnH/v37G9HR0cUe+0xFXY/iPlPnO25JFfdZvNj9Z2dnG25ubsaLL75ovPzyywZgAEZUVJSxf//+C47t3nvvNdzd3Y158+Y5ltetW9fo1q1bofKHDh0yAGPcuHEXHLeIiIhUfOqGLCIiIhfNMAzH7/n5+YwbN46rrroKDw8P3Nzc8PDwYPfu3ezatatE+7uUfWRmZrJhwwZ69+6Nl5eXY7m/vz89e/Z0vD916hRLliyhV69e+Pj4kJ+f73jdfPPNnDp1irVr1zrt+9Zbb3V636RJE06dOuWYETozM5N169Zx55134ufn5yjn6urK/fffz4EDB/jzzz+Ljb1FixZMmzaNV199lbVr15KXl3f+i3UBsW3cuJHbb78dDw8PRzk/Pz+n61ISd9xxh+P3C72Ol3KOl+Jij1vSz+LF7n/z5s3k5+fTokUL7rvvPhYtWsRLL71Eeno6119/PRkZGQCkpKQ4zUpe1MvNzY3PP/+cvLw8+vTp49Ri9FxduMvLDOEiIiJSvqgbsoiIiFyUzMxMjh07RuPGjQEYPnw4kydPZvTo0XTo0IEqVarg4uLCQw89RHZ2don2eSn7OHHiBDabjfDw8ELrzlx27Ngx8vPzee+993jvvfeK3FdKSorT++DgYKf3np6eAI6YTpw4gWEYRc5WHBkZ6ThucWbPns2rr77KJ598wvPPP4+fnx+9evXi9ddfL/J8Lia2oiahudCJac48vwu9jpdyjpfiYo9b0s/ixe5//fr1gL2yMSQkhFq1atG1a1fq1avHPffcw9q1a7nxxhvx9/fn448/Pu95Lly4kDlz5nDbbbc5PhPBwcFFfu6OHz8OQNWqVc+7XxEREal8VFkoIiIiF2X+/PlYrVY6duwIwIwZM+jXrx/jxo1zKpeSkkJQUFCJ9nkp+6hSpQoWi4XDhw8XWnfmsipVqjha/D322GNF7is2NrZE8Z65TxcXF5KSkgqtO3ToEAAhISHFbh8SEsKkSZOYNGkSiYmJfPfddzz99NMkJyezcOHCC4qlqNgsFkuR4xMWda3O5cyWaBd6Hc93jl5eXuTk5BTax9kVtxfqYq9tST+LF7v/9evXU6tWrWI/FwWVzJ6enjz00EPnPMf58+fzww8/cOeddzJr1izHGIWNGzdm1qxZ5OfnO41bWDBJS6NGjc65XxEREamc1A1ZRERELlhiYiIjRowgMDCQRx55BLBXJBW0aiswf/58Dh48WGj7s1u/FbiQfZzN19eXFi1aMHfuXE6dOuVYnp6ezvfff+947+PjQ6dOndiyZQtNmjQhLi6u0Ovs1nolOXbLli2ZO3eu0znZbDZmzJhBjRo1qFevXon2FRUVxeOPP06XLl3YvHnzBcVRXGxxcXF888035ObmOpZnZGTwww8/XPR+L+U6FnWOMTExJCcnO1Vq5ubmsmjRohLFU9xn6nzHLc7FfBYvZP/r168v1BLVMAw++eQTGjVqxFVXXXXO7c/0xhtvOCbSObNSsFevXmRkZDBnzhyn8tOnTycyMtIxuYqIiIjImdSyUERERM5p+/btjrHokpOTWblyJZ9++imurq7MmzeP0NBQAG655RamTZtGgwYNaNKkCZs2beKNN96gRo0ahfZZ0HX5nXfeoX///ri7u1O/fv0L2kdRXnnlFW666Sa6dOnCU089hdVqZeLEifj6+jq6XhYct127drRv355HH32UmJgY0tPT2bNnD99//z2//PLLBV+n8ePH06VLFzp16sSIESPw8PDggw8+YPv27cyaNavY8eFSU1Pp1KkT995
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-30 02:36:40 +02:00
{
"cell_type": "markdown",
"id": "ecbcc6c7",
"metadata": {},
"source": [
"The steps our algorithm take are the following:\n",
"- Create $I + \\beta A$: O(N) complexity\n",
"- Run the tridiagonal equation solver:\n",
" - LU decomposition: O(N) complexity\n",
" - Solve $Lu = y$: O(N) complexity\n",
" - Solve $Us = u$: O(N) complexity\n",
"\n",
"This means the whole algorithm also has O(N) complexity (the constant we multiply N by to bound the number of elementary operations is the sum of the constants used to bound the number of operation of each step)."
]
},
{
"cell_type": "markdown",
"id": "6444b82a",
"metadata": {},
"source": [
"Keeping a (N + 1) x (N + 1) tridiagonal matrix in memory requires $3N + 1$ elements. When performing LU decomposition, we keep an additional $2N + 1$ elements in memory. This yields a total memory consumption of $5N + 2$ elements, which means our spatial complexity is $O(N)$ with the constant we multiply $N$ by being $c = 6$.\n",
"\n",
"If we instead used dense matrices, we would need $(N + 1)^2$ elements for the base matrix alone, with the additional $2N + 1$ still being required for the composition. This yields a spatial complexity of $N^2 + 4N + 2$, which means we have spatial copmplexity $O(N^2)$ with the constant we multiply $N^2$ by being $d = 2$."
]
},
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-06-04 15:11:18 +02:00
"execution_count": 4,
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-30 02:36:40 +02:00
"scrolled": false
2023-03-29 03:05:50 +02:00
},
"outputs": [
{
"data": {
2023-06-04 15:11:18 +02:00
"image/png": "iVBORw0KGgoAAAANSUhEUgAABQsAAAPLCAYAAAD4+99NAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOzdd3gUVRfA4d+m9wTSSChJIBBRqvQainSQplTpvQr40aWKBlAQREFFIRaa0hREijTpHaQJEhOCEEoo6aTe7481K5tsGiQsSc77PPuQnb1z58zsbji5c4tGKaUQQgghhBBCCCGEEEIUeibGDkAIIYQQQgghhBBCCPFikMZCIYQQQgghhBBCCCEEII2FQgghhBBCCCGEEEKIf0ljoRBCCCGEEEIIIYQQApDGQiGEEEIIIYQQQgghxL+ksVAIIYQQQgghhBBCCAFIY6EQQgghhBBCCCGEEOJf0lgohBBCCCGEEEIIIYQApLFQCCGEEEIIIYQQQgjxL2ksFEKIQiAwMBCNRsPJkyd127Zt28bMmTONF1Q24vD29qZv377PNR5j0mg0etci9X0LCQl57rEsXbqUwMDAZ64n7TkZ08yZM9FoNDne70U6h8KkUaNGNGrUyNhhZGjdunW88sorWFtbo9FoOHv2bJ4e78qVK7z++uu4ublhb29PlSpVWLFiRa7Vf/DgQQYOHEi1atWwtLTM8nfPkiVLeOmll7C0tMTHx4dZs2aRmJiYrtzdu3fp27cvLi4u2NjYUKdOHXbv3m2wzt9++406depgY2ODi4sLffv25e7du7l1irkiu7+XU3/fuLm5ERUVle51b29v2rZtmysxLVq0iE6dOuHj44NGo8n0e1PQ3g8hhBB5QxoLhRCikNq2bRuzZs0ydhiZxrFp0yamTZv2nCN6cbRp04YjR47g4eHx3I+dW42FL5KBAwdy5MiRHO935MgRBg4cmAcRicwsXbqUpUuXGjsMg+7du0evXr0oU6YM27dv58iRI5QrVy7Pjvf48WNat27NiRMnmD9/Pj///DMNGjRgwIAB/PLLL7lyjN27d/Pbb79RqlQp6tatm2nZ999/n7fffptOnTqxY8cOhg8fzgcffMCIESP0ysXHx9O0aVN2797N4sWL+emnn3B3d6dly5bs379fr+z+/ftp1aoV7u7u/PTTTyxevJjffvuNpk2bEh8fnyvnaAz37t1j/vz5eXqMzz//nOvXr9OkSRNcXV0zLCfvhxBCiOwyM3YAQgghCpbY2FhsbGxypa6qVavmSj35laura6Z/+KXKzWtekJUoUYISJUrkeL/atWvnQTQiI6mf55dfftnYoWTo6tWrJCYm8tZbb+Hv758rdWb2PT5+/Dh///03X3/9ta63dePGjVm1ahW7du2iTZs2z3z8adOmMWPGDAA++ugj9u3bZ7Dc/fv3mTNnDoMGDeKDDz4AtL1AExMTeffddxkzZozuvfv666+5cOEChw8fpk6dOrq4K1euzIQJEzh27Jiu3vHjx1OuXDnWr1+PmZn2TxQfHx/q1avHihUrGDZs2DOfozG0bNmSjz/+mBEjRlCsWLE8OcalS5cwMdH2AalQoUKG5eT9EEIIkV3Ss1AIIQqhvn378tlnnwHaIZapj9RhVUopli5dSpUqVbC2tqZIkSK88cYb/P3333r1NGrUiAoVKvD7779Tt25dbGxs6N+/P6Adote8eXM8PDywtramfPnyTJo0iZiYmGzH8eQw5Hv37mFhYWGwp+Gff/6JRqPhk08+0W27ffs2Q4YMoUSJElhYWOiGySUlJWV6bTp06ICXlxcpKSnpXqtVqxavvvqq7vmPP/5IrVq1cHR0xMbGhtKlS+vOPzORkZEMGjQIZ2dn7OzsaNmyJVevXk1XztBwt8yueWRkJP/73//w8fHBwsKC4sWLM2bMGL1rDpCSksKSJUt076+TkxO1a9fm559/BrTX/eLFi+zfv1/3nnh7e+fKOQH89ddf9OjRAzc3NywtLSlfvrzuc5Bq3759aDQa1qxZw9SpU/H09MTBwYHXXnuNK1eupKtzxYoVVK5cGSsrK4oWLUrHjh25fPmyXhlDw5D37NlDo0aNcHZ2xtramlKlStG5c2diY2N1ZTIaHr53716GDRuGi4sLzs7OdOrUiVu3bunVHx8fzzvvvEOxYsWwsbGhYcOGnDp1KttD7OPj45k9ezbly5fHysoKZ2dnGjduzOHDh3VlHj9+zOTJk/Xe9xEjRvDo0SO9ulKHPW7dupWqVavqvpdbt27VnVf58uWxtbWlZs2aetMWgPb7amdnx8WLF2natCm2tra4uroycuRIvesF8Nlnn9GwYUPc3NywtbWlYsWKzJ8/P90w1cw+z4aGIS9btozKlStjZ2eHvb09L730ElOmTNErc+HCBdq3b0+RIkWwsrKiSpUqfPPNN3plcvr5Snsd6tevD0DXrl3TDfv8+eefdUM37e3tadasWboeramfxdOnT/PGG29QpEgRypQpk+ExQ0NDAfQaUKOiooiKisLc3DzTeLMrtbEpK9u3b+fx48f069dPb3u/fv1QSrF582bdtk2bNuHn56drmAIwMzPjrbfe4vjx49y8eROAmzdvcuLECXr16qVrmAKoW7cu5cqVY9OmTVnGNWvWLGrVqkXRokVxcHDg1Vdf5euvv0YppVcu9Xuwfft2Xn31VaytrXnppZcMDuk+evQo9erVw8rKCk9PTyZPnmxwqHVm5syZQ1JSUp5OZZDd9+55vh9CCCHyN+lZKIQQhdC0adOIiYlh/fr1en/Epg53HTJkCIGBgYwePZp58+bx4MEDZs+eTd26dTl37hzu7u66fcLCwnjrrbeYMGECH3zwge6Plr/++ovWrVszZswYbG1t+fPPP5k3bx7Hjx9nz5492YrjSa6urrRt25ZvvvmGWbNm6f1xtHLlSiwsLOjZsyegbSisWbMmJiYmTJ8+nTJlynDkyBHmzJlDSEgIK1euzPDa9O/fn/bt27Nnzx5ee+013fY///yT48eP6xokjxw5QteuXenatSszZ87EysqK69ev684tI0opOnTowOHDh5k+fTo1atTg0KFDtGrVKtP9nmTomsfGxuLv788///zDlClTqFSpEhcvXmT69OmcP3+e3377TddQ1rdvX77//nsGDBjA7NmzsbCw4PTp07pGyU2bNvHGG2/g6OioGwZqaWmZK+d06dIl6tatS6lSpViwYAHFihVjx44djB49mvDwcF3PplRTpkyhXr16fPXVV0RGRjJx4kTatWvH5cuXMTU1BSAgIIApU6bQvXt3AgICuH//PjNnzqROnTqcOHGCsmXLGow7JCSENm3a0KBBA1asWIGTkxM3b95k+/btJCQkZNlbc+DAgbRp04bVq1dz48YNxo8fz1tvvaX3GejXrx/r1q1jwoQJNGnShEuXLtGxY0ciIyMzrRsgKSmJVq1aceDAAcaMGUOTJk1ISkri6NGjhIaGUrduXd213717N5MnT6ZBgwb88ccfzJgxgyNHjnDkyBG99+7cuXNMnjyZqVOn4ujoyKxZs+jUqROTJ09m9+7dfPDBB2g0GiZOnEjbtm0JDg7G2tpat39iYiKtW7dmyJAhTJo0icOHDzNnzhyuX7/Oli1bdOWCgoLo0aOHrgHz3LlzvP/++/z555/pGmUy+h2S1tq1axk+fDijRo3io48+wsTEhGvXrnHp0iVdmStXrlC3bl3c3Nz45JNPcHZ25vvvv6dv377cuXOHCRMm6NWZnc9XWtOmTaNmzZqMGDGCDz74gMaNG+Pg4ADA6tWr6dmzJ82bN2fNmjXEx8czf/58GjVqxO7du3WNjKk6depEt27dGDp0aLpG/bSfBdA27MTGxnLlyhWmTJlCSkoKXbp00SubkpJi8GZHWhqNJsNzzMyFCxcAqFixot52Dw8PXFxcdK+nlm3QoEG6OipVqgTAxYsXKV68uG6f1O1pyx46dCjLuEJCQhgyZAilSpUCtA19o0aN4ubNm0yfPl2v7Llz53jnnXeYNGkS7u7ufPXVVwwYMABfX18aNmwIaH9XNW3aFG9vbwIDA7GxsWHp0qW
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",
2023-06-04 15:11:18 +02:00
" 300\n",
2023-03-29 03:05:50 +02:00
" )\n",
"\n",
" alpha = 1.5 / (1 + beta * max_eigenvalue)\n",
" x = np.zeros(n + 1)\n",
"\n",
2023-03-30 02:36:40 +02:00
" iba = td.add(\n",
" td.identity(n + 1), \n",
" td.scale(beta, A)\n",
" )\n",
"\n",
2023-03-29 03:05:50 +02:00
" for _ in range(num_iterations):\n",
2023-03-30 02:36:40 +02:00
" gradient = td.multiply_vector(*iba, x) - data\n",
2023-03-29 03:05:50 +02:00
" x = x - alpha * gradient\n",
"\n",
" return x\n",
"\n",
"fig, axs = plot.subplots(\n",
" 3,\n",
" 2,\n",
" layout=\"constrained\",\n",
" figsize=(12.8, 9.6)\n",
")\n",
"\n",
"N = 100\n",
2023-06-04 15:11:18 +02:00
"beta = 1000\n",
2023-03-29 03:05:50 +02:00
"\n",
"fig.suptitle(f\"Iterative vs direct denoising comparison for β={beta} and {N=}\")\n",
"\n",
"for i in range(len(axs.flat)):\n",
" ax = axs.flat[i]\n",
" iterations = 3**(i + 1)\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-30 02:36:40 +02:00
{
"cell_type": "markdown",
"id": "5665fb36",
"metadata": {},
"source": [
"The steps our gradient descent implementation takes for $K$ iterations are the following:\n",
"- Creating $A$ and $I + βA$: O(N) complexity\n",
"- Repeat for a fixed number of iterations (in our case):\n",
" - Perform tridiagonal matrix - vector multiplication: O(N) complexity\n",
" - Perform vector - scalar division: O(N) complexity\n",
" => The process of finding the largerst eigenvalue also has O(N) complexity\n",
"- Perform gradient descent for $K$ iterations:\n",
" - Perform tridiagonal matrix-vector multiplication and vector subtraction to compute the gradient: O(N) complexity\n",
" - Update $X$ in terms of the gradient: O(N) complexity\n",
" => The proper gradient descent step has time complexity O(N*K)\n",
"\n",
"We can then conclude the whole algorithm also has time complexity O(N*K).\n",
"\n",
"I think the complexity of the direct solver is better for this exact problem. In cases where we have to work with non tridiagonal matrices, the gradient descent approach might be better."
]
},
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-30 02:36:40 +02:00
"id": "46db17e7",
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
}
