{
"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 02 - Simple linear systems"
]
},
{
"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": "292b3b9c-9456-4645-8d7c-0ec4f5b0ba0d",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-8dd6436372d2524c",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"
\n",
"\n",
"### Simple linear systems"
]
},
{
"cell_type": "markdown",
"id": "6aeac6d6-19da-4544-a851-80aff6c469b2",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-f9b1603e57818fd9",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"Consider the following two linear systems:\n",
"\n",
"$$\n",
" A_1\\mathbf{x} = \\mathbf{b}_1, \\qquad \\text{and} \\qquad\n",
" A_1= \\begin{bmatrix}\n",
" 2 & 4 & 8 \\\\\n",
" 1 & 1 & 4 \\\\\n",
" 3 & 6 & 7\n",
" \\end{bmatrix}, \\quad\n",
" \\mathbf{b}_1 = \\begin{bmatrix} 6 \\\\ 5 \\\\ 4 \\end{bmatrix},\n",
"$$\n",
"\n",
"and\n",
"\n",
"$$\n",
" A_2\\mathbf{x} = \\mathbf{b}_2, \\qquad \\text{and} \\qquad\n",
" A_2 = \\begin{bmatrix}\n",
" 1 & 1 & 1 & 1 \\\\\n",
" 2 & 2 & 5 & 3 \\\\\n",
" 4 & 6 & 8 & 0 \\\\\n",
" 3 & 3 & 9 & 8\n",
" \\end{bmatrix}, \\quad\n",
" \\mathbf{b}_2 = \\begin{bmatrix} 1 \\\\ 2 \\\\ 5 \\\\ 0 \\end{bmatrix}.\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "82bd30d7-c89d-45eb-8d16-79c87724c63c",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-28b47eb60c1c806f",
"locked": true,
"points": 0,
"schema_version": 3,
"solution": false,
"task": true
}
},
"source": [
"\n",
" \n",
"**Exercise 1 (Theoretical):** Solve $A_1 \\mathbf{x} = \\mathbf{b}_1$ using Gaussian elimination and determine the LU factorization $A_1 = L_1 U_1$.\n",
"
\n",
"\n",
"=== BEGIN MARK SCHEME ===\n",
"\n",
"First, we append $\\mathbf{b}_1$ to $A_1$:\n",
"\n",
"$$\n",
"[A_1 \\mid \\mathbf{b}_1] = \\left[\n",
"\\begin{array}{ccc|c}\n",
"2 & 4 & 8 & 6 \\\\\n",
"1 & 1 & 4 & 5 \\\\\n",
"3 & 6 & 7 & 4\n",
"\\end{array}\n",
"\\right]\n",
"= \n",
"\\begin{bmatrix}\n",
"r_1^{(1)} \\\\\n",
"r_2^{(1)} \\\\\n",
"r_3^{(1)}\n",
"\\end{bmatrix}\n",
"$$\n",
"\n",
"To this augmented matrix, we then apply Gaussian elimination to convert it to row echelon form.\n",
"\n",
"$$\n",
"r_2^{(2)} \\leftarrow r_2^{(1)} - \\frac{1}{2} r_1^{(1)}, r_3^{(2)} \\leftarrow r_3^{(1)} - \\frac{3}{2} r_1^{(1)} \\Rightarrow\n",
"\\left[\n",
"\\begin{array}{ccc|c}\n",
"2 & 4 & 8 & 6 \\\\\n",
"0 & -1 & 0 & 2 \\\\\n",
"0 & 0 & -5 & -5\n",
"\\end{array}\n",
"\\right]\n",
"$$\n",
"\n",
"From the first three columns of the row echelon form we can read the factor $U_1$. \n",
"The entries in the $L_1$ factor are the multipliers used in the Gaussian elimination process:\n",
"\n",
"$$\n",
"l_{21}^{(1)} = \\frac{1}{2}, \\quad\n",
"l_{31}^{(1)} = \\frac{3}{2} , \\quad \\textrm{and} \\quad\n",
"l_{32}^{(1)} = 0.\n",
"$$\n",
"\n",
"Consequently, we have \n",
"\n",
"$$\n",
" L_1 = L_1^{(2)} = \\begin{bmatrix}\n",
" 1 & 0 & 0 \\\\\n",
" 1/2 & 1 & 0 \\\\\n",
" 3/2 & 0 & 1\n",
" \\end{bmatrix} \\quad \\text{and} \\quad\n",
" U_1 = U_1^{(2)} = \\begin{bmatrix}\n",
" 2 & 4 & 8 \\\\\n",
" 0 & -1 & 0 \\\\\n",
" 0 & 0 & -5\n",
" \\end{bmatrix}.\n",
"$$\n",
"\n",
"The reduced row echelon form can be obtained with further steps of Gaussian elimination:\n",
"\n",
"$$\n",
"r_1^{(3)} \\leftarrow r_1^{(2)} + 4 r_2^{(2)}, r_1^{(3)} \\leftarrow r_1^{(2)} + \\frac{8}{5} r_1^{(2)} \\Rightarrow\n",
"\\left[\n",
"\\begin{array}{ccc|c}\n",
"2 & 0 & 0 & 6 \\\\\n",
"0 & -1 & 0 & 2 \\\\\n",
"0 & 0 & -5 & -5\n",
"\\end{array}\n",
"\\right]\n",
"$$\n",
"\n",
"and \n",
"\n",
"$$\n",
"r_1^{(4)} \\leftarrow \\frac{1}{2} r_1^{(3)}, r_2^{(4)} \\leftarrow -r_2^{(3)}, r_3^{(4)} = -\\frac{1}{5} r_3^{(3)} \\Rightarrow\n",
"\\left[\n",
"\\begin{array}{ccc|c}\n",
"1 & 0 & 0 & 3 \\\\\n",
"0 & 1 & 0 & -2 \\\\\n",
"0 & 0 & 1 & 1\n",
"\\end{array}\n",
"\\right].\n",
"$$\n",
"\n",
"From this we see that the solution is \n",
"$$\n",
"\\mathbf{x} = \\begin{bmatrix} 3 \\\\ -2 \\\\ 1 \\end{bmatrix}.\n",
"$$\n",
"\n",
"=== END MARK SCHEME ==="
]
},
{
"cell_type": "markdown",
"id": "3cbfe1a7-8e5d-412f-a287-6a43281927e6",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-3b6a1d75bfd8b314",
"locked": true,
"points": 0,
"schema_version": 3,
"solution": false,
"task": true
}
},
"source": [
"\n",
" \n",
"**Exercise 2 (Theoretical):** Solve $A_2 \\mathbf{x} = \\mathbf{b}_2$ using Gaussian elimination with pivoting and determine the LU factorization of $P_2 A_2 = L_2 U_2$.\n",
"
\n",
"\n",
"=== BEGIN MARK SCHEME ===\n",
"\n",
"First, we append $\\mathbf{b}_2$ to $A_2$:\n",
"\n",
"$$\n",
"[A_2 \\mid \\mathbf{b}_2] = \\left[\n",
"\\begin{array}{cccc|c}\n",
" 1 & 1 & 1 & 1 & 1\\\\\n",
" 2 & 2 & 5 & 3 & 2\\\\\n",
" 4 & 6 & 8 & 0 & 5\\\\\n",
" 3 & 3 & 9 & 8 & 0\n",
" \\end{array}\n",
" \\right] = \n",
" \\begin{bmatrix}\n",
"r_1^{(1)} \\\\\n",
"r_2^{(1)} \\\\\n",
"r_3^{(1)} \\\\\n",
"r_4^{(1)}\n",
"\\end{bmatrix}\n",
"$$\n",
"\n",
"To this augmented matrix, we then apply Gaussian elimination to convert it to row echelon form.\n",
"\n",
"$$\n",
"r_2^{(2)} \\leftarrow r_2^{(1)} - 2 r_1^{(1)}, r_3^{(2)} \\leftarrow r_3^{(1)} - 4 r_1^{(1)}, r_4^{(2)} \\leftarrow r_4^{(1)} - 3 r_1^{(1)} \\Rightarrow\n",
" \\left[\n",
"\\begin{array}{cccc|c}\n",
" 1 & 1 & 1 & 1 & 1\\\\\n",
" 0 & \\boxed{0} & 3 & 1 & 0\\\\\n",
" 0 & 2 & 4 & -4 & 1\\\\\n",
" 0 & 0 & 6 & 5 & -3\n",
" \\end{array}\n",
" \\right]\n",
"$$\n",
"\n",
"Since the pivot is $0$, we cannot continue the algorithm. We have to swap two rows:\n",
"\n",
"$$\n",
"r_2^{(3)} \\leftrightarrow r_3^{(2)}, r_3^{(3)} \\leftrightarrow r_2^{(2)} \\Rightarrow\n",
" \\left[\n",
"\\begin{array}{cccc|c}\n",
" 1 & 1 & 1 & 1 & 1\\\\\n",
" 0 & 2 & 4 & -4 & 1\\\\\n",
" 0 & 0 & 3 & 1 & 0\\\\\n",
" 0 & 0 & 6 & 5 & -3\n",
" \\end{array}\n",
" \\right].\n",
"$$\n",
"\n",
"Swapping these two rows is equivalent to multiplying the matrix $A_2$ from the left with the permutation matrix\n",
"$$\n",
"P_2=\\begin{bmatrix}\n",
" 1 & 0 & 0 & 0 \\\\\n",
" 0 & 0 & 1 & 0 \\\\\n",
" 0 & 1 & 0 & 0 \\\\\n",
" 0 & 0 & 0 & 1\n",
" \\end{bmatrix}.\n",
"$$\n",
"and we also need to swap the multipliers in the $L_2$-factor to account for this permutation. Then we can continue with the usual row eliminations:\n",
"\n",
"$$\n",
"r_4^{(4)} \\leftrightarrow r_4^{(3)} - 2 r_3^{(3)} \\Rightarrow\n",
" \\left[\n",
"\\begin{array}{cccc|c}\n",
" 1 & 1 & 1 & 1 & 1\\\\\n",
" 0 & 2 & 4 & -4 & 1\\\\\n",
" 0 & 0 & 3 & 1 & 0\\\\\n",
" 0 & 0 & 0 & 3 & -3\n",
" \\end{array}\n",
" \\right]\n",
"$$\n",
"\n",
"Thus, we get to the following factors:\n",
"\n",
"$$\n",
" L_2 = L_2^{(4)} = \\begin{bmatrix}\n",
" 1 & 0& 0& 0 \\\\\n",
" 4 & 1 & 0& 0\\\\\n",
" 2 & 0 & 1 & 0\\\\\n",
" 3 & 0 & 2 & 1\n",
" \\end{bmatrix}, \\qquad\n",
" U_2 = U_2^{(4)} = \\begin{bmatrix}\n",
" 1 & 1 & 1 & 1 \\\\\n",
" 0 & 2 & 4 &-4 \\\\\n",
" 0 & 0 & 3 & 1 \\\\\n",
" 0 & 0 & 0 & 3\n",
" \\end{bmatrix}.\n",
"$$\n",
"\n",
"Continuing the elimination to reduced row echelon form as in Exercise 1 gives\n",
"$$\n",
"\\mathbf{x} = \\begin{bmatrix} \\frac{23}{6} \\\\ \\frac{13}{6} \\\\ \\frac{1}{3} \\\\ -1 \\end{bmatrix}.\n",
"$$\n",
"\n",
"=== END MARK SCHEME ==="
]
},
{
"cell_type": "markdown",
"id": "dfddca72-a5a7-48e6-9cfd-9378d43a487c",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-218f10fd8ec74a7d",
"locked": true,
"points": 0,
"schema_version": 3,
"solution": false,
"task": true
}
},
"source": [
"\n",
" \n",
"**Exercise 3 (Theoretical):** Compute the determinant of the matrix $A_2$ using its LU factorization.\n",
"\n",
"*Hint:* We know that\n",
"$$\n",
"\\text{det}(A_2)=\\text{det}(P_2^{-1}L_2U_2)=\\displaystyle\\frac{\\text{det}(L_2) \\text{det}(U_2)}{\\text{det}(P_2)}.\n",
"$$\n",
"
\n",
"\n",
"=== BEGIN MARK SCHEME ===\n",
"\n",
"We use the relation $\\text{det}(A_2)=\\text{det}(P_2^{-1}L_2U_2)=\\displaystyle\\dfrac{1}{\\text{det}(P_2)}\\text{det}(L_2) \\text{det}(U_2)$. Noticing \n",
"that the determinant of a triangular matrix is given by \n",
"the product of its diagonal elements, we find that\n",
"\n",
"$$\n",
" \\textrm{det}(A_2) = -1 \\cdot 1 \\cdot 18 = -18 .\n",
"$$\n",
"\n",
"Additionally, note that the determinant of a permutation matrix between two\n",
"rows is always $-1$, and therefore the determinant of $P_2$ is $+1$\n",
"if we performed an even number of permutations and\n",
"$-1$ if we performed an odd number of permutations:\n",
"in the case of this exercise, we have $\\text{det}(P_2)=-1$.\n",
"\n",
"=== END MARK SCHEME ==="
]
},
{
"cell_type": "markdown",
"id": "81d2b00c-45cd-4d89-ab1c-cc4e306b666f",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-b9890a011292fd00",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"
\n",
"\n",
"### LU factorization for dense and sparse matrix "
]
},
{
"cell_type": "markdown",
"id": "34866e1d-3d2e-488e-8db1-be7314e6fb54",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-ac16c788bf6da4d7",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"\n",
" \n",
"**Exercise 4:** Complete the function `binomial_matrix` which creates a matrix $A \\in \\mathbb{R}^{m \\times m}$ whose entries $a_{ij}$ for $i,j = 1, 2, \\dots, m$ independently follow a binomial distribution with $n \\in \\mathbb{N}$ trials and a success probability of $p \\in [0, 1]$. Use your implementation to create the $(n=10, p=0.5)$-binomial matrix $A \\in \\mathbb{R}^{400 \\times 400}$. Visualize its non-zero entries using the function `plt.spy(A)`, which plots the non-zero entries of a matrix in black.\n",
"\n",
"*Hint:* The NumPy function `np.random.binomial(n, p, (m, m))` outputs an $m \\times m$ NumPy array whose entries are independently $(n, p)$-binomial distributed. The seed of NumPy's random number generator is fixed to `seed=0`, such that you will always get the same random matrix for the same parameters.\n",
"
"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "d240c37e-6e7a-456f-bf74-7402eb39c146",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-69ad0d0c4384592d",
"locked": false,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"def binomial_matrix(m, n, p, seed=0):\n",
" np.random.seed(seed)\n",
" ### BEGIN SOLUTION\n",
" A = np.random.binomial(n, p, (m, m))\n",
" return A\n",
" ### END SOLUTION\n",
"\n",
"### BEGIN SOLUTION\n",
"n = 10\n",
"p = 0.5\n",
"A = binomial_matrix(400, n, p)\n",
"plt.spy(A)\n",
"plt.show()\n",
"### END SOLUTION"
]
},
{
"cell_type": "markdown",
"id": "f88dbe2d-9085-483a-ab64-8150e87d34cf",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-bce9533bc6dcdf29",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"\n",
" \n",
"**Exercise 5:** Compute the LU factorization using the function `sp.linalg.lu`, which takes as an input a matrix $A$ and returns the three matrices $P$, $L$, and $U$ which characterize the LU factorization of $A$. Visualize the non-zero entries of $L$, $U$, and $P$ (again use `plt.spy`). Have row-permutations been performed? How do the number of non-zero elements in $L$ and $U$ compare with what we saw for $A$?\n",
"
"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "fc19547c-088f-4e1b-816f-ecbb47f71a26",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-467b7459d7bcdec0",
"locked": false,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"### BEGIN SOLUTION\n",
"P, L, U = sp.linalg.lu(A)\n",
"\n",
"fig, ax = plt.subplots(1, 3, figsize=(10, 3))\n",
"ax[0].spy(P)\n",
"ax[0].set_title(\"$P$\")\n",
"ax[1].spy(L)\n",
"ax[1].set_title(\"$L$\")\n",
"ax[2].spy(U)\n",
"ax[2].set_title(\"$U$\")\n",
"plt.show()\n",
"\n",
"# We notice that the matrix P is not the identity matrix, this means that we can be sure\n",
"# that at least one permutation was performed. We see that A has almost only non-zero elements\n",
"# whereas in L and U, only around half of the elements are non-zero.\n",
"# The non-zero structure of L and U displays the triangular patterns typical of LU factorization. \n",
"\n",
"### END SOLUTION"
]
},
{
"cell_type": "markdown",
"id": "acb5d7e4-78fc-4462-abf2-a3fdc3ee6d50",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-c63d4ea99b0a801a",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"\n",
" \n",
"**Exercise 6:** For the binomial matrix $A$ with dimensions $m=20^2, 21^2,\\ldots,35^2$, compute the number of non-zero entries $\\text{nnz}(A)$ in $A$ and compare it to the number of non-zero elements $\\text{nnz}(L) + \\text{nnz}(U)$ in its LU factors. Plot the results in a logarithmic plot, and add the number of entries $m^2$ of the full matrix for comparison.\n",
"\n",
"*Hint:* Use `np.count_nonzero` to count the number of non-zero entries in a NumPy array.\n",
"
"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "24724b39-0ac7-4fc8-b30a-455567d87933",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-b98a075242a67026",
"locked": false,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"### BEGIN SOLUTION\n",
"m_list = np.arange(20, 35) ** 2\n",
"\n",
"nnz_list_A = []\n",
"nnz_list_A_LU = []\n",
"for m in m_list:\n",
" A = binomial_matrix(m, n, p)\n",
" P, L, U = sp.linalg.lu(A)\n",
" nnz_list_A.append(np.count_nonzero(A))\n",
" nnz_list_A_LU.append(np.count_nonzero(L) + np.count_nonzero(U))\n",
"\n",
"plt.loglog(m_list, nnz_list_A, label=\"non-zero elements in $A$\")\n",
"plt.loglog(m_list, nnz_list_A_LU, label=\"non-zero elements in $L$ and $U$\")\n",
"plt.loglog(m_list, m_list**2, ls=\"--\", c=\"k\", label=\"$m^2$\")\n",
"plt.ylabel(\"non-zero elements\")\n",
"plt.xlabel(\"matrix dimension $m$\")\n",
"plt.grid(True, which=\"major\", linestyle=\"-\")\n",
"plt.grid(True, which=\"minor\", linestyle=\"--\")\n",
"plt.legend()\n",
"plt.show()\n",
"\n",
"# Since A is almost a dense matrix, the number of non-zero elements are very close to m².\n",
"# Similarly, the number of non-zero elements of both the LU factors together are also approximately m².\n",
"### END SOLUTION"
]
},
{
"cell_type": "markdown",
"id": "e59d4257-2b03-4377-8760-de1d6de9277f",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-f7e107d62b3543c1",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"Below, we provide you with a function which computes the runtime of the LU factorization algorithm `sp.linalg.lu` for any matrix `A` by executing the function on the matrix multiple times, and averaging the individual runtimes. The number of repetitions of the function execution can be controlled with the argument `repeats`."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "0133d01f-7bf7-499f-8d77-2da5cd7f0402",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-55565969b2825125",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"outputs": [],
"source": [
"def lu_runtime(A, repeats):\n",
" t = timeit.timeit(lambda: sp.linalg.lu(A), number=repeats) / repeats\n",
" return t"
]
},
{
"cell_type": "markdown",
"id": "b328a03b-a1cd-4696-bc28-2dccaecd62b5",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-79c88eb1fc0ecc84",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"\n",
" \n",
"**Exercise 7:** Measure the approximate runtime for computing the LU factorization of the binomial matrix $A$ of dimensions $m=20^2, 21^2,\\ldots,35^2$. Set `repeats=50` to average over $50$ executions of the algorithm to get a more stable estimate of the true run-time.\n",
"
"
]
},
{
"cell_type": "markdown",
"id": "64c64a98-b216-4f07-a554-8c76e029a6a8",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-4de1828a62d21c86",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"\n",
"\n",
"**Note:** It may happen that the runtime of the first $m$ is significantly higher. This is due to the caching which Python performes in the background. If it bothers you, you may simply ignore the first $m$.\n",
"
"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "7d8978a9-20a3-4a8e-bf7b-c98490e0ece4",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-af07e1ccd799913c",
"locked": false,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"### BEGIN SOLUTION\n",
"t_list_A = []\n",
"for m in m_list:\n",
" A = binomial_matrix(m, n, p)\n",
" t = lu_runtime(A, repeats=50)\n",
" t_list_A.append(t)\n",
"\n",
"plt.loglog(m_list, t_list_A, label=\"LU factorization of $A$\")\n",
"plt.plot(m_list, 1e-8 * m_list ** 2, ls=\"--\", c=\"k\", label=\"$m^2$\")\n",
"plt.ylabel(\"runtime $t$\")\n",
"plt.xlabel(\"matrix dimension $m$\")\n",
"plt.grid(True, which=\"major\", linestyle=\"-\")\n",
"plt.grid(True, which=\"minor\", linestyle=\"--\")\n",
"plt.legend()\n",
"plt.show()\n",
"\n",
"# We notice that the run-time grows slightly stronger than m²\n",
"### END SOLUTION"
]
},
{
"cell_type": "markdown",
"id": "b0c6bc08-e93e-483f-8522-f7ba6f3c2249",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-c06bf4068c06f6f7",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"The following function `laplacian_matrix` generates a Laplacian matrix $B \\in \\mathbb{R}^{m \\times m}$."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "1e9a9525-fdd6-4c7f-9752-14b944987bc4",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-768b7fce1ea61dab",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"outputs": [],
"source": [
"def laplacian_matrix(m):\n",
" n = int(np.sqrt(m))\n",
" d = np.ones(n ** 2)\n",
" mat = sp.sparse.spdiags([d, -2 * d, d], [-1, 0, 1], n, n)\n",
" I = sp.sparse.eye(n)\n",
" A = sp.sparse.kron(I, mat) + sp.sparse.kron(mat, I) \n",
" return A.todense()"
]
},
{
"cell_type": "markdown",
"id": "e0612d6d-6a5a-4200-b29c-6dc2630fc6ed",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-0a38a19bf51d0502",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"\n",
" \n",
"**Exercise 8:** Repeat all of the above exercises in this section for the Laplacian matrix $B$ with the same dimensions. Compare to the results you've obtained for the binomial matrix $A$.\n",
"
"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "d21bb3a2-1600-43d5-8097-b81f63691dbd",
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-84676c464c1bde58",
"locked": false,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"### BEGIN SOLUTION\n",
"B = laplacian_matrix(400)\n",
"plt.spy(B)\n",
"plt.show()\n",
"\n",
"P, L, U = sp.linalg.lu(B)\n",
"\n",
"fig, ax = plt.subplots(1, 3, figsize=(10, 3))\n",
"ax[0].spy(P)\n",
"ax[0].set_title(\"$P$\")\n",
"ax[1].spy(L)\n",
"ax[1].set_title(\"$L$\")\n",
"ax[2].spy(U)\n",
"ax[2].set_title(\"$U$\")\n",
"plt.show()\n",
"\n",
"nnz_list_B = []\n",
"nnz_list_B_LU = []\n",
"for m in m_list:\n",
" B = laplacian_matrix(m)\n",
" P, L, U = sp.linalg.lu(B)\n",
" nnz_list_B.append(np.count_nonzero(B))\n",
" nnz_list_B_LU.append(np.count_nonzero(L) + np.count_nonzero(U))\n",
"\n",
"plt.loglog(m_list, nnz_list_B, label=\"non-zero elements in $B$\")\n",
"plt.loglog(m_list, nnz_list_B_LU, label=\"non-zero elements in $L$ and $U$\")\n",
"plt.loglog(m_list, m_list ** (3/2), ls=\"--\", c=\"k\", label=\"$m^{3/2}$\")\n",
"plt.ylabel(\"non-zero elements\")\n",
"plt.xlabel(\"matrix dimension $m$\")\n",
"plt.grid(True, which=\"major\", linestyle=\"-\")\n",
"plt.grid(True, which=\"minor\", linestyle=\"--\")\n",
"plt.legend()\n",
"plt.show()\n",
"\n",
"t_list_B = []\n",
"for m in m_list:\n",
" B = laplacian_matrix(m)\n",
" t = lu_runtime(B, repeats=50)\n",
" t_list_B.append(t)\n",
"\n",
"plt.loglog(m_list, t_list_B, label=\"LU factorization of $B$\")\n",
"plt.plot(m_list, 1e-8 * m_list ** 2, ls=\"--\", c=\"k\", label=\"$m^2$\")\n",
"plt.ylabel(\"runtime $t$\")\n",
"plt.xlabel(\"matrix dimension $m$\")\n",
"plt.grid(True, which=\"major\", linestyle=\"-\")\n",
"plt.grid(True, which=\"minor\", linestyle=\"--\")\n",
"plt.legend()\n",
"plt.show()\n",
"\n",
"# The LU factorization of a sparse matrix A often results in L and U matrices that are denser\n",
"# than A. This illustrates how factorization affects matrix sparsity, which is especially relevant\n",
"# for applications where memory and computational efficiency are key considerations.\n",
"# We notice that the computing time, which increases as m³ for the binomial case, in this case\n",
"# only increases proportionally to m². The memory consumption, which increases as m² for the\n",
"# binomial case, is here proportional to m^(3/2). This is due to the fact that the LU factorization\n",
"# of A will fill the band of the matrix, which has a width of √m. The matrices L and U therefore\n",
"# contain O(m√m) = O(m^(3/2)) non-zero elements, hence the factorization needs O(m√m²) = O(m²) operations.\n",
"### END SOLUTION"
]
},
{
"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 second 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
}