{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Numerical Methods\n",
"\n",
"## Lecture 6: Numerical Linear Algebra II\n",
"\n",
"### Exercise solutions"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Exercise 6.1: lower-triangular matrices\n",
"\n",
"Convince yourselves of the following facts:\n",
"\n",
"\n",
"* The multiplication of arbitrary lower-triangular square matrices is also lower-triangular.\n",
"\n",
"\n",
"* $L_2(L_1(L_0A)) = U \\implies A = L_0^{-1}(L_1^{-1}(L_2^{-1}U))$\n",
"\n",
"\n",
"* and hence that $A=LU$ where $U$ is the upper-triangular matrix found at the end of Guassian elimination, and where $L$ is the \n",
"following matrix\n",
"$$ L = L_0^{-1}L_1^{-1}L_2^{-1} $$\n",
"\n",
"\n",
"* Finally, compute this product of these lower-triangular matrices to show that \n",
"$$L = \n",
" \\begin{bmatrix}\n",
" {\\color{black}1} & {\\color{black}0} & {\\color{black}0} & {\\color{black}0}\\\\\n",
" {\\color{black}{1}} & {\\color{black}1} & {\\color{black}0} & {\\color{black}0}\\\\\n",
" {\\color{black}{4}} & {\\color{black}7} & {\\color{black}1} & {\\color{black}0}\\\\\n",
" {\\color{black}{5}} & {\\color{black}8} & {\\color{black}2} & {\\color{black}1}\\\\ \n",
" \\end{bmatrix}\n",
"$$\n",
"i.e. that the multiplication of these individual atomic matrices (importantly in this order) simply merges the entries from the non-zero columns of each atomic matrix, and hence is both lower-triangular, as well as trivial to compute."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"array([[ 5., 7., 5., 9.],\n",
" [ 5., 14., 7., 10.],\n",
" [20., 77., 41., 48.],\n",
" [25., 91., 55., 67.]])\n",
"array([[ 1, 0, 0, 0],\n",
" [-1, 1, 0, 0],\n",
" [-4, -7, 1, 0],\n",
" [-5, -8, 0, 1]])\n",
"array([[ 1, 0, 0, 0],\n",
" [ 0, 1, 0, 0],\n",
" [ 0, -7, 1, 0],\n",
" [ 0, 6, -2, 1]])\n",
"array([[5., 7., 5., 9.],\n",
" [0., 7., 2., 1.],\n",
" [0., 0., 7., 5.],\n",
" [0., 0., 0., 4.]])\n",
"array([[ 5., 7., 5., 9.],\n",
" [ 5., 14., 7., 10.],\n",
" [20., 77., 41., 48.],\n",
" [25., 91., 55., 67.]])\n",
"array([[1., 0., 0., 0.],\n",
" [1., 1., 0., 0.],\n",
" [4., 7., 1., 0.],\n",
" [5., 8., 2., 1.]])\n"
]
},
{
"data": {
"text/plain": [
"True"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import numpy as np\n",
"from scipy import linalg\n",
"from pprint import pprint\n",
"\n",
"# as above: A matrix,\n",
"A=np.array([[5., 7. , 5., 9.],\n",
" [5., 14., 7., 10.],\n",
" [20., 77., 41., 48.],\n",
" [25., 91., 55., 67.]])\n",
"# lower triangular matrices,\n",
"L0 = np.array([[1,0,0,0],[-1,1,0,0],[-4,0,1,0],[-5,0,0,1]])\n",
"L1 = np.array([[1,0,0,0],[0,1,0,0],[0,-7,1,0],[0,-8,0,1]])\n",
"L2 = np.array([[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,-2,1]])\n",
"# and their inverse matrices\n",
"L0_ = linalg.inv(L0)\n",
"L1_ = linalg.inv(L1)\n",
"L2_ = linalg.inv(L2)\n",
"\n",
"# our base matrix\n",
"pprint(A)\n",
"\n",
"# ad: The multiplication of arbitrary lower-triangular square matrices is also lower-triangular\n",
"pprint(np.dot(L0, L1))\n",
"pprint(np.dot(L2, L1))\n",
"\n",
"# ad: L2(L1(L0 A)) = U\n",
"U = np.dot(L2, np.dot(L1, np.dot(L0, A)))\n",
"pprint(U)\n",
"\n",
"# ad: A = L0^-1(L1^-1(L2^-1 U))\n",
"A = np.dot(L0_, np.dot(L1_, np.dot(L2_, U))) \n",
"pprint(A)\n",
"\n",
"# ad: L = L0^-1 L1^-1 L2^-1\n",
"L = np.dot(L0_,np.dot(L1_,L2_))\n",
"pprint(L)\n",
"\n",
"# ad: A = LU\n",
"# numpy.allclose takes two numpy arrays and checks whether all elements are identical to some tolerance\n",
"np.allclose(A, np.dot(L, U))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Exercise 6.2: LU decomposition\n",
"\n",
"Starting from your Gaussian elimination code produce a new code to compute the LU decomposition of a matrix. First, store $L$ and $U$ in two different matrices; you could then try a version where you store the entries of both $L$ and $U$ in $A$ as described above."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"array([[ 5., 7., 5., 9.],\n",
" [ 5., 14., 7., 10.],\n",
" [20., 77., 41., 48.],\n",
" [25., 91., 55., 67.]])\n",
"array([[1., 0., 0., 0.],\n",
" [1., 1., 0., 0.],\n",
" [4., 7., 1., 0.],\n",
" [5., 8., 2., 1.]])\n",
"array([[5., 7., 5., 9.],\n",
" [0., 7., 2., 1.],\n",
" [0., 0., 7., 5.],\n",
" [0., 0., 0., 4.]])\n"
]
},
{
"data": {
"text/plain": [
"True"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import numpy as np\n",
"import scipy.linalg as sl\n",
"from pprint import pprint\n",
"\n",
"def LU_dec(A):\n",
" # upper triangular matrix contains gaussian elimination result\n",
" # we won't change A in-place but create a local copy\n",
" A = A.copy()\n",
" m, n = A.shape\n",
" # lower triangular matrix has identity diagonal and scaling factors\n",
" # start from the identity:\n",
" L = np.identity(n)\n",
" # Loop over each pivot row.\n",
" for k in range(n):\n",
" # Loop over each equation below the pivot row.\n",
" for i in range(k+1, n):\n",
" # Define the scaling factor outside the innermost\n",
" # loop otherwise its value gets changed.\n",
" s = (A[i,k]/A[k,k])\n",
" for j in range(k, n):\n",
" A[i,j] = A[i,j] - s*A[k,j]\n",
" # scaling factors make up the lower matrix \n",
" L[i,k] = s\n",
" # A now is the upper triangular matrix U\n",
" return L, A\n",
"\n",
"\n",
"A = np.array([[ 5., 7., 5., 9.],\n",
" [ 5., 14., 7., 10.],\n",
" [20., 77., 41., 48.],\n",
" [25., 91. ,55., 67.]])\n",
"\n",
"L, U = LU_dec(A)\n",
"\n",
"pprint(A)\n",
"pprint(L)\n",
"pprint(U)\n",
"np.allclose(np.dot(L,U), A)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Exercise 6.3: Partial pivoting\n",
"\n",
"Implement partial pivoting."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"array([[ 5., 7., 5., 9.],\n",
" [ 5., 14., 7., 10.],\n",
" [20., 77., 41., 48.],\n",
" [25., 91., 55., 67.]])\n",
"array([[ 1. , 0. , 0. , 0. ],\n",
" [ 0.2 , 1. , 0. , 0. ],\n",
" [ 0.8 , -0.375 , 1. , 0. ],\n",
" [ 0.2 , 0.375 , 0.33333333, 1. ]])\n",
"array([[ 25. , 91. , 55. , 67. ],\n",
" [ 0. , -11.2 , -6. , -4.4 ],\n",
" [ 0. , 0. , -5.25 , -7.25 ],\n",
" [ 0. , 0. , 0. , 0.66666667]])\n",
"array([[ 5., 7., 5., 9.],\n",
" [ 5., 14., 7., 10.],\n",
" [20., 77., 41., 48.],\n",
" [25., 91., 55., 67.]])\n"
]
},
{
"data": {
"text/plain": [
"True"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import numpy as np\n",
"from pprint import pprint\n",
"\n",
"# Part of Ex. 3.3. This function is not necessary but makes the LU_dec_pp function less cluttered\n",
"# Making general operations into functions is good way of improving readability and reducing your workload later!\n",
"def swap_rows(A,j,k):\n",
" B = np.copy(A[j,:])\n",
" C = np.copy(A[k,:])\n",
" # second label, BUT NOT CREATING A 2ND ARRAY. Try it for yourself if you want.\n",
" A[k,:] = B\n",
" A[j,:] = C\n",
" return A\n",
"\n",
"# Ex. 3.3 A function to perform LU decomposition with partial pivoting\n",
"def LU_dec_pp(A):\n",
" m, n = A.shape\n",
" A = A.copy() # we won't modify in place but create local copy\n",
" P_ = np.identity(m)\n",
" L = np.identity(m)\n",
" for k in range(m):\n",
" j = np.argmax(abs(A[k:,k])) \n",
" # Find the index of the largest ABSOLUTE value. np.argmax will return \n",
" # the index of the largest element in an array\n",
" j+= k # A[1+2,2] = A[3,2] = 3! \n",
" A = swap_rows(A,j,k)\n",
" P_ = swap_rows(P_,j,k)\n",
" for j in range(k+1,m):\n",
" s = A[j,k]/A[k,k]\n",
" A[j,k:] -= A[k,k:]*s \n",
" L[j,k] = s\n",
" U = A \n",
" return P_.T, L, U\n",
"\n",
"A = np.array([[ 5., 7., 5., 9.],\n",
" [ 5., 14., 7., 10.],\n",
" [20., 77., 41., 48.],\n",
" [25., 91. ,55., 67.]])\n",
"\n",
"P, L, U = LU_dec_pp(A)\n",
"pprint(A)\n",
"pprint(L)\n",
"pprint(U)\n",
"pprint(np.dot(P, np.dot(L,U)))\n",
"np.allclose(np.dot(P, np.dot(L,U)), A)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Note the solution was given at the end of the solutions to last week's lecture."
]
}
],
"metadata": {
"celltoolbar": "Slideshow",
"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.7.6"
},
"toc": {
"base_numbering": 1,
"nav_menu": {},
"number_sections": false,
"sideBar": false,
"skip_h1_title": false,
"title_cell": "Table of Contents",
"title_sidebar": "Contents",
"toc_cell": false,
"toc_position": {
"height": "calc(100% - 180px)",
"left": "10px",
"top": "150px",
"width": "213.809px"
},
"toc_section_display": false,
"toc_window_display": false
}
},
"nbformat": 4,
"nbformat_minor": 1
}