{ "cells": [ { "cell_type": "markdown", "id": "56889fcd-f08a-4475-9928-af817748f5c8", "metadata": { "nbgrader": { "grade": false, "grade_id": "cell-5e2b58ff19f0997c", "locked": true, "schema_version": 3, "solution": false, "task": false } }, "source": [ "## Numerical Analysis - Fall semester 2025\n", "# Serie 07 - Conjugate gradient and curve fitting" ] }, { "cell_type": "markdown", "id": "5d0c899b-58da-49e8-ace0-828c8d5b60ca", "metadata": { "nbgrader": { "grade": false, "grade_id": "cell-e56f5ef08590e711", "locked": true, "schema_version": 3, "solution": false, "task": false } }, "source": [ "First, we will need to import some of the usual packages. You will have to run this cell every time you restart your notebook." ] }, { "cell_type": "code", "execution_count": null, "id": "5afb622e-c357-4a77-9356-e75ad0400743", "metadata": { "nbgrader": { "grade": false, "grade_id": "cell-44ae0a875bc8bfaa", "locked": true, "schema_version": 3, "solution": false, "task": false } }, "outputs": [], "source": [ "import numpy as np\n", "import scipy as sp\n", "import matplotlib.pyplot as plt\n", "np.random.seed(0); np.savetxt(\"data.txt\", np.vstack((np.linspace(0, 10, 20), np.sin(np.linspace(0, 10, 20)) + np.linspace(0, 10, 20) + 0.1 * np.random.randn(20))).T); np.savetxt(\"big_data.txt\", np.vstack((np.linspace(0, 10, 20000), np.sin(np.linspace(0, 10, 20000)) + np.linspace(0, 10, 20000) + 0.1 * np.random.randn(20000))).T)" ] }, { "cell_type": "markdown", "id": "8204bba6-b272-42e0-be57-c8a331a09ab2", "metadata": { "nbgrader": { "grade": false, "grade_id": "cell-d801e88e8b693094", "locked": true, "schema_version": 3, "solution": false, "task": false } }, "source": [ "
\n", "\n", "### Preconditioned conjugate gradient" ] }, { "cell_type": "markdown", "id": "cb042136-9122-4f9d-b6dd-309374d86365", "metadata": { "nbgrader": { "grade": false, "grade_id": "cell-046942c1a0b208cf", "locked": true, "schema_version": 3, "solution": false, "task": false } }, "source": [ "The below function generates the symmetric matrix $A \\in \\mathbb{R}^{m \\times m}$, such that the diagonal entries are $a_{jj} = 0.5 + \\sqrt{j}, j=1, 2, \\dots, m$, and the first and $\\sqrt{m}-th$ sub- and superdiagonal are $-1$. This matrix has quite a high condition number, meaning, it is *ill-conditioned*." ] }, { "cell_type": "code", "execution_count": null, "id": "97e162ba-c9aa-49ad-a5d8-27c0054a2d4c", "metadata": { "nbgrader": { "grade": false, "grade_id": "cell-d72ac70f9fc29f42", "locked": true, "schema_version": 3, "solution": false, "task": false } }, "outputs": [], "source": [ "def ill_conditioned_matrix(m):\n", " n = int(np.sqrt(m))\n", " e = np.ones(n ** 2)\n", " v = np.sqrt(np.arange(n ** 2))\n", " A = sp.sparse.spdiags([-e, -e, 0.5*e + v, -e, -e], [-n, -1, 0, 1, n])\n", " return A.toarray()\n", "\n", "A = ill_conditioned_matrix(100)\n", "plt.title(r\"$100 \\times 100$ ill-conditioned matrix\")\n", "plt.spy(A)\n", "plt.show()\n", "\n", "print(\"condition number for the 2-norm: κ(A) = {:.3f}\".format(np.linalg.cond(A)))" ] }, { "cell_type": "markdown", "id": "e52f7253-aac5-498b-914a-959350abcab5", "metadata": { "nbgrader": { "grade": false, "grade_id": "cell-39476650e6021d29", "locked": true, "schema_version": 3, "solution": false, "task": false } }, "source": [ "Since the speed with which the conjugate gradient method convergences depends on the condition number for the 2-norm $\\kappa(A) = \\lambda_{\\max}(A) / \\lambda_{\\min}(A)$, we can use a matrix $M$ such that the new condition number for the 2-norm $\\kappa(M A) = \\lambda_{\\max}(M A) / \\lambda_{\\min}(M A)$ is smaller, i.e. $M$ is a good approximation to the inverse of $A$. To measure how fast the conjugate gradient method converges for different choices of preconditioning $M$, we can use the following function `precond_conjugate_gradient_niter`, which takes as inputs a linear system with matrix `A`, a right-hand side `b`, and a preconditioner $M$, and uses the preconditioned conjugate gradient method with starting vector `x_0`, maximum number of iterations `n_max`, and tolerance `tol` to solve the system. The function returns the number of iterations it took to find a solution which satisfies the stopping criterion with tolerance `tol`." ] }, { "cell_type": "code", "execution_count": null, "id": "94e34d33-220e-43fe-81c6-7bd3d574124a", "metadata": { "nbgrader": { "grade": false, "grade_id": "cell-f23cc2ea3a078130", "locked": true, "schema_version": 3, "solution": false, "task": false } }, "outputs": [], "source": [ "def precond_conjugate_gradient_niter(A, b, M, x_0, n_max, tol):\n", " global niter\n", " niter = 0\n", " def counter_conjugate(arr):\n", " global niter\n", " niter = niter + 1\n", " sp.sparse.linalg.cg(A, b, x_0, maxiter=n_max, rtol=tol, M=M, callback=counter_conjugate)\n", " return niter" ] }, { "cell_type": "markdown", "id": "33a34512-7338-4be6-bfa4-e09a88447b03", "metadata": { "nbgrader": { "grade": false, "grade_id": "cell-1bf09cc809140a45", "locked": true, "schema_version": 3, "solution": false, "task": false } }, "source": [ "We now want to study the dependence of the number of iterations it takes the algorithm to converge for different choices of preconditioning. \n", "\n", "* No preconditioning:\n", "$$\n", "M = I_{m}\n", "$$\n", "where $I_{m}$ is the $m \\times m$ identity matrix.\n", "\n", "* Jacobi preconditioning:\n", "$$\n", "M = \\operatorname{diag}(A)^{-1},\n", "$$\n", "i.e. the matrix with the diagonal of $A$, and all other entries being zero.\n", "\n", "* Tridiagonal preconditioning:\n", "$$\n", "M = \\operatorname{tridiag}(A)^{-1},\n", "$$\n", "i.e. the matrix which just contains the diagonal and both off-diagonals of $A$, and all other entries being zero." ] }, { "cell_type": "markdown", "id": "668253ed-454f-445a-b07f-1d15b7a91f04", "metadata": { "nbgrader": { "grade": false, "grade_id": "cell-31332c77ca2b2807", "locked": true, "schema_version": 3, "solution": false, "task": false } }, "source": [ "
\n", " \n", "**Warning:** The NumPy function `np.diag` called on an $m \\times m$ matrix $A$ returns its diagonal as a length `m` vector, and not the matrix $\\operatorname{diag}(A)$ which is everywhere zero except for the diagonal entries, which coincide with the ones in $A$. However, `np.diag` called on a length $m$ vector will generate an $m \\times m$ diagonal matrix with this vector on its diagonal. Hence, you can generate $\\operatorname{diag}(A)$ by calling the `np.diag` function twice, i.e. `A_diag = np.diag(np.diag(A))`.\n", "
" ] }, { "cell_type": "markdown", "id": "314ea833-2801-4eb5-9784-22b9ce1ded95", "metadata": { "nbgrader": { "grade": false, "grade_id": "cell-fee4a2e3d7bf23e3", "locked": true, "schema_version": 3, "solution": false, "task": false } }, "source": [ "
\n", "\n", "**Exercise 1:** For the system $A \\mathbf{x} = \\boldsymbol{1}$ with $A \\in \\mathbb{R}^{100 \\times 100}$ generated with `ill_conditioned_matrix` and $\\boldsymbol{1} \\in \\mathbb{R}^{100}$ the vector of ones, use the function `precond_conjugate_gradient_niter` to determine how many iterations the preconditioned conjugate gradient method requires to satisfy the stopping criterion with tolerance `tol = 1e-16` for the different choices of the preconditioner $M$ mentioned above. Use a random normal vector $\\mathbf{x}^{(0)}$ (generated with `np.random.randn(m)`) as starting vector and choose `n_max = 100` which won't be reached before the stopping criterion is satisfied. Can you explain the difference?\n", "\n", "*Hint:* For a matrix $A$, you can extract the sub-diagonal with `np.diag(A, k=-1)` and super-diagonal with `np.diag(A, k=1)`. Same goes for generating a sub-diagonal and super-diagonal matrix from a vector.\n", "
" ] }, { "cell_type": "code", "execution_count": null, "id": "bdc5f04a-472e-4d4e-b9fa-eda335e50b22", "metadata": { "nbgrader": { "grade": false, "grade_id": "cell-3cfc35a0429b5bf9", "locked": false, "schema_version": 3, "solution": true, "task": false } }, "outputs": [], "source": [ "### BEGIN SOLUTION\n", "m = 100\n", "A = ill_conditioned_matrix(m)\n", "b = np.ones(m)\n", "x_0 = np.random.randn(m)\n", "tol = 1e-16\n", "n_max = 100\n", "\n", "M = np.eye(m)\n", "niter = precond_conjugate_gradient_niter(A, b, M, x_0, n_max, tol)\n", "print(\"no preconditioning: niter = {:d}\".format(niter))\n", "print(\"no preconditioning: κ(MA) = {:2f}\".format(np.linalg.cond(M @ A)))\n", "\n", "M = np.diag(1 / np.diag(A))\n", "niter = precond_conjugate_gradient_niter(A, b, M, x_0, n_max, tol)\n", "print(\"Jacobi preconditioning: niter = {:d}\".format(niter))\n", "print(\"Jacobi preconditioning: κ(MA) = {:2f}\".format(np.linalg.cond(M @ A)))\n", "\n", "A_tridiag = np.diag(np.diag(A)) + np.diag(np.diag(A, 1), 1) + np.diag(np.diag(A, -1), -1)\n", "M = np.linalg.inv(A_tridiag)\n", "niter = precond_conjugate_gradient_niter(A, b, M, x_0, n_max, tol)\n", "print(\"Tridiagonal preconditioning: niter = {:d}\".format(niter))\n", "print(\"Tridiagonal preconditioning: κ(MA) = {:2f}\".format(np.linalg.cond(M @ A)))\n", "\n", "# With the Jacobi preconditioning, the number of iterations needed is around 1/4\n", "# less than without preconditioning. This is reflected in the lower condition\n", "# number for the 2-norm κ(MA) for Jacobi preconditioning. With tridiagonal\n", "# preconditioning, the condition number for the 2-norm is reduced even more,\n", "# hence it converges even faster.\n", "\n", "### END SOLUTION" ] }, { "cell_type": "markdown", "id": "0d83815f-d787-470e-9653-72a0a2d02831", "metadata": { "nbgrader": { "grade": false, "grade_id": "cell-30fc4f46fd194a6d", "locked": true, "schema_version": 3, "solution": false, "task": false } }, "source": [ "
\n", "\n", "### Least squares fitting" ] }, { "cell_type": "markdown", "id": "4e0c5290-4cbe-4a99-9688-f0f1224629d6", "metadata": { "nbgrader": { "grade": false, "grade_id": "cell-469b041ade29a41c", "locked": true, "schema_version": 3, "solution": false, "task": false } }, "source": [ "Let us consider the function\n", "\n", "$$ h(x) = \\frac{x}{3+x}, \\quad x\\in [-1, 1], $$\n", "\n", "and the points $x_1=-1$, $x_2=0$ et $x_3 =1$. The polynomial $p(x) = a_0 + a_1 x$ which minimizes the function\n", "\n", "$$\n", "\\Phi(a_0, a_1) = \\sum_{i=1}^3 [h(x_i) - p(x_i)]^2= \\sum_{i=1}^3 [h(x_i) - a_0 - a_1 x_i]^2\n", "$$\n", "\n", "is the polynomial of least squares error of degree 1, also referred to as the regression line." ] }, { "cell_type": "markdown", "id": "d6449b56-6340-45b8-99b6-0ff1b4031093", "metadata": { "nbgrader": { "grade": false, "grade_id": "cell-f460765126d18da7", "locked": true, "points": 0, "schema_version": 3, "solution": false, "task": true } }, "source": [ "
\n", " \n", "**Exercise 2 (Theoretical):** Write and solve the system of equations to compute the coefficients $a_0$ and $a_1$ of the least squares polynomial $p(x) = a_0 + a_1 x$ for the data points $(x_i,h(x_i))$, $i=1,2,3$.\n", "
\n", "\n", "=== BEGIN MARK SCHEME ===\n", "\n", "The system of equations is given by\n", "\n", "$$\n", "V^\\top V\\mathbf{a}=V^\\top\\mathbf{y}\n", "$$\n", "\n", "where $\\mathbf{a}=[a_0 \\,\\, a_1]^\\top$ is the unknowns vector, $\\mathbf{y}=[h(x_1) \\,\\, h(x_2) \\,\\,h(x_3) ]^\\top$ is the data vector and $V$\n", "is the Vandermonde's matrix\n", "\n", "$$\n", "V = \n", "\\begin{bmatrix}\n", "1 & x_1 \\\\ \n", "1 & x_2 \\\\\n", "1 & x_3\n", "\\end{bmatrix}\n", "=\n", "\\begin{bmatrix}\n", "1 & -1 \\\\ \n", "1 & 0 \\\\\n", "1 & 1\n", "\\end{bmatrix}.\n", "$$\n", "\n", "We have\n", "\n", "$$\n", "V^\\top V = \n", "\\begin{bmatrix} \n", " 1 & 1 & 1 \\\\ \n", " -1 & 0 & 1\n", "\\end{bmatrix}\n", "\\begin{bmatrix} \n", " 1 & -1 \\\\ \n", " 1 & 0 \\\\ \n", " 1 & 1\n", "\\end{bmatrix}\n", "=\n", "\\begin{bmatrix} \n", " 3 & 0 \\\\ \n", " 0 & 2\n", "\\end{bmatrix}\n", "$$\n", "\n", "and\n", "\n", "$$\n", "V^\\top y=\n", "\\begin{bmatrix} \n", " 1 & 1 & 1 \\\\ \n", " -1 & 0 & 1 \n", "\\end{bmatrix}\n", "\\begin{bmatrix} \n", " -1/2 \\\\ \n", " 0 \\\\ \n", " 1/4 \n", "\\end{bmatrix}\n", "=\n", "\\begin{bmatrix} \n", " -1/4 \\\\ \n", " 3/4\n", "\\end{bmatrix}.\n", "$$\n", "\n", "The $2\\times 2$ system to solve is therefore\n", "\n", "$$\n", "\\begin{bmatrix} \n", " 3 & 0 \\\\ \n", " 0 & 2\n", "\\end{bmatrix}\n", "\\begin{bmatrix} \n", " a_0 \\\\ \n", " a_1\n", "\\end{bmatrix} =\n", "\\begin{bmatrix} \n", " -1/4 \\\\ \n", " 3/4\n", "\\end{bmatrix},\n", "$$\n", "\n", "whose solution is\n", "\n", "$$\n", "a_0 = -\\frac{1}{12}\\quad \\text{and} \\quad a_1 = \\frac{3}{8}.\n", "$$\n", "\n", "\n", "=== END MARK SCHEME ===" ] }, { "cell_type": "markdown", "id": "a1489c37-0ce8-413c-b8ce-141ad0e49b28", "metadata": { "nbgrader": { "grade": false, "grade_id": "cell-6abd888c1463b7e3", "locked": true, "schema_version": 3, "solution": false, "task": false } }, "source": [ "
\n", " \n", "**Exercise 3:** Visualize the function $h$ and the least squares polynomial $p$ you have found in Exercise 2 by plotting their values at $100$ uniformly spaced points on the interval $[-1, 1]$.\n", "
" ] }, { "cell_type": "code", "execution_count": null, "id": "f382012a-01fc-4bd7-aa89-b2de6315c67f", "metadata": { "nbgrader": { "grade": false, "grade_id": "cell-6cfcc832e82f1cdd", "locked": false, "schema_version": 3, "solution": true, "task": false } }, "outputs": [], "source": [ "### BEGIN SOLUTION\n", "\n", "def h(x):\n", " return x / (3 + x)\n", "\n", "def p(x):\n", " return - 1 / 12 + 3 / 8 * x\n", "\n", "x_lin = np.linspace(-1, 1)\n", "plt.plot(x_lin, h(x_lin), label=r\"$h(x)$\")\n", "plt.plot(x_lin, p(x_lin), label=r\"$p(x)$\")\n", "plt.xlabel(r\"$x$\")\n", "plt.legend()\n", "plt.grid()\n", "plt.show()\n", "\n", "### END SOLUTION" ] }, { "cell_type": "markdown", "id": "00d8690d-1187-436b-be6f-83862d6bb44b", "metadata": { "nbgrader": { "grade": false, "grade_id": "cell-dcb85adf51dc30cd", "locked": true, "schema_version": 3, "solution": false, "task": false } }, "source": [ "
\n", "\n", "### Least squares approximation of data" ] }, { "cell_type": "markdown", "id": "73a664d6-35bd-4981-a0b6-1c6af72e28e6", "metadata": { "nbgrader": { "grade": false, "grade_id": "cell-5924054146106def", "locked": true, "schema_version": 3, "solution": false, "task": false } }, "source": [ "Let us consider the data $(x_i,y_i)$, $i=1,\\ldots,20$ contained in the file `data.txt` which has been generated for you in the same directory as this notebook is stored in. The data corresponds to the evaluation of the function \n", "\n", "$$\n", " f(x) = \\sin(x) + x, \\quad x\\in I=[0, 10],\n", "$$\n", "\n", "on $n+1=20$ equally distributed nodes (therefore $x_1=0$ and $x_{20}=10$).\n", "The evaluations are perturbed by Gaussian noise as follows:\n", "\n", "$$\n", "y_i=f(x_i) + \\varepsilon_i, \\quad i=1,\\ldots,n+1.\n", "$$\n", "\n", "Let us suppose that the terms $\\varepsilon_i$ are stochastically independent, with zero mean and standard deviation $\\sigma=0.1$." ] }, { "cell_type": "markdown", "id": "b5357719-5686-47b4-afac-82656c7c8974", "metadata": { "nbgrader": { "grade": false, "grade_id": "cell-cbc6cd8f874059f2", "locked": true, "schema_version": 3, "solution": false, "task": false } }, "source": [ "
\n", " \n", "**Exercise 4:** Load the data using the function `np.loadtxt('[FILENAME HERE].txt')` and visualize it in a plot. Also add a visualization of the function $f$ in your plot.\n", "
" ] }, { "cell_type": "code", "execution_count": null, "id": "5f1af077-1888-454b-bb62-18e3f1015b1e", "metadata": { "nbgrader": { "grade": false, "grade_id": "cell-c45265f255105923", "locked": false, "schema_version": 3, "solution": true, "task": false } }, "outputs": [], "source": [ "def f(x):\n", " ### BEGIN SOLUTION\n", " return np.sin(x) + x\n", " ### END SOLUTION\n", "\n", "### BEGIN SOLUTION\n", "data = np.loadtxt(\"data.txt\")\n", "x = data[:, 0]\n", "y = data[:, 1]\n", "x_lin = np.linspace(0, 10, 100)\n", "plt.scatter(x, y, label=r\"data $(x_i, y_i)$\")\n", "plt.plot(x_lin, f(x_lin), color=\"black\", label=r\"function $f(x)$\")\n", "plt.legend()\n", "plt.show()\n", "### END SOLUTION" ] }, { "cell_type": "markdown", "id": "0895c0c8-328d-49b0-8c7d-26239bf1d1ef", "metadata": { "nbgrader": { "grade": false, "grade_id": "cell-63e41a2916c6a825", "locked": true, "schema_version": 3, "solution": false, "task": false } }, "source": [ "
\n", " \n", "**Exercise 5:** For $m = 4, 7, 15$, compute and plot the least-squares polynomial $p_{m}^{LS}$ of degree $m$ of the data. Do this in a separate figure for each $m$ where you also plot the function and the data. What do you observe?\n", "
" ] }, { "cell_type": "code", "execution_count": null, "id": "66874b6b-339b-4cc9-ba47-7b9e887bedaa", "metadata": { "nbgrader": { "grade": false, "grade_id": "cell-cb6ad7a07390cd65", "locked": false, "schema_version": 3, "solution": true, "task": false } }, "outputs": [], "source": [ "### BEGIN SOLUTION\n", "m_list = [4, 7, 15]\n", "for m in m_list:\n", " coeff = np.polyfit(x, y, m)\n", " y_leastsquares = np.polyval(coeff, x_lin)\n", " plt.scatter(x, y, label=r\"data $(x_i, y_i)$\")\n", " plt.plot(x_lin, f(x_lin), color=\"black\", label=r\"function $f(x)$\")\n", " plt.plot(x_lin, y_leastsquares, color=\"red\", label=f\"least squares polynomial ($m = {m}$)\")\n", " plt.legend()\n", " plt.show()\n", "\n", "# We notice that when the degree increases, the precision of the interpolation\n", "# decreases. In fact, the least squares polynomial of degree m=7 gives in this\n", "# case a better approximation than the one of degree m=15.\n", "### END SOLUTION" ] }, { "cell_type": "markdown", "id": "7d23c852-40cb-408b-a125-7d61ac21d9ff", "metadata": { "nbgrader": { "grade": false, "grade_id": "cell-0aa7223141fd86a9", "locked": true, "schema_version": 3, "solution": false, "task": false } }, "source": [ "
\n", " \n", "**Exercise 6:** Evaluate the approximation error\n", "\n", "$$\n", "E(p_m^{LS},f)=\\max_{x \\in [0,10]} \\vert f(x) - p_m^{LS}(x) \\vert\n", "$$\n", "\n", "for $m=1,\\ldots,15$. (Do this on a very fine grid.) Then plot the error in terms of $m$ in semi-logarithmic scale. What behavior do you observe for high values of $m$? For which value of $m$ is the approximation error the smallest?\n", "
" ] }, { "cell_type": "code", "execution_count": null, "id": "1b6200f9-9d8a-42a9-a44d-1b7ac6a48627", "metadata": { "nbgrader": { "grade": false, "grade_id": "cell-e4ec6d9523788126", "locked": false, "schema_version": 3, "solution": true, "task": false } }, "outputs": [], "source": [ "### BEGIN SOLUTION\n", "errors = []\n", "m_list = np.arange(1, 16)\n", "for m in m_list:\n", " coeff = np.polyfit(x, y, m)\n", " x_lin_fine = np.linspace(0, 10, 10000)\n", " y_leastsquares = np.polyval(coeff, x_lin_fine)\n", " y_true = f(x_lin_fine)\n", " errors.append(np.max(np.abs(y_true - y_leastsquares)))\n", "\n", "plt.semilogy(m_list, errors)\n", "plt.xlabel(r\"$m$\")\n", "plt.ylabel(r\"approximation error\")\n", "plt.grid(True, which=\"major\",ls=\"-\")\n", "plt.grid(True, which=\"minor\",ls=\"--\")\n", "plt.ylim(0.09, 1.5)\n", "plt.show()\n", "\n", "# We notice that the polynomial of degree m=7 gives a better approximation than\n", "# m=4. However, the polynomial of degree m=15 gives a bad approximation at the\n", "# extremities of the interval. This is due to the fact that the least squares\n", "# approximation becomes more and more unstable (for equally distributed nodes)\n", "# when m approaches n.\n", "\n", "# We can notice that the approximation error decreases up to the value m=6, \n", "# but afterwards the error increases as the number of data points (n=19) is\n", "# not big enough.\n", "### END SOLUTION" ] }, { "cell_type": "markdown", "id": "c3288548-f7e6-4c9a-ac8c-ac672ecedd9b", "metadata": { "nbgrader": { "grade": false, "grade_id": "cell-156bcd0b077ab985", "locked": true, "schema_version": 3, "solution": false, "task": false } }, "source": [ "
\n", " \n", "**Exercise 7:** Repeat all of the above exercises with the data contained in the file `big_data.txt`, which includes $n+1=20000$ data points following the same model as the data in `data.txt`. Compare for which $m$ the smallest approximation error is obtained for the two data sets and how big it is respectively.\n", "
" ] }, { "cell_type": "code", "execution_count": null, "id": "d14f7c0d-d042-4657-a2c9-ceff0d07447f", "metadata": { "nbgrader": { "grade": false, "grade_id": "cell-12f0869c68e6e15e", "locked": false, "schema_version": 3, "solution": true, "task": false } }, "outputs": [], "source": [ "### BEGIN SOLUTION\n", "def f(x):\n", " return np.sin(x) + x\n", "\n", "# Plot data\n", "data = np.loadtxt(\"big_data.txt\")\n", "x = data[:, 0]\n", "y = data[:, 1]\n", "x_lin = np.linspace(0, 10, 100)\n", "plt.scatter(x, y, label=r\"data $(x_i, y_i)$\")\n", "plt.plot(x_lin, f(x_lin), color=\"black\", label=r\"function $f(x)$\")\n", "plt.legend()\n", "plt.show()\n", "\n", "# Plot least squares approximants\n", "m_list = [4, 7, 15]\n", "for m in m_list:\n", " coeff = np.polyfit(x, y, m)\n", " y_leastsquares = np.polyval(coeff, x_lin)\n", " plt.scatter(x, y, label=r\"data $(x_i, y_i)$\")\n", " plt.plot(x_lin, f(x_lin), color=\"black\", label=r\"function $f(x)$\")\n", " plt.plot(x_lin, y_leastsquares, color=\"red\", label=f\"least squares polynomial ($m = {m}$)\")\n", " plt.legend()\n", " plt.show()\n", "\n", "# Plot errors\n", "errors = []\n", "m_list = np.arange(1, 16)\n", "for m in m_list:\n", " coeff = np.polyfit(x, y, m)\n", " x_lin_fine = np.linspace(0, 10, 1000)\n", " y_leastsquares = np.polyval(coeff, x_lin_fine)\n", " y_true = f(x_lin_fine)\n", " errors.append(np.max(np.abs(y_true - y_leastsquares)))\n", "\n", "plt.semilogy(m_list, errors)\n", "plt.xlabel(r\"$m$\")\n", "plt.ylabel(r\"interpolation error\")\n", "plt.grid(True, which=\"major\",ls=\"-\")\n", "plt.grid(True, which=\"minor\",ls=\"--\")\n", "plt.show()\n", "\n", "# Compute variance/noise estimators\n", "n = len(x) - 1\n", "m = 4\n", "coeff = np.polyfit(x, y, m)\n", "y_leastsquares = np.polyval(coeff, x)\n", "\n", "# With n+1=20000 nodes, the least squares polynomials of degree m ≤ 15 are\n", "# accurate since the number of data is very high. Therefore, the approximation\n", "# error decreases still for larger values of m.\n", "### END SOLUTION" ] }, { "cell_type": "markdown", "id": "c2566fd9-0e88-418c-8d92-07c47de0286d", "metadata": { "nbgrader": { "grade": false, "grade_id": "cell-a26fae9678ff662c", "locked": true, "schema_version": 3, "solution": false, "task": false } }, "source": [ "
\n", "\n", "## The end\n", "\n", "Congratulations! You have made it to the end of the seventh exercise notebook. " ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.12.3" } }, "nbformat": 4, "nbformat_minor": 5 }