{ "cells": [ { "cell_type": "markdown", "id": "56889fcd-f08a-4475-9928-af817748f5c8", "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "8a937467843c19736f57c1141ce8931c", "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": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "10edd7f8d5f42ef44d458651f4e462e8", "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": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "code", "checksum": "037ff64536fbf78f10bab91b57a07726", "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": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "67dbe5349f8ff0b6f04846a0fe91ed38", "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": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "2406efd7aa00c7ce33332260d2e4fad0", "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": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "code", "checksum": "0a4c87a4860e908c141ed0281edde75c", "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": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "1aa23c460af2bb6143918dc6607b1933", "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": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "code", "checksum": "79572ce891917abec83692f4862e83f1", "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": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "6571422f2442b26e90d730b9d3cfe3b2", "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": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "cedb655df9987c5d2dea46390b184ca1", "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": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "56d954aa02fd845829faf347cab0e735", "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": { "deletable": false, "nbgrader": { "cell_type": "code", "checksum": "2e25d48b8047199adc909a8b2f6e4c09", "grade": false, "grade_id": "cell-3cfc35a0429b5bf9", "locked": false, "schema_version": 3, "solution": true, "task": false } }, "outputs": [], "source": [ "# YOUR CODE HERE\n", "raise NotImplementedError()" ] }, { "cell_type": "markdown", "id": "0d83815f-d787-470e-9653-72a0a2d02831", "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "b4385622810b4b10e3d9cdf58c4186cf", "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": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "4c223277428c0e10741952b93dd61e30", "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": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "aeb77bc687d4655332990cb21539a548", "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" ] }, { "cell_type": "markdown", "id": "a1489c37-0ce8-413c-b8ce-141ad0e49b28", "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "f03ae142788aad18c6579edcda72ddd1", "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": { "deletable": false, "nbgrader": { "cell_type": "code", "checksum": "e684fc9a01af66cc02be2a137f0baccd", "grade": false, "grade_id": "cell-6cfcc832e82f1cdd", "locked": false, "schema_version": 3, "solution": true, "task": false } }, "outputs": [], "source": [ "# YOUR CODE HERE\n", "raise NotImplementedError()" ] }, { "cell_type": "markdown", "id": "00d8690d-1187-436b-be6f-83862d6bb44b", "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "e3ae3276ea4494c0479271987dc88d30", "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": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "7b233b50cbb9077f39daf151c385b916", "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": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "330bc2a014d483fb9bac49c7bac749fb", "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": { "deletable": false, "nbgrader": { "cell_type": "code", "checksum": "ef002d074b9a41036b98da7c5e951577", "grade": false, "grade_id": "cell-c45265f255105923", "locked": false, "schema_version": 3, "solution": true, "task": false } }, "outputs": [], "source": [ "def f(x):\n", " # YOUR CODE HERE\n", " raise NotImplementedError()\n", "\n", "# YOUR CODE HERE\n", "raise NotImplementedError()" ] }, { "cell_type": "markdown", "id": "0895c0c8-328d-49b0-8c7d-26239bf1d1ef", "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "bbbe3da13e8fa92191606af6b72b4e80", "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": { "deletable": false, "nbgrader": { "cell_type": "code", "checksum": "26bd2c15dfc9a18f3727f84b11154b14", "grade": false, "grade_id": "cell-cb6ad7a07390cd65", "locked": false, "schema_version": 3, "solution": true, "task": false } }, "outputs": [], "source": [ "# YOUR CODE HERE\n", "raise NotImplementedError()" ] }, { "cell_type": "markdown", "id": "7d23c852-40cb-408b-a125-7d61ac21d9ff", "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "217ae0cb95e7609cfa6cf4f43bf0064f", "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": { "deletable": false, "nbgrader": { "cell_type": "code", "checksum": "1843c9ad5b4a1ef36fbfee6461d99863", "grade": false, "grade_id": "cell-e4ec6d9523788126", "locked": false, "schema_version": 3, "solution": true, "task": false } }, "outputs": [], "source": [ "# YOUR CODE HERE\n", "raise NotImplementedError()" ] }, { "cell_type": "markdown", "id": "c3288548-f7e6-4c9a-ac8c-ac672ecedd9b", "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "ad627e191f931987b9001a187cdbc9da", "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": { "deletable": false, "nbgrader": { "cell_type": "code", "checksum": "79afd66a2fa60af7c4ea14d4e618b340", "grade": false, "grade_id": "cell-12f0869c68e6e15e", "locked": false, "schema_version": 3, "solution": true, "task": false } }, "outputs": [], "source": [ "# YOUR CODE HERE\n", "raise NotImplementedError()" ] }, { "cell_type": "markdown", "id": "c2566fd9-0e88-418c-8d92-07c47de0286d", "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "markdown", "checksum": "a30417f3f4997b0ff1da11017c847b98", "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 }