{
"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 03 - Linear systems and least squares "
]
},
{
"cell_type": "markdown",
"id": "e59533b3-b808-4126-b177-3390877b24ad",
"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",
"import matplotlib\n",
"import timeit\n",
"import sys"
]
},
{
"cell_type": "markdown",
"id": "c6a11316-90ea-4626-9f28-a2e771931fb6",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-84aa00fc84169745",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"
\n",
"\n",
"### Error indicators for linear systems"
]
},
{
"cell_type": "markdown",
"id": "96ada285-fd70-4e32-8922-9698937feaef",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-a8820993a4136c45",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"We consider the Vandermonde matrix\n",
"\n",
"$$\n",
"A = \\begin{bmatrix}\n",
" 1 & t_1 & t_1^2 & \\cdots & t_1^{n-1}\\\\\n",
" 1 & t_2 & t_2^2 & \\cdots & t_2^{n-1} \\\\\n",
" \\vdots & \\vdots & \\vdots & \\ddots & \\vdots \\\\\n",
" 1 & t_n & t_n^2 & \\cdots & t_n^{n-1} \n",
"\\end{bmatrix} \\in \\mathbb{R}^{n\\times n},\n",
"$$\n",
"\n",
"where the points $t_1, \\ldots, t_n$ are equispaced in the interval $[0, 1]$. For a set of points `t`, the Vandermonde matrix can be generated with `A = np.vander(t, increasing=True)`\n",
"\n",
"We want to solve the linear system $A\\mathbf{x}=\\mathbf{b}$ for $\\mathbf{x}$ where\n",
"$$\n",
"\\mathbf{b} = \\begin{bmatrix}\n",
" 1+t_1^2 \\\\\n",
" 1+t_2^2 \\\\\n",
" \\vdots \\\\\n",
" 1+t_n^2 \n",
"\\end{bmatrix} \\in \\mathbb{R}^{n}.\n",
"$$\n",
"Clearly, the exact solution is $\\mathbf{x}=(1,0,1,0\\ldots,0)^\\top \\in \\mathbb{R}^{n}$."
]
},
{
"cell_type": "markdown",
"id": "4dd7a1a8-6d06-445e-b79e-1c8a06b60168",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-c4b90ad935697c8b",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"\n",
" \n",
"**Exercise 1:** For $n=4$, solve the linear system numerically using the command `np.linalg.solve`. Let $\\mathbf{x_c}$ be the obtained solution. Compute the relative error $\\varepsilon=\\|\\mathbf{x_c}-\\mathbf{x}\\| / \\|\\mathbf{x}\\|$. \n",
"\n",
"*Hint:* The Euclidean norm of a vector can be computed with the NumPy command `np.linalg.norm`.\n",
"
\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "b39c4e6b-b1b0-49b1-af62-1c3fe8cc0de0",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-25b5edf479728cc5",
"locked": false,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"### BEGIN SOLUTION\n",
"n = 4\n",
"t = np.linspace(0, 1, n)\n",
"A = np.vander(t, increasing=True)\n",
"b = 1 + t ** 2\n",
"\n",
"x_c = np.linalg.solve(A, b)\n",
"\n",
"x_ex = np.zeros(n)\n",
"x_ex[[0, 2]] = 1\n",
"err = np.linalg.norm(x_c - x_ex) / np.linalg.norm(x_ex)\n",
"\n",
"print(r\"relative error of obtained solution: ε = {:.4e}\".format(err))\n",
"### END SOLUTION"
]
},
{
"cell_type": "markdown",
"id": "c869583c-d712-46dc-88bb-efaee81612de",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-082097f7af39bc91",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"\n",
" \n",
"**Exercise 2:** For the same linear system but for matrix sizes $n=4, 6, 8, \\ldots, 20$, visualize the relative error $\\varepsilon$ and the normalized residual $r = \\|\\mathbf{b}-A\\mathbf{x_c}\\|/\\|\\mathbf{b}\\|$ in terms of $n$, in two graphs, one in logarithmic scale (`plt.loglog`) and one in semi-logarithmic scale (`plt.semilogy`). What type of growth do we see for the error $\\varepsilon$? Is the residual $r$ a good indicator of the error $\\varepsilon$?\n",
"
"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "453aed5a-d2e3-4d05-8a32-dc45bf98a76d",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-eabe04df715c7588",
"locked": false,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"### BEGIN SOLUTION\n",
"n_list = 2 * np.arange(2, 11)\n",
"\n",
"err_list = []\n",
"res_list = []\n",
"\n",
"for n in n_list:\n",
"\n",
" x = np.linspace(0, 1, n)\n",
" A = np.vander(x, increasing=\"True\") \n",
" b = 1 + x ** 2\n",
" x_ex = np.zeros(n)\n",
" x_ex[[0, 2]] = 1\n",
" \n",
" cond_A = np.linalg.cond(A) \n",
" x_c = np.linalg.solve(A, b)\n",
" \n",
" res_list.append(np.linalg.norm(b - A @ x_c) / np.linalg.norm(b))\n",
" err_list.append(np.linalg.norm(x_c - x_ex) / np.linalg.norm(x_ex))\n",
"\n",
"plt.loglog(n_list, err_list, label=r\"relative error $\\varepsilon$\")\n",
"plt.loglog(n_list, res_list, color=\"black\", linestyle=\"--\", label=r\"residual $r$\")\n",
"plt.legend()\n",
"plt.grid(True, which=\"major\", linestyle=\"-\")\n",
"plt.grid(True, which=\"minor\", linestyle=\"--\")\n",
"plt.xlabel(r\"matrix size $n$\")\n",
"plt.show()\n",
"\n",
"plt.semilogy(n_list, err_list, label=r\"relative error $\\varepsilon$\")\n",
"plt.semilogy(n_list, res_list, color=\"black\", linestyle=\"--\", label=r\"residual $r$\")\n",
"plt.legend()\n",
"plt.grid(True, which=\"major\", linestyle=\"-\")\n",
"plt.grid(True, which=\"minor\", linestyle=\"--\")\n",
"plt.xlabel(r\"matrix size $n$\")\n",
"plt.show()\n",
"\n",
"# We notice that the relative error as well as the upper bound grow linearly in\n",
"# the graph in semi-logarithmic scale, which shows that the growth of the error\n",
"# is exponential: ε ≈ exp(αn) for some α > 0.\n",
"\n",
"# However, the residual is not at all a good indicator of the error as it is always\n",
"# on the order of 10⁻¹⁵ while the true relative error becomes quite large, of\n",
"# the order of a few percents for n = 20. This is due to the fact that the matrix\n",
"# is ill-conditioned.\n",
"### END SOLUTION"
]
},
{
"cell_type": "markdown",
"id": "3a449edc-be9a-48f3-b0ac-f73faea9de57",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-a0ae4caad8c3c188",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"We now consider a different matrix\n",
"$$\n",
"B = \\begin{bmatrix} 2 & -1 & & & \\\\ -1 & 2 & \\ddots & & \\\\ & \\ddots\n",
"& \\ddots & \\ddots & \\\\ & &\\ddots & \\ddots & -1 \\\\ &&&-1&2 \\end{bmatrix} \\in \\mathbb{R}^{n \\times n}\n",
"$$\n",
"which can be generated with the command `B = np.diag(2 * np.ones(n), 0)-np.diag(np.ones(n - 1), 1) - np.diag(np.ones(n - 1), -1)`.\n",
"\n",
"We want to solve the linear system $B \\mathbf{y}=\\mathbf{c}$ where $\\mathbf{c}=(2,2,\\cdots,2)^\\top \\in \\mathbb{R}^{n}$. The exact solution in this case is\n",
"$$\n",
"\\mathbf{y} =\n",
"\\begin{bmatrix}\n",
"1\\cdot(n) \\\\\n",
"2\\cdot(n-1)\\\\\n",
"3\\cdot(n-2)\\\\\n",
"\\vdots \\\\\n",
"n\\cdot 1 \n",
"\\end{bmatrix} \\in \\mathbb{R}^{n}.\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "7fe8c0a4-4177-4764-a252-37fe0eecb25c",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-c568133a88fa9b70",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"\n",
" \n",
"**Exercise 3:** For $n=5, 10, \\ldots, 100$, repeat the previous exercise for the matrix $B$. Comment on the result.\n",
"
"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "261e974b-2fea-4b8f-ad53-7ce4a6cbcc5e",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-06fbc99ea2e15816",
"locked": false,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"### BEGIN SOLUTION\n",
"n_list = 5 * np.arange(1, 21)\n",
"\n",
"err_list = []\n",
"res_list = []\n",
"\n",
"for n in n_list:\n",
"\n",
" y = np.linspace(0, 1, n)\n",
" B = np.diag(2 * np.ones(n), 0) - np.diag(np.ones(n - 1), 1) - np.diag(np.ones(n - 1), -1)\n",
" c = 2 * np.ones(n)\n",
" y_ex = np.arange(1, n + 1) * np.arange(n, 0, step=-1)\n",
" \n",
" cond_B = np.linalg.cond(B) \n",
" y_c = np.linalg.solve(B, c)\n",
" \n",
" res_list.append(np.linalg.norm(c - B @ y_c) / np.linalg.norm(c))\n",
" err_list.append(np.linalg.norm(y_c - y_ex) / np.linalg.norm(y_ex))\n",
"\n",
"plt.loglog(n_list, err_list, label=r\"relative error $\\varepsilon$\")\n",
"plt.loglog(n_list, res_list, color=\"black\", linestyle=\"--\", label=r\"residual $r$\")\n",
"plt.legend()\n",
"plt.grid(True, which=\"major\", linestyle=\"-\")\n",
"plt.grid(True, which=\"minor\", linestyle=\"--\")\n",
"plt.xlabel(r\"matrix size $n$\")\n",
"plt.show()\n",
"\n",
"plt.semilogy(n_list, err_list, label=r\"relative error $\\varepsilon$\")\n",
"plt.semilogy(n_list, res_list, color=\"black\", linestyle=\"--\", label=r\"residual $r$\")\n",
"plt.legend()\n",
"plt.grid(True, which=\"major\", linestyle=\"-\")\n",
"plt.grid(True, which=\"minor\", linestyle=\"--\")\n",
"plt.xlabel(r\"matrix size $n$\")\n",
"plt.show()\n",
"\n",
"# In this case, the error increases much more slowly. Therefore, we\n",
"# can conclude from the upper graph that the relative error increases as ε ≈ n².\n",
"# Moreover, the normalized residual provides this time a reasonable bound on the\n",
"# relative error, better than the bound η based on the condition number of the matrix.\n",
"### END SOLUTION"
]
},
{
"cell_type": "markdown",
"id": "4f21d81b-8a71-47b6-83bc-0748ed5bdf53",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-32ac17238c11e7c0",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"
\n",
"\n",
"### Least squares solution"
]
},
{
"cell_type": "markdown",
"id": "4b2c97a9-957d-4f05-b8f1-babf80072401",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-35894ac7feb6c7ef",
"locked": true,
"points": 0,
"schema_version": 3,
"solution": false,
"task": true
}
},
"source": [
"\n",
" \n",
"**Exercise 4 (Theoretical):** Derive the least-squares solution to $A \\mathbf{x} = \\mathbf{b}$ for a matrix $A$ with linearly independent columns. That is, show that $\\mathbf{x}^{\\ast} = (A^{\\top} A)^{-1} A^{\\top} \\mathbf{b}$ is the minimizer of the least-squares problem\n",
"$$\n",
"\\mathbf{x}^{\\ast} = \\arg\\min_{\\mathbf{x} \\in \\mathbb{R}^m} ~\\lVert A \\mathbf{x} - \\mathbf{b} \\rVert^2\n",
"$$\n",
"
\n",
"\n",
"=== BEGIN MARK SCHEME ===\n",
"You can prove this by showing that $\\mathbf{x}^{\\ast}$ is the orthogonal projection of $\\mathbf{b}$ onto the column span of $A$ as you have seen in the lecture.\n",
"\n",
"It is also possible to give an algebraic proof. For this, we write \n",
"$$\n",
"\\phi(\\mathbf{x}) = \\lVert A \\mathbf{x} - \\mathbf{b} \\rVert^2 = (A \\mathbf{x} - \\mathbf{b})^{\\top} (A \\mathbf{x} - \\mathbf{b}) = \\mathbf{x}^{\\top} A^{\\top} A \\mathbf{x} - \\mathbf{x}^{\\top} A^{\\top} \\mathbf{b} - \\mathbf{b}^{\\top} A \\mathbf{x} + \\mathbf{b}^{\\top}\\mathbf{b}.\n",
"$$\n",
"\n",
"We can replace $\\mathbf{x}^{\\top} A^{\\top} \\mathbf{b} = \\mathbf{b}^{\\top} A \\mathbf{x} $ because both are scalars (and transposing a scalar leaves it invariant). The gradient of $\\phi$ is \n",
"$$\n",
"\\nabla \\phi(\\mathbf{x}) = 2 \\mathbf{x}^{\\top} A^{\\top} A - 2 \\mathbf{b}^{\\top} A,\n",
"$$\n",
"and the Hessian is \n",
"$$\n",
"\\nabla^2 \\phi(\\mathbf{x}) = 2 A^{\\top} A.\n",
"$$\n",
"The Hessian is positive definite, since $A$ has linearly independent columns. Therefore, $\\phi$ is a strictly convex function, and assumes its unique minimum when $\\nabla \\phi(\\mathbf{x}) = 0$, i.e.,\n",
"$$\n",
"\\mathbf{x}^{\\top} A^{\\top} A = \\mathbf{b}^{\\top} \\mathbf{A} \\iff A^{\\top} A \\mathbf{x} = \\mathbf{A}^{\\top} \\mathbf{b} \\iff \\mathbf{x}^{\\ast} = ( A^{\\top} A )^{-1} \\mathbf{A}^{\\top} \\mathbf{b}\n",
"$$\n",
"\n",
"=== END MARK SCHEME ==="
]
},
{
"cell_type": "markdown",
"id": "4064c37b-f1b9-46ef-bc51-1a69c4215a37",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-f940899dd63121ec",
"locked": true,
"points": 0,
"schema_version": 3,
"solution": false,
"task": true
}
},
"source": [
"\n",
" \n",
"**Exercise 5 (Theoretical):** Compute the least-squares solution $\\mathbf{x}^{\\ast}$ for $A \\mathbf{x} = \\mathbf{b}$ \n",
"$$\n",
"A = \\begin{bmatrix} 1 & -1 \\\\ 1 & 0 \\\\ 1 & 1 \\end{bmatrix} \\quad \\text{and} \\quad \\mathbf{b} = \\begin{bmatrix} 1 \\\\ 1/2 \\\\ 1/3 \\end{bmatrix}\n",
"$$\n",
"
\n",
"\n",
"=== BEGIN MARK SCHEME ===\n",
"\n",
"We first compute\n",
"\n",
"$$\n",
"A^\\top A =\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",
"A^\\top \\mathbf{b} =\n",
"\\begin{bmatrix} \n",
" 1 & 1 & 1 \\\\ \n",
" -1 & 0 & 1 \n",
"\\end{bmatrix}\n",
"\\begin{bmatrix} \n",
" 1 \\\\ \n",
" 1/2 \\\\ \n",
" 1/3 \n",
"\\end{bmatrix}\n",
"=\n",
"\\begin{bmatrix} \n",
" 11/6 \\\\ \n",
" -2/3 \n",
"\\end{bmatrix}\n",
"$$\n",
"\n",
"Finally, the solution to \n",
"$$\n",
"\\mathbf{x}^{\\ast} = (A^{\\top} A)^{-1} A^{\\top} \\mathbf{b} = \\begin{bmatrix} \n",
" 3 & 0\\\\ \n",
" 0 & 2 \n",
"\\end{bmatrix}^{-1} \\begin{bmatrix} \n",
" 11/6 \\\\ \n",
" -2/3 \n",
"\\end{bmatrix} = \\begin{bmatrix} \n",
" 11/18 \\\\ \n",
" -1/3 \n",
"\\end{bmatrix}\n",
"$$\n",
"\n",
"=== END MARK SCHEME ==="
]
},
{
"cell_type": "markdown",
"id": "2443aa45-0fca-4113-995c-70e9adc40d89",
"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 third 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
}