{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Numerical Methods\n",
"\n",
"\n",
"# Lecture 5: Numerical Linear Algebra I\n",
"\n",
"## Exercise solutions"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import scipy.linalg as sl"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"### Exercise 5.1: Solving a linear system \n",
"Formulate and solve the linear system and check you get the answer quoted above."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We need to form $A$ and $\\pmb{b}$ and perform the following\n",
"\n",
"\\begin{align}\n",
"A\\pmb{x} & = \\pmb{b}\\\\\n",
"\\implies A^{-1}A\\pmb{x} & = A^{-1}\\pmb{b}\\\\\n",
"\\implies I\\pmb{x} & = A^{-1}\\pmb{b}\\\\\n",
"\\implies \\pmb{x} & = A^{-1}\\pmb{b}\n",
"\\end{align}\n",
"\n",
"so we can find the solution $\\pmb{x}$ by multiplying the inverse of $A$ with the RHS vector $\\pmb{b}$."
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"-11.0\n"
]
}
],
"source": [
"A=np.array([[2., 3.],[1., -4.]])\n",
"\n",
"# check first whether the determinant of A is non-zero - see below for reasons why.\n",
"print(sl.det(A)) "
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[3.36363636 0.09090909]\n"
]
}
],
"source": [
"b=np.array([7., 3.])\n",
"\n",
"# compute A inverse and multiply by b\n",
"print(sl.inv(A) @ b)"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[3.36363636 0.09090909]\n",
"Check it against the solution we calculated by hand earlier: 3.3636363636363638 0.09090909090909091\n"
]
}
],
"source": [
"# or solve the linear system using scipy - actually does \n",
"# the same thing as line above using LU decomposition (see later in lectures for the details)\n",
"print(sl.solve(A,b))\n",
"\n",
"print(\"Check it against the solution we calculated by hand earlier: \", 37./11., 1./11.)"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"True\n"
]
}
],
"source": [
"# this is a more rigorous way of checking the result - using numpy.allclose:\n",
"print(np.allclose(np.array([37./11., 1./11.]), sl.solve(A,b)))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Exercise 5.2: Matrix manipulation in python \n",
"\n",
"Let\n",
"$$\n",
"A = \\left(\n",
" \\begin{array}{ccc}\n",
" 1 & 2 & 3 \\\\\n",
" 4 & 5 & 6 \\\\\n",
" 7 & 8 & 9 \\\\\n",
" \\end{array}\n",
"\\right)\n",
"\\mathrm{\\quad\\quad and \\quad\\quad}\n",
"b = \\left(\n",
" \\begin{array}{c}\n",
" 2 \\\\\n",
" 4 \\\\\n",
" 6 \\\\\n",
" \\end{array}\n",
"\\right)\n",
"$$\n",
"\n",
"- Store $A$ and $b$ in NumPy array structures. Print them.\n",
"- Print their shape and size. What is the difference ?\n",
"- Create a NumPy array $I$ containing the identity matrix $I_3$. Do $A = A+I$.\n",
"- Substitute the 3rd column of $A$ with $b$.\n",
"- Solve $Ax=b$."
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"A = [[1 2 3]\n",
" [4 5 6]\n",
" [7 8 9]]\n",
"b = [2 4 6]\n",
"Size of A: 9 and shape of A: (3, 3)\n",
"Size of b: 3 and shape of b: (3,)\n",
"I = [[1. 0. 0.]\n",
" [0. 1. 0.]\n",
" [0. 0. 1.]]\n",
"A = [[ 2. 2. 3.]\n",
" [ 4. 6. 6.]\n",
" [ 7. 8. 10.]]\n",
"A = [[2. 2. 2.]\n",
" [4. 6. 4.]\n",
" [7. 8. 6.]]\n",
"x = [0. 0. 1.]\n"
]
}
],
"source": [
"A = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])\n",
"b = np.array([2, 4, 6])\n",
"print(\"A =\", A)\n",
"print(\"b = \",b)\n",
"\n",
"print(\"Size of A: \", A.size,\" and shape of A: \",A.shape)\n",
"print(\"Size of b: \", b.size,\" and shape of b: \",b.shape)\n",
"\n",
"I = np.eye(3)\n",
"print(\"I = \",I)\n",
"A = A + I\n",
"print(\"A = \",A)\n",
"\n",
"A[:, 2] = b\n",
"print(\"A = \",A)\n",
"\n",
"x = sl.solve(A,b)\n",
"print(\"x = \", x)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Exercise 5.3: Gaussian elimination $3 \\times 3$ example (by hand) \n",
"\n",
"Consider the system of linear equations\n",
"\n",
"\\begin{align*}\n",
" 2x + 3y - 4z &= 5 \\\\\n",
" 6x + 8y + 2z &= 3 \\\\\n",
" 4x + 8y - 6z &= 19\n",
"\\end{align*}\n",
"\n",
"write this in matrix form, form the corresponding augmented system and perform row operations until you get to upper-triangular form, find the solution using back substitution (**do this all with pen and paper**).\n",
"\n",
"Write some code to check your answer using `sl.inv(A) @ b`.\n",
"\n",
"You should find $x=-6$, $y=5$, $z=-1/2$."
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[-6. 5. -0.5]\n"
]
}
],
"source": [
"A = np.array([[2.,3.,-4.],[6.,8.,2.],[4.,8.,-6.]])\n",
"b = np.array([5.,3.,19.])\n",
"print(sl.inv(A) @ b)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Exercise 5.4: Gaussian elimination\n",
"\n",
"Write some code that takes a matrix $A$ and a vector $\\pmb{b}$ and converts it into upper-triangular form using the above algorithm. For the $2 \\times 2$ and $3\\times 3$ examples from above compare the resulting $A$ and $\\pmb{b}$ you obtain following elimination.\n"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Our A matrix following row operations to transform it into upper-triangular form:\n",
"array([[ 2. , 3. ],\n",
" [ 0. , -5.5]])\n",
"The correspondingly updated b vector:\n",
"array([ 7. , -0.5])\n",
"\n",
"Our A matrix following row operations to transform it into upper-triangular form:\n",
"array([[ 2., 3., -4.],\n",
" [ 0., -1., 14.],\n",
" [ 0., 0., 30.]])\n",
"The correspondingly updated b vector:\n",
"array([ 5., -12., -15.])\n"
]
}
],
"source": [
"def upper_triangle(A, b):\n",
" \"\"\" A function to covert A into upper triangluar form through row operations.\n",
" The same row operations are performed on the vector b.\n",
" \n",
" Note that this implementation does not use partial pivoting which is introduced below.\n",
" \n",
" Also note that A and b are overwritten, and hence we do not need to return anything\n",
" from the function.\n",
" \"\"\"\n",
" n = np.size(b)\n",
" rows, cols = np.shape(A)\n",
" # check A is square\n",
" assert(rows == cols)\n",
" # and check A has the same numner of rows as the size of the vector b\n",
" assert(rows == n)\n",
"\n",
" # Loop over each pivot row - all but the last row which we will never need to use as a pivot\n",
" for k in range(n-1):\n",
" # Loop over each row below the pivot row, including the last row which we do need to update\n",
" for i in range(k+1, n):\n",
" # Define the scaling factor for this row outside the innermost loop otherwise \n",
" # its value gets changed as you over-write A!!\n",
" # There's also a performance saving from not recomputing things when not strictly necessary\n",
" s = (A[i, k] / A[k, k])\n",
" # Update the current row of A by looping over the column j\n",
" # start the loop from k as we can assume the entries before this are already zero\n",
" for j in range(k, n):\n",
" A[i, j] = A[i, j] - s*A[k, j]\n",
" # and update the corresponding entry of b\n",
" b[i] = b[i] - s*b[k]\n",
"\n",
"\n",
"# Test our code on our 2x2 and 3x3 examples from above\n",
"\n",
"# 2x2 example from above\n",
"A = np.array([[2., 3.], [1., -4.]])\n",
"b = np.array([7., 3.])\n",
"\n",
"upper_triangle(A, b)\n",
"\n",
"\n",
"# Here is a new trick for you - \"pretty print\"\n",
"from pprint import pprint\n",
"\n",
"print('Our A matrix following row operations to transform it into upper-triangular form:')\n",
"pprint(A)\n",
"print('The correspondingly updated b vector:')\n",
"pprint(b)\n",
"\n",
"# 3x3 example from homework exercise\n",
"A = np.array([[2., 3., -4.],\n",
" [6., 8., 2.],\n",
" [4., 8., -6.]])\n",
"b = np.array([5., 3., 19.])\n",
"\n",
"upper_triangle(A, b)\n",
"\n",
"print('\\nOur A matrix following row operations to transform it into upper-triangular form:')\n",
"pprint(A)\n",
"print('The correspondingly updated b vector:')\n",
"pprint(b)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Exercise 5.5: Back substitution\n",
"\n",
"Extend your code to perform back substitution and hence to obtain the final solution $\\pmb{x}$. Check against the solutions found earlier. Come up with some random $n\\times n$ matrices (you can use `np.random.rand` for that) and check your code against `sl.inv(A)@b` (remember to use the original $A$ and $\\pmb{b}$ here of course!)\n"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Our solution: [-6. 5. -0.5]\n",
"SciPy solution: [-6. 5. -0.5]\n",
"Success: True\n"
]
}
],
"source": [
"# This function assumes that A is already an upper triangular matrix, \n",
"# e.g. we have already run our upper_triangular function if needed.\n",
"\n",
"def back_substitution(A, b):\n",
" \"\"\" Function to perform back subsitution on the system Ax=b.\n",
" \n",
" Returns the solution x.\n",
" \n",
" Assumes that A is on upper triangular form.\n",
" \"\"\"\n",
" n = np.size(b)\n",
" # Check A is square and its number of rows and columns same as size of the vector b\n",
" rows, cols = np.shape(A)\n",
" assert(rows == cols)\n",
" assert(rows == n)\n",
" # We can/should check that A is upper triangular using np.triu which is the \n",
" # upper triangular part of a matrix - if A is already upper triangular, then\n",
" # it should of course match the upper-triangular component of A!!\n",
" assert(np.allclose(A, np.triu(A)))\n",
" \n",
" x = np.zeros(n)\n",
" # start at the end (row n-1) and work backwards\n",
" for k in range(n-1, -1, -1):\n",
" # note that we could do this update in a single vectorised line \n",
" # using np.dot or @ - this could also speed things up\n",
" s = 0.\n",
" for j in range(k+1, n):\n",
" s = s + A[k, j]*x[j]\n",
" x[k] = (b[k] - s)/A[k, k]\n",
"\n",
" return x\n",
"\n",
"\n",
"# This A is the upper triangular matrix carried forward\n",
"# from the Python box above, and b the correspondingly updated b vector.\n",
"A = np.array([[ 2., 3., -4.],\n",
" [ 0., -1., 14.],\n",
" [ 0., 0., 30.]])\n",
"b = np.array([ 5., -12., -15.])\n",
"\n",
"# print the solution using our codes\n",
"x = back_substitution(A, b)\n",
"print('Our solution: ',x) \n",
"\n",
"# Reinitialise A and b !\n",
"# remember our functions overwrote them\n",
"A = np.array([[2., 3., -4.],\n",
" [6., 8., 2.],\n",
" [4., 8., -6.]])\n",
"b = np.array([5., 3., 19.])\n",
"\n",
"# check our answer against what SciPy gives us by multiplying b by A inverse \n",
"print('SciPy solution: ',sl.inv(A) @ b)\n",
"\n",
"print('Success: ', np.allclose(x, sl.inv(A) @ b))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Exercise 5.6: Matrix inversion\n",
"\n",
"Try updating your code to construct the inverse matrix. Check your answer against the inverse matrix given by a built-in function.\n",
"\n",
"Hint: Once you have performed your Gaussian elimination to transform A into an upper triangular matrix, perform another elimination \"from bottom to top\" to transform A into a diagonal matrix."
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Upper triangular transformed A = \n",
"array([[ 2. , 3. , -4. ],\n",
" [ 0. , -5.5 , 8. ],\n",
" [ 0. , 0. , 4.18181818]])\n",
"\n",
"and following application of our lower triangular function = \n",
"array([[ 2. , 0. , 0. ],\n",
" [ 0. , -5.5 , 0. ],\n",
" [ 0. , 0. , 4.18181818]])\n",
"\n",
"Our final transformed A = \n",
"array([[ 1., 0., 0.],\n",
" [-0., 1., -0.],\n",
" [ 0., 0., 1.]])\n",
"\n",
"and the correspondingly transformed B = \n",
"array([[ 0.13043478, 0.30434783, -0.04347826],\n",
" [-0.04347826, -0.43478261, 0.34782609],\n",
" [-0.2173913 , -0.17391304, 0.23913043]])\n",
"\n",
"SciPy computes the inverse as:\n",
"array([[ 0.13043478, 0.30434783, -0.04347826],\n",
" [-0.04347826, -0.43478261, 0.34782609],\n",
" [-0.2173913 , -0.17391304, 0.23913043]])\n",
"\n",
"Success: True\n"
]
}
],
"source": [
"# This updated version of the upper_triangular function now\n",
"# assumes that a matrix, B, is in the old vector location (what was b)\n",
"# in the augmented system, and applies the same operations to\n",
"# B as to A - only a minor difference\n",
"\n",
"def upper_triangle2(A, B):\n",
" n, m = np.shape(A)\n",
" assert(n == m) # this is designed to work for a square matrix\n",
"\n",
" # Loop over each pivot row.\n",
" for k in range(n-1):\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 as you are\n",
" # over-writing A\n",
" s = (A[i, k]/A[k, k])\n",
" for j in range(n):\n",
" A[i, j] = A[i, j] - s*A[k, j]\n",
" # replace the old b update with the same update as A\n",
" B[i, j] = B[i, j] - s*B[k, j]\n",
"\n",
"\n",
"# and this is a version which transforms the matrix into lower\n",
"# triangular form - the point here is that if you give it a\n",
"# matrix that is already in upper triangular form, then the\n",
"# result will be a diagonal matrix\n",
"def lower_triangle2(A, B):\n",
" n, m = np.shape(A)\n",
" assert(n == m) # this is designed to work for a square matrix\n",
"\n",
" # now it's basically just the upper triangular algorithm \n",
" # applied backwards\n",
" for k in range(n-1, -1, -1):\n",
" for i in range(k-1, -1, -1):\n",
" s = (A[i, k]/A[k, k])\n",
" for j in range(n):\n",
" A[i, j] = A[i, j] - s*A[k, j]\n",
" B[i, j] = B[i, j] - s*B[k, j]\n",
"\n",
"\n",
"# Let's redefine A as our matrix above\n",
"A = np.array([[2., 3., -4.], [3., -1., 2.], [4., 2., 2.]])\n",
"\n",
"# and B is the identity of the corresponding size\n",
"B = np.eye(np.shape(A)[0])\n",
"\n",
"# transform A into upper triangular form \n",
"# (and perform the same operations on B)\n",
"upper_triangle2(A, B)\n",
"print('Upper triangular transformed A = ')\n",
"pprint(A)\n",
"\n",
"# now make this updated A lower triangular as well \n",
"# (the result should be diagonal)\n",
"lower_triangle2(A, B)\n",
"print('\\nand following application of our lower triangular function = ')\n",
"pprint(A)\n",
"\n",
"# The final step to achieve the identity is just to divide each row through by the value \n",
"# of the diagonal to end up with 1's on the main diagonal and 0 everywhere else.\n",
"for i in range(np.shape(A)[0]):\n",
" B[i, :] = B[i, :]/A[i, i]\n",
" A[i, :] = A[i, :]/A[i, i]\n",
"\n",
"# the final A should be the identity\n",
"print('\\nOur final transformed A = ')\n",
"pprint(A)\n",
"\n",
"# the final B should therefore be the inverse of the original B\n",
"print('\\nand the correspondingly transformed B = ')\n",
"pprint(B)\n",
"\n",
"# let's compute the inverse using built-in functions and check\n",
"# we get the same answer (we need to reinitialise A)\n",
"A = np.array([[2., 3., -4.], [3., -1., 2.], [4., 2., 2.]])\n",
"print('\\nSciPy computes the inverse as:')\n",
"pprint(sl.inv(A))\n",
"\n",
"# B should now store the inverse of the original A - let's check\n",
"print('\\nSuccess: ', np.allclose(B, sl.inv(A)))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Exercise 5.7: Zeros on the diagonal\n",
"\n",
"You may have noticed above that we have no way of guaranteeing that the $A_{kk}$ we divide through by in the Guassian elimination or back substitution algorithms is non-zero (or not very small which will also lead to computational problems).\n",
"Note also that we commented that we are free to exchange two rows in our augmented system - how could you use this fact to build robustness into our algorithms in order to deal with matrices for which our algorithms do lead to very small or zero $A_{kk}$ values? \n",
"\n",
"See if you can figure out how to do this - more on this next week!"
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"swapping rows 2 and 0\n",
"\n",
"A and b with row swaps: \n",
"array([[ 4. , 2. , 2. ],\n",
" [ 0. , -2.5, 0.5],\n",
" [ 0. , 0. , -4.6]])\n",
"array([ 8. , -3. , 3.6])\n",
"\n",
"A and b without any row swaps: \n",
"array([[ 2. , 3. , -4. ],\n",
" [ 0. , -5.5 , 8. ],\n",
" [ 0. , 0. , 4.18181818]])\n",
"array([ 10. , -12. , -3.27272727])\n",
"\n",
"These two upper triangular systems are equivalent (i.e. have the same solution): True\n"
]
}
],
"source": [
"# This function swaps rows in matrix A\n",
"# (and remember that we need to do likewise for the vector b \n",
"# we are performing the same operations on)\n",
"\n",
"def swap_row(A, b, i, j):\n",
" \"\"\" Swap rows i and j of the matrix A and the vector b.\n",
" \"\"\" \n",
" if i == j:\n",
" return\n",
" print('swapping rows', i,'and', j)\n",
" # If we are swapping two values, we need to take a copy of one of them first otherwise\n",
" # we will lose it when we make the first swap and will not be able to use it for the second.\n",
" # We need to make sure it is a real copy - not just a copy of a reference to the data!\n",
" # use np.copy to do this. \n",
" iA = np.copy(A[i, :])\n",
" ib = np.copy(b[i])\n",
"\n",
" A[i, :] = A[j, :]\n",
" b[i] = b[j]\n",
"\n",
" A[j, :] = iA\n",
" b[j] = ib\n",
"\n",
" \n",
"# This is a new version of the upper_triangular function\n",
"# with the added step of swapping rows so the largest\n",
"# magnitude number is always our pivot/\n",
"# pp stands for partial pivoting which will be explained\n",
"# in more detail below.\n",
"\n",
"def upper_triangle_pp(A, b):\n",
" \"\"\" A function to covert A into upper triangluar form through row operations.\n",
" The same row operations are performed on the vector b.\n",
" \n",
" This version uses partial pivoting.\n",
" \n",
" Note that A and b are overwritten, and hence we do not need to return anything\n",
" from the function.\n",
" \"\"\"\n",
" n = np.size(b)\n",
" # check A is square and its number of rows and columns same as size of the vector b\n",
" rows, cols = np.shape(A)\n",
" assert(rows == cols)\n",
" assert(rows == n)\n",
"\n",
" # Loop over each pivot row - all but the last row\n",
" for k in range(n-1):\n",
" # Swap rows so we are always dividing through by the largest number.\n",
" # initiatise kmax with the current pivot row (k)\n",
" kmax = k\n",
" # loop over all entries below the pivot and select the k with the largest abs value\n",
" for i in range(k+1, n):\n",
" if abs(A[kmax, k]) < abs(A[i, k]):\n",
" kmax = i\n",
" # and swap the current pivot row (k) with the row with the largest abs value below the pivot\n",
" swap_row(A, b, kmax, k)\n",
"\n",
" for i in range(k+1, n):\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",
" b[i] = b[i] - s*b[k]\n",
"\n",
"\n",
"# Apply the new code with row swaps to our matrix problem from above\n",
"A = np.array([[2., 3., -4.],\n",
" [3., -1., 2.],\n",
" [4., 2., 2.]])\n",
"b = np.array([10., 3., 8.])\n",
"\n",
"upper_triangle_pp(A, b)\n",
"\n",
"print('\\nA and b with row swaps: ')\n",
"pprint(A)\n",
"pprint(b)\n",
"# compute the solution from these using our back substitution code\n",
"# could also have used SciPy of course\n",
"x1 = back_substitution(A, b)\n",
"\n",
"# compare with our first function with no row swaps\n",
"A = np.array([[2., 3., -4.],\n",
" [3., -1., 2.],\n",
" [4., 2., 2.]])\n",
"b = np.array([10., 3., 8.])\n",
"\n",
"upper_triangle(A, b)\n",
"\n",
"print('\\nA and b without any row swaps: ')\n",
"pprint(A)\n",
"pprint(b)\n",
"x2 = back_substitution(A, b)\n",
"\n",
"# check these two systems are equivalent\n",
"print('\\nThese two upper triangular systems are equivalent (i.e. have the same solution): ',np.allclose(x1, x2))"
]
}
],
"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": {},
"toc_section_display": false,
"toc_window_display": false
}
},
"nbformat": 4,
"nbformat_minor": 1
}