Software tutorial/Matrix operations

LU decomposition

Let \(A = \left(\begin{array}{cccc} 4 & 0 & 1 \\ 0 & 6 & 2 \\ 3.2 & 7 & 4 \end{array} \right)\) and let \(b = \left(\begin{array}{cccc} 9.1 \\ 2 \\ 0 \end{array} \right) \)

MATLAB

Python

A = [4, 0, 1; 0, 6, 2; 3.2, 7, 4]; % LU decomposition of A [L,U,P] = lu(A) L = 1.0000 0 0 0.8000 1.0000 0 0 0.8571 1.0000 U = 4.0000 0 1.0000 0 7.0000 3.2000 0 0 -0.7429 P = 1 0 0 0 0 1 0 1 0 The P matrix shows that rows 2 and rows 3 were interchanged during the Gauss elimination pivoting steps.

import numpy as np from scipy.linalg import * A = np.array([[4, 0, 1],[0, 6, 2],[3.2, 7, 4]]) LU, P = lu_factor(A) >>> LU array([[ 4. , 0. , 1. ], [ 0.8 , 7. , 3.2 ], [ 0. , 0.85714286, -0.74285714]]) >>> P array([0, 2, 2]) The LU matrices are overlayed in the output. Notice that the results are the same as MATLAB's, except we do not need to store the diagonal set of ones from the L matrix. The P vector is also a bit cryptic - it should be read from left to right, in order. The way to read it is that row \(i\) of matrix \(A\) was interchanged with row P[i] in \(A\), where i is the index into vector P. row 0 in \(A\) was interchanged with row 0: i.e. no pivoting for the first row row 1 in \(A\) was interchanged with row 2: i.e. rows 2 and 3 in \(A\) were swapped row 2 in \(A\) was interchanged with row 2: i.e. no further pivoting.

Once we have \(L, U\) and \(P\) we can solve the system of equations for a given vector, \(b\). It takes a bit more work to do this in MATLAB than with Python. Recall from class that \(A = LU\), but with a permutation or pivot matrix, it becomes \(PA = LU\), or alternatively: \(A = P^{-1}LU\)

\[ \begin{align} Ax &= b \\ P^{-1}LUx &= b \\ \text{Let}\qquad Ux &= y \\ \text{Then solve}\qquad \left(P^{-1}L\right)y &= b \qquad \text{for}\,\, y \\ \text{Finally, solve}\qquad Ux &= y \qquad \text{for}\,\, x \end{align} \]

MATLAB	Python
b = [9.1, 2, 0]'; [L,U,P] = lu(A); y = mldivide(inv(P)*L, b); x = mldivide(U, y); >> x x = 5.0481 4.0308 -11.0923 % Verify the solution: >> A*x - b ans = 1.0e-14 * 0.1776 0 0 The small values of the solution error indicate that our LU solution is accurate. You should compare the LU solution with solution from mldivide(A, b).	LU, P = lu_factor(A) b = np.array([9.1, 2, 0]) x = lu_solve((LU,P), b) >>> x array([ 5.04807692, 4.03076923, -11.09230769]) # Verify the solution: >>> np.dot(A, x) - b array([ 1.77635684e-15, 0.00000000e+00, 0.00000000e+00])

Matrix determinant

One way to test if a matrix is invertable is to calculate the determinant. If it is non-zero, then \(A\) is invertable.

Let \(A = \left(\begin{array}{cccc} 4 & 0 & 1 \\ 0 & 6 & 2 \\ 3.2 & 7 & 4 \end{array} \right)\)

MATLAB	Python
A = [4, 0, 1; 0, 6, 2; 3.2, 7, 4]; >> det(A) ans = 20.8000	import numpy as np from scipy.linalg import * A = np.array([[4, 0, 1],[0, 6, 2],[3.2, 7, 4]]) >>> det(A) 20.800000000000011

Matrix inverse

The matrix inverse is defined so that \(A A^{-1} = I\). Let \(A = \left(\begin{array}{cccc} 4 & 0 & 1 \\ 0 & 6 & 2 \\ 3.2 & 7 & 4 \end{array} \right)\)

MATLAB

Python

A = [4, 0, 1; 0, 6, 2; 3.2, 7, 4]; >> inv(A) ans = 0.4808 0.3365 -0.2885 0.3077 0.6154 -0.3846 -0.9231 -1.3462 1.1538 >> A * inv(A) ans = 1.0000 0 0 -0.0000 1.0000 0.0000 0 0 1.0000

import numpy as np from scipy.linalg import * A = np.array([[4, 0, 1],[0, 6, 2],[3.2, 7, 4]]) >>> inv(A) array([[ 0.48076923, 0.33653846, -0.28846154], [ 0.30769231, 0.61538462, -0.38461538], [-0.92307692, -1.34615385, 1.15384615]]) >>> np.dot(A, inv(A)) array([[ 1.00000000e+00, 0.00000000e+00, 0.00000000e+00], [ -2.22044605e-16, 1.00000000e+00, 4.44089210e-16], [ 0.00000000e+00, 0.00000000e+00, 1.00000000e+00]]) Notice the small, essentially zero elements in the off-diagonal elements