{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "ChEn-3170: Computational Methods in Chemical Engineering Fall 2018 UMass Lowell; Prof. V. F. de Almeida **24Sep2018**\n", "\n", "# 05. Computational Linear Algebra Fundamentals\n", "$ \n", " \\newcommand{\\Amtrx}{\\boldsymbol{\\mathsf{A}}}\n", " \\newcommand{\\Bmtrx}{\\boldsymbol{\\mathsf{B}}}\n", " \\newcommand{\\Cmtrx}{\\boldsymbol{\\mathsf{C}}}\n", " \\newcommand{\\Dmtrx}{\\boldsymbol{\\mathsf{D}}}\n", " \\newcommand{\\Mmtrx}{\\boldsymbol{\\mathsf{M}}}\n", " \\newcommand{\\Imtrx}{\\boldsymbol{\\mathsf{I}}}\n", " \\newcommand{\\Pmtrx}{\\boldsymbol{\\mathsf{P}}}\n", " \\newcommand{\\Qmtrx}{\\boldsymbol{\\mathsf{Q}}}\n", " \\newcommand{\\Lmtrx}{\\boldsymbol{\\mathsf{L}}}\n", " \\newcommand{\\Umtrx}{\\boldsymbol{\\mathsf{U}}}\n", " \\newcommand{\\xvec}{\\boldsymbol{\\mathsf{x}}}\n", " \\newcommand{\\avec}{\\boldsymbol{\\mathsf{a}}}\n", " \\newcommand{\\bvec}{\\boldsymbol{\\mathsf{b}}}\n", " \\newcommand{\\cvec}{\\boldsymbol{\\mathsf{c}}}\n", " \\newcommand{\\rvec}{\\boldsymbol{\\mathsf{r}}}\n", " \\newcommand{\\norm}[1]{\\bigl\\lVert{#1}\\bigr\\rVert}\n", " \\DeclareMathOperator{\\rank}{rank}\n", "$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "## Table of Contents\n", "* [Theory](#theory)\n", "* [Matrix-vector and matrix-matrix product operations](#product)\n", "* [NumPy and SciPy Linear Algebra](#pylinalg)\n", " + [Matrix solve](#pysolve)\n", " + [$\\Lmtrx$ forward solve](#pyl)\n", " + [$\\Umtrx$ backward solve](#pyu)\n", " + [$\\Amtrx = \\Pmtrx\\,\\Lmtrx\\,\\Umtrx$ factorization](#[pyplu)\n", "* [Course Linear Algebra](#courselinalg)\n", " + [$\\Lmtrx$ forward solve](#l)\n", " + [$\\Umtrx$ forward solve](#u)\n", " + [$\\Amtrx = \\Lmtrx\\,\\Umtrx$ factorization](#lu)\n", " + [$\\Pmtrx\\,\\Amtrx = \\Lmtrx\\,\\Umtrx$ factorization](#plu)\n", " + [$\\Pmtrx\\,\\Amtrx\\,\\Qmtrx = \\Lmtrx\\,\\Umtrx$ factorization](#pqlu)\n", "---" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Theory\n", "The course notes (OneNote [ChEn-3170-linalg](https://studentuml-my.sharepoint.com/:o:/g/personal/valmor_dealmeida_uml_edu/EhHdq8qIW_JNiMBu60tKHi8BV6Nv8jonllo__NOtX7JgkQ?e=jA5xj8)) cover basic elements of lineary system of algebraic equations. Their representation in row and column formats and the insight gained from them when searching for a solution of the system.\n", "\n", "Some basic theoretical aspects of solving for $\\xvec$ in the matrix equation $\\Amtrx\\,\\xvec = \\bvec$ are covered. $\\Amtrx$ is a matrix, $\\bvec$ and $\\xvec$ are vectors." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Matrix-vector and matrix-matrix product operations\n", "The following operations between vectors and matrices are obtained directly from the buil-in functions in the `numpy` package." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "'''Import the NumPy package as usual'''\n", "\n", "import numpy as np" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Scalar product of two vectors: $\\avec \\cdot \\bvec$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "'''Vector scalar product or dot product of vectors'''\n", "\n", "a_vec = np.array( np.random.random(3) )\n", "b_vec = np.array( np.random.random(3) )\n", "\n", "a_vec_dot_b_vec = np.dot( a_vec, b_vec ) # clear linear algebra operation\n", "print('a.b =', a_vec_dot_b_vec)\n", "\n", "a_vec_x_b_vec = a_vec @ b_vec # consistent linear algebra multiplication\n", "print('a@b =', a_vec_x_b_vec )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Matrix vector product: $\\Amtrx\\,\\bvec$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "'''Matrix-vector product'''\n", "\n", "a_mtrx = np.array( [ [ 2., 1., 1.], # per course notes (p 04)\n", " [ 4., -6., 0.],\n", " [-2., 7., 2.]])\n", "\n", "b_vec = np.array( [5., -2., 9.]) # per course notes\n", "\n", "a_mtrx_x_b_vec = a_mtrx @ b_vec # linear algebra matrix-vector product\n", "\n", "print('A x b =', a_mtrx_x_b_vec)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Matrix-vector product: $\\Imtrx\\,\\bvec = \\bvec$. Note: $\\begin{pmatrix}\n", "1 & 0 & 0 \\\\\n", "0 & 1 & 0 \\\\\n", "0 & 0 & 1\n", "\\end{pmatrix} \\, \\begin{pmatrix} b_1\\\\b_2\\\\b_3 \\end{pmatrix} = \\begin{pmatrix} b_1\\\\b_2\\\\b_3 \\end{pmatrix} $." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "'''Identity-matrix vector product'''\n", "\n", "i_mtrx = np.eye(3)\n", "\n", "i_mtrx_x_b_vec = i_mtrx @ b_vec # linear algebra matrix-vector product\n", "\n", "print('I x b =', i_mtrx_x_b_vec)\n", "print('b =', b_vec)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Matrix-matrix product: $\\Imtrx\\,\\Amtrx = \\Amtrx$. Note: $\\begin{pmatrix}\n", "1 & 0 & 0 \\\\\n", "0 & 1 & 0 \\\\\n", "0 & 0 & 1\n", "\\end{pmatrix} \\, \n", "\\begin{pmatrix} \n", "A_{1,1} & A_{1,2} & A_{1,3} \\\\\n", "A_{2,1} & A_{2,2} & A_{2,3} \\\\\n", "A_{3,1} & A_{3,2} & A_{3,3}\n", "\\end{pmatrix} = \n", "\\begin{pmatrix} \n", "A_{1,1} & A_{1,2} & A_{1,3} \\\\\n", "A_{2,1} & A_{2,2} & A_{2,3} \\\\\n", "A_{3,1} & A_{3,2} & A_{3,3}\n", "\\end{pmatrix}\n", "$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "'''Matrix-matrix product IA = A'''\n", "\n", "i_mtrx_x_a_mtrx = i_mtrx @ a_mtrx # linear algebra matrix-matrix product\n", "\n", "print('I x A =\\n', i_mtrx_x_a_mtrx)\n", "print('A =\\n', a_mtrx)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Matrix-matrix product: $\\Amtrx\\,\\Bmtrx = \\Cmtrx$. Note: \n", "$\\begin{pmatrix}\n", "A_{1,1} & A_{1,2} & A_{1,3} \\\\\n", "A_{2,1} & A_{2,2} & A_{2,3} \\\\\n", "A_{3,1} & A_{3,2} & A_{3,3}\n", "\\end{pmatrix} \\, \n", "\\begin{pmatrix} \n", "B_{1,1} & B_{1,2} & B_{1,3} \\\\\n", "B_{2,1} & B_{2,2} & B_{2,3} \\\\\n", "B_{3,1} & B_{3,2} & B_{3,3}\n", "\\end{pmatrix} = \n", "\\begin{pmatrix} \n", "C_{1,1} & C_{1,2} & C_{1,3} \\\\\n", "C_{2,1} & C_{2,2} & C_{2,3} \\\\\n", "C_{3,1} & C_{3,2} & C_{3,3}\n", "\\end{pmatrix}\n", "$ where each $C_{i,j}$ is a vector product of the $i$th row of $\\Amtrx$ and the $j$th column of $\\Bmtrx$, *i.e.* \n", "$C_{i,j} = \\sum\\limits_{k=1}^3 A_{i,k}\\, B_{k,j}$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "'''Matrix-matrix product AB = C'''\n", "\n", "b_mtrx = np.array( [[5. , 5. , 5.],\n", " [-2., -2., -2.],\n", " [9. , 9. , 9.]]\n", " )\n", "c_mtrx = a_mtrx @ b_mtrx # linear algebra matrix-matrix product\n", "\n", "print('A =\\n', a_mtrx)\n", "print('B =\\n', b_mtrx)\n", "print('C =\\n', c_mtrx)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "The matrix-matrix product: $\\Bmtrx\\,\\Amtrx = \\Dmtrx \\ne \\Cmtrx$, does not commute in general. Note \n", "$D_{i,j} = \\sum\\limits_{k=1}^3 B_{i,k}\\, A_{k,j}$.\n", "
" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "'''Matrix-matrix product BA = D'''\n", "\n", "b_mtrx = np.array( [[5. , 5. , 5.],\n", " [-2., -2., -2.],\n", " [9. , 9. , 9.]]\n", " )\n", "d_mtrx = b_mtrx @ a_mtrx # linear algebra matrix-matrix product\n", "\n", "print('A =\\n', a_mtrx)\n", "print('B =\\n', b_mtrx)\n", "print('D =\\n', d_mtrx)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## NumPy and SciPy Linear Algebra\n", "[NumPy](http://www.numpy.org/) has extensive support for [linear algebra](https://docs.scipy.org/doc/numpy/reference/routines.linalg.html?highlight=linear%20algebra) arrays. We collect here the relevant operations for this course.\n", "However additional resources are instead added to [SciPy](https://docs.scipy.org/doc/scipy-1.1.0/reference/) for general scientific computing including [linear algebra](https://docs.scipy.org/doc/scipy-1.1.0/reference/tutorial/linalg.html).\n", "\n", "Linear algebra operations are obtained from the `linalg` sub-package of the `numpy` package, and the `linalg` sub-package of `scipy`." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "'''Import the NumPy linear algebra sub-package as usual'''\n", "\n", "#import numpy.linalg as linalg\n", "#from numpy import linalg # often used alternative\n", "'''or leave it commented since the usage of np.linalg is self-documenting'''" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The 2-norm or norm (magnitude) of a vector $\\bvec$ is indicated as $\\norm{\\bvec}$ and computed as follows:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "'''Vector norm (or magnitude)'''\n", "\n", "norm_b_vec = np.linalg.norm( b_vec ) # default norm is the 2-norm \n", "print('||b|| =', norm_b_vec) # same as magnitude" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Matrix solve \n", "\n", "*Solve* for $\\xvec$ in the matrix equation $\\Amtrx\\,\\xvec = \\bvec$, where $\\Amtrx = \n", "\\begin{pmatrix}\n", "2 & 1 & 1 \\\\\n", "4 & -6 & 0 \\\\\n", "-2 & 7 & 2\n", "\\end{pmatrix}\n", "$\n", "and $\\bvec = \\begin{pmatrix} 5\\\\ -2\\\\ 9 \\end{pmatrix}$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "'''Matrix solver (this is short for solution of a linear algebraic system of equations)'''\n", "\n", "x_vec = np.linalg.solve( a_mtrx, b_vec ) # solve linear system for A, b; per course notes 04\n", "print('solution x =', x_vec)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The residual vector defined as $\\rvec = \\bvec - \\Amtrx\\,\\xvec$ is of importance. So is its norm $\\norm{\\rvec}$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "'''Verify the accuracy of the solution'''\n", "\n", "res_vec = b_vec - a_mtrx @ x_vec\n", "print('b - A x =',res_vec)\n", "print('||b - A x|| =',np.linalg.norm( res_vec ))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The rank of a matrix of coefficients, $\\rank(\\Amtrx)$, of a linear algebraic system of equations determines weather the solution is unique or singular." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "'''Matrix rank'''\n", "\n", "k = np.linalg.matrix_rank( a_mtrx ) # rank; per course notes 14\n", "print('rank(A) =',k)\n", "print('shape(A) =',a_mtrx.shape)\n", "\n", "if k == a_mtrx.shape[0] and k == a_mtrx.shape[1]: # flow control\n", " print('A is non-singular; solution is unique ')\n", "else: \n", " print('A is singular')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Why is this matrix $\\Bmtrx$ singular?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "b_mtrx = np.array( [ [ 2., 1., 3.], # singular\n", " [ 4., -6., -2.],\n", " [-2., 7., 5.]])\n", "\n", "k = np.linalg.matrix_rank( b_mtrx ) # rank \n", "print('rank(B) =',k)\n", "print('shape(B) =',b_mtrx.shape)\n", "\n", "if k == b_mtrx.shape[0] and k == b_mtrx.shape[1]: # flow control\n", " print('B is non-singular; solution is unique ')\n", "else: \n", " print('B is singular')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "'''Matrix determinant'''\n", "\n", "det_a_mtrx = np.linalg.det( a_mtrx ) # determinant; course notes 16\n", "print('det(A) =', det_a_mtrx)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The inverse matrix is denoted as $\\Amtrx^{-1}$ and is computed as the matrix that multiplies $\\bvec$ and produces the solution $\\xvec$, that is, $\\xvec = \\Amtrx^{-1}\\,\\bvec$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "'''Matrix inverse'''\n", "\n", "a_mtrx_inv = np.linalg.inv( a_mtrx ) # matrix inverse; per course notes 17\n", "print('A^-1 =\\n', a_mtrx_inv)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Recall $\\Amtrx^{-1}\\,\\Amtrx = \\Imtrx$ where $\\Imtrx$ is the identity matrix." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "'''Identity matrix'''\n", "\n", "i_mtrx = a_mtrx_inv @ a_mtrx # identity matrix; per course notes 17\n", "print('A^-1 A =\\n',i_mtrx)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Using the inverse, the same solution will be found: $\\xvec = \\Amtrx^{-1}\\,\\bvec$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "'''Solution using the inverse'''\n", "\n", "x_vec_again = a_mtrx_inv @ b_vec # matrix-vector multiply; per course notes 17\n", "print('solution x =', x_vec_again)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This is the element-by-element reciprocal of the matrix $(\\Amtrx)^{-1}$, which is very different than the inverse." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "'''Inverse power of a matrix'''\n", "\n", "#a_mtrx_to_negative_1 = a_mtrx**(-1) # this will cause an error (division by zero)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's look at the determinant of a larger matrix, say $\\Mmtrx$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "'''Generate a larger matrix from an image'''\n", "\n", "from matplotlib import pyplot as plt # import the pyplot function of the matplotlib package\n", "plt.rcParams['figure.figsize'] = [20, 4] # extend the figure size on screen output\n", "mtrx = plt.imread('https://raw.githubusercontent.com/dpploy/chen-3170/master/notebooks/images/cermet.png',format='png')\n", "m_mtrx = mtrx[:,:]\n", "plt.figure(1) # create a figure placeholder\n", "plt.imshow( m_mtrx,cmap='gray')\n", "plt.title('Cermet',fontsize=14)\n", "plt.xlabel('x pixels',fontsize=12)\n", "plt.ylabel('y pixels',fontsize=12)\n", "print('M shape =', m_mtrx.shape)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "'''Larger matrix determinant'''\n", "\n", "det_m_mtrx = np.linalg.det( m_mtrx ) # determinant\n", "print('max(M) =',m_mtrx.max())\n", "print('min(M) =',m_mtrx.min())\n", "\n", "print('det(M) = %10.3e (not an insightful number)'%det_m_mtrx) # formatting numeric output\n", "\n", "print('rank(M) = ',np.linalg.matrix_rank( m_mtrx, tol=1e-5 ) )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's *solve* for this matrix with $\\cvec$ as the right side vector, that is, $\\Mmtrx\\,\\xvec = \\cvec$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "'''Solve M x = c and plot x'''\n", "\n", "c_vec = np.random.random(mtrx.shape[0]) # any c will do it\n", "\n", "sol = np.linalg.solve( m_mtrx, c_vec ) # solve linear system for A, b\n", "\n", "plt.figure(2)\n", "plt.plot(range(c_vec.size),sol,'k')\n", "plt.title('M x = c',fontsize=20)\n", "plt.xlabel('n',fontsize=18)\n", "plt.ylabel('$c_j$',fontsize=18)\n", "print('')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "res_vec = c_vec - m_mtrx @ sol\n", "#print('c - M x =',res_vec)\n", "print('||c - M x|| =%12.4e'%np.linalg.norm( res_vec ))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### $\\Lmtrx$ forward solve\n", "A lower triangular matrix allows for a forward solve." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "'''L forward solve'''\n", "\n", "l_mtrx = np.array( [[1., 0., 0.], # per course notes \n", " [2., 3., 0.],\n", " [4., 5., 6.]] )\n", "\n", "b_vec = np.array( [1.,2.,3.] )\n", "\n", "x_vec = np.linalg.solve( l_mtrx, b_vec )\n", "\n", "np.set_printoptions(precision=3) # one way to control printing of numpy arrays\n", "print('x = ',x_vec)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### $\\Umtrx$ backward solve\n", "An upper triangular matrix allows for a backward solve." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "'''U backward solve'''\n", "\n", "u_mtrx = np.array( [[1., 2., 3.], # per course notes\n", " [0, 4., 5.],\n", " [0., 0., 6.]] )\n", "\n", "b_vec = np.array( [1.,2.,3.] )\n", "\n", "x_vec = np.linalg.solve( u_mtrx, b_vec )\n", "\n", "np.set_printoptions(precision=3) # one way to control printing of numpy arrays\n", "print('x = ',x_vec)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### $\\Pmtrx\\,\\Lmtrx\\,\\Umtrx$ factorization\n", "The factors: $\\Pmtrx$, $\\Lmtrx$, and $\\Umtrx$ where $\\Pmtrx\\,\\Lmtrx\\,\\Umtrx = \\Amtrx$ can be obtained from the SciPy linear algebra package. $\\Pmtrx$ is a permutation matrix if the underlying Gaussian elimination used to construct the $\\Lmtrx$ and $\\Umtrx$ factors." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "'''Import only the linear algebra package'''\n", "\n", "import scipy.linalg\n", "import numpy as np" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false }, "outputs": [], "source": [ "'''P L U factors of A'''\n", "\n", "a_mtrx = np.array( [[1, 2, 3],\n", " [4, 5, 6],\n", " [7, 8, 10]] )\n", "\n", "(p_mtrx, l_mtrx, u_mtrx) = scipy.linalg.lu( a_mtrx )\n", "\n", "print('P =\\n',p_mtrx)\n", "print('L =\\n',l_mtrx)\n", "print('U =\\n',u_mtrx)\n", "print('Checking...')\n", "print('PLU - A =\\n', p_mtrx @ l_mtrx @ u_mtrx - a_mtrx)\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false }, "outputs": [], "source": [ "'''P^-1 = P^T (i.e. the transpose of a permutation matrix is its inverse)'''\n", "\n", "pinv_mtrx = np.linalg.inv(p_mtrx)\n", "print('P^-1 =\\n', pinv_mtrx)\n", "print('Checking...')\n", "print('P^-1 - P^T =\\n', pinv_mtrx - p_mtrx.transpose())\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "'''PLU x = b; that is: Forward: L y = P^-1 b, Backward: U x = y '''\n", "\n", "b_vec = np.array([1.,2.,3.])\n", "y_vec = scipy.linalg.solve(l_mtrx, p_mtrx.transpose() @ b_vec) # L y = P^T b\n", "\n", "x_vec = scipy.linalg.solve(u_mtrx, y_vec) # U x = y\n", "\n", "print('x =', x_vec)\n", "\n", "x_vec_gold = scipy.linalg.solve( a_mtrx, b_vec ) # solution using A x = b\n", "print('||x - x_gold|| =',scipy.linalg.norm(x_vec-x_vec_gold))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "'''Deterninant of U or L: product of the diagonal'''\n", "\n", "det_u = np.linalg.det(u_mtrx)\n", "print('det(U) = %8.3e'%det_u)\n", "\n", "diag_vec = np.diagonal(u_mtrx)\n", "prod = np.prod(diag_vec)\n", "print('diag(U) product = %8.3e'%prod )" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "'''Determinant of P (always +1 or -1)'''\n", "\n", "det_p = np.linalg.det(p_mtrx)\n", "print('det(P) = %8.3e'%det_p)\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "'''Determinant of A = det(PLU)'''\n", "\n", "det_l = np.prod( np.diagonal(l_mtrx) )\n", "det_plu = det_p * det_l * det_u # last term is det of L\n", "\n", "print('det(PLU) = %8.3e'%det_plu)\n", "print('det(A) = %8.3e'%np.linalg.det(a_mtrx))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## ChEn-3170 Linear Algebra\n", "In this course various algorithms need to be programmed. These should be compared to `SciPy` and/or `NumPy`." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### $\\Lmtrx$ forward solve\n", "A lower triangular matrix allows for a forward solve.\n", "The algorithm for $\\Lmtrx\\,\\xvec=\\bvec$ is as follows: \n", "\n", "\\begin{equation*}\n", "x_i = \\Bigl(b_i - \\sum\\limits_{j=1}^{i-1} L_{i,j}\\,x_j \\Bigr)\\,L^{-1}_{i,i} \\quad\\ \\forall \\quad\\ i=1,\\ldots,m\n", "\\end{equation*}\n", "\n", "for $i$ and $j$ with offset 1. Recall that `NumPy` and `Python` have offset 0 for their sequence data types." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "'''L forward solve'''\n", "\n", "l_mtrx = np.array( [[1., 0., 0.], # per course notes \n", " [2., 3., 0.],\n", " [4., 5., 6.]] )\n", "\n", "b_vec = np.array( [1.,2.,3.] )\n", "\n", "from chen_3170.help import forward_solve \n", "\n", "x_vec = forward_solve( l_mtrx, b_vec )\n", "\n", "np.set_printoptions(precision=3) # one way to control printing of numpy arrays\n", "print('x = ',x_vec)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### $\\Umtrx$ backward solve\n", "A upper triangular matrix allows for a backward solve.\n", "The algorithm for $\\Umtrx\\,\\xvec=\\bvec$ is as follows: \n", "\n", "\\begin{equation*} x_i = \\Bigl(b_i - \\sum\\limits_{j=i+1}^{m} U_{i,j}\\,x_j \\Bigr)\\,U^{-1}_{i,i} \\quad\\ \\forall \\quad\\ i=m,\\ldots,1\n", "\\end{equation*}\n", "\n", "for $i$ and $j$ with offset 1. Recall that `NumPy` and `Python` have offset 0 for their sequence data types." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "'''U backward solve'''\n", "\n", "u_mtrx = np.array( [[1., 2., 3.], # per course notes\n", " [0, 4., 5.],\n", " [0., 0., 6.]] )\n", "\n", "b_vec = np.array( [1.,2.,3.] )\n", "\n", "try: \n", " from chen_3170.toolkit import backward_solve \n", "except ModuleNotFoundError:\n", " assert False, 'You need to provide your own backward_solve function here. Bailing out.'\n", "\n", "x_vec = backward_solve( u_mtrx, b_vec )\n", "\n", "np.set_printoptions(precision=3) # one way to control printing of numpy arrays\n", "print('x = ',x_vec)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### $\\Amtrx = \\Lmtrx\\,\\Umtrx$ factorization\n", "$\\Lmtrx\\,\\Umtrx$ factorization algorithm (without using pivoting) for a square matrix $\\overset{(m \\times m)}{\\Amtrx}$ computes the $\\Lmtrx\\,\\Umtrx$ factors. The factorization is obtained by elimination steps $k = 1,\\ldots,m-1$ so that \n", "\n", "\\begin{equation*}\n", " A^{(k+1)}_{i,j} = A^{(k)}_{i,j} - A^{(k)}_{k,j}\\, m_{i,k} \\quad\\ \\forall\\ i=k+1,\\ldots,m \\quad\\ \\text{and} \\quad\\ j=k,\\ldots,m\n", "\\end{equation*}\n", "\n", "where the multipliers $m_{i,k}$ are given by $m_{i,k} = \\frac{A^{(k)}_{i,k}}{A^{(k)}_{k,k}}$. When $k = m-1$, $A^{(m)}_{i,j}$, is upper triangular, that is, $U_{i,j} = A^{(m)}_{i,j}$ . The lower triangular matrix is obtained using the multipliers $m_{i,k}$, that is $L_{i,j} = m_{i,j} \\ \\forall \\ i>j$, $L_{i,i}=1$, and $L_{i,j}=0 \\ \\forall \\ i factorization\n", "The factors: $\\Pmtrx$, $\\Lmtrx$, and $\\Umtrx$ where $\\Lmtrx\\,\\Umtrx = \\Pmtrx\\,\\Amtrx$ can be obtained from **partial pivoting** strategy to the $\\Lmtrx\\,\\Umtrx$ factorization algorithm shown above. $\\Pmtrx$ is a row permutation matrix obtained by the underlying Gaussian elimination used to construct the $\\Lmtrx$ and $\\Umtrx$ factors.\n", "\n", "Program a $\\Lmtrx\\,\\Umtrx$ factorization algorithm (using partial pivoting) for a square matrix $\\overset{(m \\times m)}{\\Amtrx}$ and compute the $\\Pmtrx\\,\\Lmtrx\\,\\Umtrx$ factors. \n", "\n", "The factorization is obtained by elimination steps $k = 1,\\ldots,m-1$ so that $A^{(k+1)}_{i,j} = A^{(k)}_{i,j} - A^{(k)}_{k,j}\\, m_{i,k} \\ \\forall\\ i=k+1,\\ldots,m \\ \\text{and}\\ j=k,\\ldots,m$ where the multipliers $m_{i,k}$ are given by $m_{i,k} = \\frac{A^{(k)}_{i,k}}{A^{(k)}_{k,k}}$. When $k = m-1$, $A^{(m)}_{i,j}$, is upper triangular, that is, $U_{i,j} = A^{(m)}_{i,j}$ . The lower triangular matrix is obtained using the multipliers $m_{i,k}$, that is $L_{i,j} = m_{i,j} \\ \\forall \\ i>j$, $L_{i,i}=1$, and $L_{i,j}=0 \\ \\forall \\ i factorization\n", "The factors: $\\Pmtrx$, $\\Qmtrx$, $\\Lmtrx$, and $\\Umtrx$ where $\\Lmtrx\\,\\Umtrx = \\Pmtrx\\,\\Amtrx\\,\\Qmtrx$ can be obtained from a user-developed algorithm with **complete pivoting**. $\\Pmtrx$ is a row permutation matrix and $\\Qmtrx$ is a column permutation matrix, obtained by the underlying Gaussian elimination used to construct the $\\Lmtrx$ and $\\Umtrx$ factors.\n", "\n", "Program a $\\Lmtrx\\,\\Umtrx$ factorization algorithm (using complete pivoting) for a square matrix $\\overset{(m \\times m)}{\\Amtrx}$ and compute the $\\Pmtrx\\,\\Qmtrx\\,\\Lmtrx\\,\\Umtrx$ factors. The factorization is obtained by elimination steps $k = 1,\\ldots,m-1$ so that $A^{(k+1)}_{i,j} = A^{(k)}_{i,j} - A^{(k)}_{k,j}\\, m_{i,k} \\ \\forall\\ i=k+1,\\ldots,m \\ \\text{and}\\ j=k,\\ldots,m$ where the multipliers $m_{i,k}$ are given by $m_{i,k} = \\frac{A^{(k)}_{i,k}}{A^{(k)}_{k,k}}$. When $k = m-1$, $A^{(m)}_{i,j}$, is upper triangular, that is, $U_{i,j} = A^{(m)}_{i,j}$ . The lower triangular matrix is obtained using the multipliers $m_{i,k}$, that is $L_{i,j} = m_{i,j} \\ \\forall \\ i>j$, $L_{i,i}=1$, and $L_{i,j}=0 \\ \\forall \\ i