Upload New File

a0765e66 · Benjamin Berkels · 68eae05c · a0765e66
Commit a0765e66 authored 5 years ago by Benjamin Berkels
--- a/exp-357.ipynb
+++ b/exp-357.ipynb
+{
+  "cells": [
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "### IPython notebook for Example 3.5.7 from the lecture"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "import numpy as np\n",
+        "import matplotlib.pyplot as plt\n",
+        "from IPython.display import set_matplotlib_formats, display, Math\n",
+        "\n",
+        "\n",
+        "def k(s, t, d):\n",
+        "    return d / np.power(d**2 + (s-t)**2, 3/2)\n",
+        "\n",
+        "\n",
+        "d = 0.1\n",
+        "n = 80\n",
+        "h = 1/n\n",
+        "\n",
+        "display(Math(r'\\text{Create }A\\text{ from Example 1.7.1, }\\xi_j=\\begin{cases}2&1\\leq j\\leq40\\\\1&41\\leq j\\leq80\\end{cases},\\ y=A\\xi'))\n",
+        "A = np.zeros((n, n))\n",
+        "for i in range(n):\n",
+        "    for j in range(n):\n",
+        "        A[i, j] = h * k((i+0.5)*h, (j+0.5)*h, d)\n",
+        "\n",
+        "xi = np.ones(n)\n",
+        "xi[:n//2] = 2\n",
+        "y = np.matmul(A, xi)\n",
+        "\n",
+        "display(Math(r'\\tilde y= y+\\delta y'))\n",
+        "np.random.seed(0)\n",
+        "y_tilde = y + 0.02*(np.random.rand(n)-0.5)"
+      ],
+      "outputs": [],
+      "execution_count": null
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "discrepancy_tilde = np.zeros(100)\n",
+        "xi_tilde_norm = np.zeros(100)\n",
+        "discrepancy = np.zeros(100)\n",
+        "xi_norm = np.zeros(100)\n",
+        "for i in range(1, 101):\n",
+        "    alpha = 0.005*i\n",
+        "    A_Tikh = np.concatenate((A, alpha*np.eye(n)), axis=0)\n",
+        "    y_tilde_Tikh = np.concatenate((y_tilde, np.zeros(n)), axis=0)\n",
+        "    xi_tilde_alpha = np.linalg.lstsq(A_Tikh, y_tilde_Tikh, rcond=0)[0]\n",
+        "    discrepancy_tilde[i-1] = np.linalg.norm(np.matmul(A, xi_tilde_alpha)-y_tilde)\n",
+        "    xi_tilde_norm[i-1] = np.linalg.norm(xi_tilde_alpha)\n",
+        "    y_Tikh = np.concatenate((y, np.zeros(n)), axis=0)\n",
+        "    xi_alpha = np.linalg.lstsq(A_Tikh, y_Tikh, rcond=0)[0]\n",
+        "    discrepancy[i-1] = np.linalg.norm(np.matmul(A, xi_alpha)-y)\n",
+        "    xi_norm[i-1] = np.linalg.norm(xi_alpha)\n",
+        "\n",
+        "set_matplotlib_formats('svg')\n",
+        "plt.title(\"Tikhonov L-curve\")\n",
+        "plt.xlabel(r'$||A\\tilde \\xi_{\\alpha}- \\tilde y||_2^2$')\n",
+        "plt.ylabel(r'$||\\xi_{\\alpha}||_2^2$')\n",
+        "plt.plot(discrepancy_tilde**2, xi_tilde_norm**2, 'k+')\n",
+        "plt.show()\n",
+        "plt.close()"
+      ],
+      "outputs": [],
+      "execution_count": null
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "plt.title(\"Tikhonov L-curves, zoom\")\n",
+        "plt.plot(discrepancy_tilde[:50]**2, xi_tilde_norm[:50]**2, 'k+', label=r'$(\\|A \\tilde \\xi_\\alpha - \\tilde y\\|_2^2, \\|\\tilde \\xi_\\alpha\\|_2^2)$')\n",
+        "plt.plot(discrepancy[:50]**2, xi_norm[:50]**2, 'k.', label=r'$(\\|A  \\xi_\\alpha -  y\\|_2^2, \\| \\xi_\\alpha\\|_2^2)$')\n",
+        "plt.legend()\n",
+        "plt.show()\n",
+        "plt.close()"
+      ],
+      "outputs": [],
+      "execution_count": null
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "plt.title(\"Tikhonov L-curves, zoom\")\n",
+        "plt.plot(discrepancy_tilde[:50], xi_tilde_norm[:50], 'k+', label=r'$(\\|A \\tilde \\xi_\\alpha - \\tilde y\\|_2, \\|\\tilde \\xi_\\alpha\\|_2)$')\n",
+        "plt.plot(discrepancy[:50], xi_norm[:50], 'k.', label=r'$(\\|A  \\xi_\\alpha -  y\\|_2, \\| \\xi_\\alpha\\|_2)$')\n",
+        "plt.legend()\n",
+        "plt.show()\n",
+        "plt.close()"
+      ],
+      "outputs": [],
+      "execution_count": null
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "def CGNE(A, y_tilde):\n",
+        "    xi_tilde = np.zeros(n)\n",
+        "    # xi_tilde is zero, so the initial residuum is the input data.\n",
+        "    r = y_tilde.copy()\n",
+        "    d = np.matmul(A.T, r)\n",
+        "    p = d.copy()\n",
+        "    p_norm = np.linalg.norm(p)\n",
+        "    discrepancy_tilde = np.zeros(201)\n",
+        "    discrepancy_tilde[0] = np.linalg.norm(np.matmul(A, xi_tilde)-y_tilde)\n",
+        "    xi_tilde_norm = np.zeros(201)\n",
+        "    for k in range(1, 201):\n",
+        "        q = np.matmul(A, d)\n",
+        "        beta = (np.linalg.norm(p)/np.linalg.norm(q))**2\n",
+        "        xi_tilde += beta * d\n",
+        "        r += -beta*q\n",
+        "        p = np.matmul(A.T, r)\n",
+        "        p_norm_new = np.linalg.norm(p)\n",
+        "        gamma = (p_norm_new/p_norm)**2\n",
+        "        p_norm = p_norm_new\n",
+        "        d = p + gamma*d\n",
+        "        discrepancy_tilde[k] = np.linalg.norm(np.matmul(A, xi_tilde)-y_tilde)\n",
+        "        xi_tilde_norm[k] = np.linalg.norm(xi_tilde)\n",
+        "    return discrepancy_tilde, xi_tilde_norm\n",
+        "\n",
+        "\n",
+        "discrepancy_tilde, xi_tilde_norm = CGNE(A, y_tilde)\n",
+        "\n",
+        "plt.title(\"CGNE L-curve\")\n",
+        "plt.xlabel(r'$||A\\tilde \\xi_{k}- \\tilde y||_2$')\n",
+        "plt.ylabel(r'$||\\xi_{k}||_2$')\n",
+        "plt.plot(discrepancy_tilde[2:], xi_tilde_norm[2:], 'k+')\n",
+        "plt.show()\n",
+        "plt.close()\n"
+      ],
+      "outputs": [],
+      "execution_count": null
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "plt.title(\"CGNE L-curve log-log scale\")\n",
+        "plt.xlabel(r'$||A\\tilde \\xi_{k}- \\tilde y||_2$')\n",
+        "plt.ylabel(r'$||\\xi_{k}||_2$')\n",
+        "plt.loglog(discrepancy_tilde[2:], xi_tilde_norm[2:], 'k+')\n",
+        "plt.show()\n",
+        "plt.close()"
+      ],
+      "outputs": [],
+      "execution_count": null
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "discrepancy, xi_norm = CGNE(A, y)\n",
+        "\n",
+        "plt.title(\"CGNE L-curve log-log scale, zoom\")\n",
+        "plt.loglog(discrepancy_tilde[10:80], xi_tilde_norm[10:80], 'k+', label=r'$(\\|A \\tilde \\xi_\\alpha - \\tilde y\\|_2, \\|\\tilde \\xi_\\alpha\\|_2)$')\n",
+        "plt.loglog(discrepancy[10:80], xi_norm[10:80], 'k.', label=r'$(\\|A  \\xi_\\alpha -  y\\|_2, \\| \\xi_\\alpha\\|_2)$')\n",
+        "plt.legend()\n",
+        "plt.show()\n",
+        "plt.close()\n"
+      ],
+      "outputs": [],
+      "execution_count": null
+    }
+  ],
+  "metadata": {
+    "anaconda-cloud": {},
+    "kernelspec": {
+      "display_name": "Python 3",
+      "language": "python",
+      "name": "python3"
+    },
+    "language_info": {
+      "codemirror_mode": {
+        "name": "ipython",
+        "version": 3
+      },
+      "file_extension": ".py",
+      "mimetype": "text/x-python",
+      "name": "python",
+      "nbconvert_exporter": "python",
+      "pygments_lexer": "ipython3",
+      "version": "3.6.1"
+    }
+  },
+  "nbformat": 4,
+  "nbformat_minor": 1
+}
\ No newline at end of file
+%% Cell type:markdown id: tags:
+
+### IPython notebook for Example 3.5.7 from the lecture
+
+%% Cell type:code id: tags:
+
+``` python
+import numpy as np
+import matplotlib.pyplot as plt
+from IPython.display import set_matplotlib_formats, display, Math
+
+
+def k(s, t, d):
+    return d / np.power(d**2 + (s-t)**2, 3/2)
+
+
+d = 0.1
+n = 80
+h = 1/n
+
+display(Math(r'\text{Create }A\text{ from Example 1.7.1, }\xi_j=\begin{cases}2&1\leq j\leq40\\1&41\leq j\leq80\end{cases},\ y=A\xi'))
+A = np.zeros((n, n))
+for i in range(n):
+    for j in range(n):
+        A[i, j] = h * k((i+0.5)*h, (j+0.5)*h, d)
+
+xi = np.ones(n)
+xi[:n//2] = 2
+y = np.matmul(A, xi)
+
+display(Math(r'\tilde y= y+\delta y'))
+np.random.seed(0)
+y_tilde = y + 0.02*(np.random.rand(n)-0.5)
+```
+
+%% Cell type:code id: tags:
+
+``` python
+discrepancy_tilde = np.zeros(100)
+xi_tilde_norm = np.zeros(100)
+discrepancy = np.zeros(100)
+xi_norm = np.zeros(100)
+for i in range(1, 101):
+    alpha = 0.005*i
+    A_Tikh = np.concatenate((A, alpha*np.eye(n)), axis=0)
+    y_tilde_Tikh = np.concatenate((y_tilde, np.zeros(n)), axis=0)
+    xi_tilde_alpha = np.linalg.lstsq(A_Tikh, y_tilde_Tikh, rcond=0)[0]
+    discrepancy_tilde[i-1] = np.linalg.norm(np.matmul(A, xi_tilde_alpha)-y_tilde)
+    xi_tilde_norm[i-1] = np.linalg.norm(xi_tilde_alpha)
+    y_Tikh = np.concatenate((y, np.zeros(n)), axis=0)
+    xi_alpha = np.linalg.lstsq(A_Tikh, y_Tikh, rcond=0)[0]
+    discrepancy[i-1] = np.linalg.norm(np.matmul(A, xi_alpha)-y)
+    xi_norm[i-1] = np.linalg.norm(xi_alpha)
+
+set_matplotlib_formats('svg')
+plt.title("Tikhonov L-curve")
+plt.xlabel(r'$||A\tilde \xi_{\alpha}- \tilde y||_2^2$')
+plt.ylabel(r'$||\xi_{\alpha}||_2^2$')
+plt.plot(discrepancy_tilde**2, xi_tilde_norm**2, 'k+')
+plt.show()
+plt.close()
+```
+
+%% Cell type:code id: tags:
+
+``` python
+plt.title("Tikhonov L-curves, zoom")
+plt.plot(discrepancy_tilde[:50]**2, xi_tilde_norm[:50]**2, 'k+', label=r'$(\|A \tilde \xi_\alpha - \tilde y\|_2^2, \|\tilde \xi_\alpha\|_2^2)$')
+plt.plot(discrepancy[:50]**2, xi_norm[:50]**2, 'k.', label=r'$(\|A  \xi_\alpha -  y\|_2^2, \| \xi_\alpha\|_2^2)$')
+plt.legend()
+plt.show()
+plt.close()
+```
+
+%% Cell type:code id: tags:
+
+``` python
+plt.title("Tikhonov L-curves, zoom")
+plt.plot(discrepancy_tilde[:50], xi_tilde_norm[:50], 'k+', label=r'$(\|A \tilde \xi_\alpha - \tilde y\|_2, \|\tilde \xi_\alpha\|_2)$')
+plt.plot(discrepancy[:50], xi_norm[:50], 'k.', label=r'$(\|A  \xi_\alpha -  y\|_2, \| \xi_\alpha\|_2)$')
+plt.legend()
+plt.show()
+plt.close()
+```
+
+%% Cell type:code id: tags:
+
+``` python
+def CGNE(A, y_tilde):
+    xi_tilde = np.zeros(n)
+    # xi_tilde is zero, so the initial residuum is the input data.
+    r = y_tilde.copy()
+    d = np.matmul(A.T, r)
+    p = d.copy()
+    p_norm = np.linalg.norm(p)
+    discrepancy_tilde = np.zeros(201)
+    discrepancy_tilde[0] = np.linalg.norm(np.matmul(A, xi_tilde)-y_tilde)
+    xi_tilde_norm = np.zeros(201)
+    for k in range(1, 201):
+        q = np.matmul(A, d)
+        beta = (np.linalg.norm(p)/np.linalg.norm(q))**2
+        xi_tilde += beta * d
+        r += -beta*q
+        p = np.matmul(A.T, r)
+        p_norm_new = np.linalg.norm(p)
+        gamma = (p_norm_new/p_norm)**2
+        p_norm = p_norm_new
+        d = p + gamma*d
+        discrepancy_tilde[k] = np.linalg.norm(np.matmul(A, xi_tilde)-y_tilde)
+        xi_tilde_norm[k] = np.linalg.norm(xi_tilde)
+    return discrepancy_tilde, xi_tilde_norm
+
+
+discrepancy_tilde, xi_tilde_norm = CGNE(A, y_tilde)
+
+plt.title("CGNE L-curve")
+plt.xlabel(r'$||A\tilde \xi_{k}- \tilde y||_2$')
+plt.ylabel(r'$||\xi_{k}||_2$')
+plt.plot(discrepancy_tilde[2:], xi_tilde_norm[2:], 'k+')
+plt.show()
+plt.close()
+```
+
+%% Cell type:code id: tags:
+
+``` python
+plt.title("CGNE L-curve log-log scale")
+plt.xlabel(r'$||A\tilde \xi_{k}- \tilde y||_2$')
+plt.ylabel(r'$||\xi_{k}||_2$')
+plt.loglog(discrepancy_tilde[2:], xi_tilde_norm[2:], 'k+')
+plt.show()
+plt.close()
+```
+
+%% Cell type:code id: tags:
+
+``` python
+discrepancy, xi_norm = CGNE(A, y)
+
+plt.title("CGNE L-curve log-log scale, zoom")
+plt.loglog(discrepancy_tilde[10:80], xi_tilde_norm[10:80], 'k+', label=r'$(\|A \tilde \xi_\alpha - \tilde y\|_2, \|\tilde \xi_\alpha\|_2)$')
+plt.loglog(discrepancy[10:80], xi_norm[10:80], 'k.', label=r'$(\|A  \xi_\alpha -  y\|_2, \| \xi_\alpha\|_2)$')
+plt.legend()
+plt.show()
+plt.close()
+```