Linear Algebra Review and SciPy Basics
Contents
6.1. Linear Algebra Review and SciPy Basics#
Screenshots from Section 2.2 Vectors and Matrices in Biegler (2010)
Recommended Links
SciPy Lecture Notes (Especially 1.1 - 1.6 for everyone and 2.1 - 2.2 for advanced users)
Tip: We will mostly use SciPy for linear algebra. (It has more sophisticated capabilities than NumPy.) We will use NumPy if/when SciPy does not offer a specific command.
# Load required Python libraries.
import matplotlib.pyplot as plt
import numpy as np
from scipy import linalg
6.1.1. Notation (for textbook)#
Scalar constants and variables - Greek letters
Vectors - Lower case Roman letters
Matrices - Upper case Roman letters
The notation for lectures will (hopefully) be obvious from context. I sometimes use Greek letters for constants and Roman letters for variables… unless the Greek letters have engineering or scientific meaning.
6.1.2. Matrix Operations#
# Generate random matrices
n = 2
m = 3
Ao = np.random.rand(n,m)
Bo = np.random.rand(m,n)
print("Ao = ")
print(Ao)
print("\nBo = ")
print(Bo)
# What is the dimension of A*B? Try to answer this without the computer.
# Calculate A*B using Python
print("\nAo*Bo =")
print(Ao.dot(Bo))
Ao =
[[0.90471685 0.17317304 0.40300532]
[0.63393766 0.90146479 0.61500922]]
Bo =
[[0.34621818 0.9393221 ]
[0.67416793 0.95718324]
[0.61310823 0.40727659]]
Ao*Bo =
[[0.67706301 1.17971349]
[1.20428661 1.7088175 ]]
Activity
Calculate \(A_o^{T}\)
Create a square matrix and check if it is symmetric
Create a 3x3 identity matrix
# Transpose
print("transpose(Ao) = \n",Ao.transpose(),"\n")
# Create a square matrix. Is it symmetric?
# Identify matrix
print("I = \n",np.identity(3),"\n")
transpose(Ao) =
[[0.90471685 0.63393766]
[0.17317304 0.90146479]
[0.40300532 0.61500922]]
I =
[[1. 0. 0.]
[0. 1. 0.]
[0. 0. 1.]]
6.1.3. Determinant#
![Book]../../media/det.png)
# Generate random matrices
nd = 2
Ad = np.random.rand(nd,nd)
Bd = np.random.rand(nd,nd)
Cd = np.array([(1, 0, 0),(-2, 0, 0),(0, 1, 1)])
print("Ad = ")
print(Ad)
print("\nBd = ")
print(Bd)
print("\nCd = ")
print(Cd)
Ad =
[[0.35144304 0.40331999]
[0.98819481 0.83989719]]
Bd =
[[0.24359748 0.02035267]
[0.32634249 0.99856976]]
Cd =
[[ 1 0 0]
[-2 0 0]
[ 0 1 1]]
Activity
Verify the properties listed under Some properties for the determinant include
If \(A_d\), \(B_d\) or \(C_d\) singular?
# Property 1
print("det(A*B) = ",linalg.det(Ad.dot(Bd)))
print("det(A)*det(B) = ",linalg.det(Ad)*linalg.det(Bd))
det(A*B) = -0.024461085784763897
det(A)*det(B) = -0.02446108578476388
# Property 2
# Property 3
# Property 4
Another Example (for at home)
Below is an example from Linear Algebra and Its Applications by David C. Lay, 3rd edition that shows how to calculate a determinant by hand.
6.1.4. Rank#
Activity
Predict the rank of \(A_o^T\), \(B_o^T\), \(A_d\), \(B_d\), \(C_d\)
Calculate np.linalg.matrix_rank(). Where you correct?
print("rank(A_o) = ",np.linalg.matrix_rank(Ao),"\n")
# Fill in remainder here.
rank(A_o) = 2
6.1.5. Inverse#
Activity
Calculate the inverse of \(A_o\) and verity that \(A_o^{-1} A_o = I\)
Verify that det(\(A^{-1}\)) = 1/det(\(A\))
Is \(Q\) an orthogonal matrix?
# Task 1
# Task 2
# Task 3
# Is this an orthogonal matrix? Yes, no, or need more information?
theta = 0.2
Q = np.array([(np.cos(theta), -np.sin(theta)),(np.sin(theta),np.cos(theta))])
print("Q = ")
print(Q)
Q =
[[ 0.98006658 -0.19866933]
[ 0.19866933 0.98006658]]
6.1.6. Solving Linear Systems#
Consider the linear system \(A_l x = b_l\)
Al = np.array([(4,3),(6,3)])
bl = np.array([1,0])
print("Al = \n",Al)
print("\nbl = \n",bl)
Al =
[[4 3]
[6 3]]
bl =
[1 0]
6.1.6.1. Explicit Inverse#
Calculate \(x\) by explicitly using \(A_l^{-1}\). Hint: Use linalg.inv().
Ainv = linalg.inv(Al)
print("Ainv = \n",Ainv,"\n")
# Now calculate and print x
Ainv =
[[-0.5 0.5 ]
[ 1. -0.66666667]]
6.1.6.2. LU Decomposition#
Perform LU decomposition on \(A_l\). What structures do the \(P\), \(L\), and \(U\) matrices have?
# Create test linear systems
A = np.random.rand(4,4)
print("A = \n",A,"\n")
b = np.random.rand(4,1)
print("b = \n",b,"\n")
A =
[[0.68331246 0.29422836 0.22419151 0.41620186]
[0.05988483 0.58705665 0.07246852 0.45441375]
[0.6024749 0.36997327 0.13872522 0.45388322]
[0.34311843 0.47101812 0.51826738 0.25525007]]
b =
[[0.00949348]
[0.82027783]
[0.90273651]
[0.37915841]]
# Perform LU decomposition
(P, L, U) = linalg.lu(A)
# Permutation matrix
print("P = \n",P)
# Lower diagonal matrix
print("L = \n",L)
# Upper diagonal matrix
print("U = \n",U)
# Verify result
print("P*L*U = \n",P.dot(L.dot(U)),"\n")
P =
[[1. 0. 0. 0.]
[0. 1. 0. 0.]
[0. 0. 0. 1.]
[0. 0. 1. 0.]]
L =
[[ 1. 0. 0. 0. ]
[ 0.08763901 1. 0. 0. ]
[ 0.50213986 0.5759686 1. 0. ]
[ 0.88169751 0.19696885 -0.18479522 1. ]]
U =
[[ 0.68331246 0.29422836 0.22419151 0.41620186]
[ 0. 0.56127076 0.05282059 0.41793823]
[ 0. 0. 0.37526888 -0.19446077]
[ 0. 0. 0. -0.03133716]]
P*L*U =
[[0.68331246 0.29422836 0.22419151 0.41620186]
[0.05988483 0.58705665 0.07246852 0.45441375]
[0.6024749 0.36997327 0.13872522 0.45388322]
[0.34311843 0.47101812 0.51826738 0.25525007]]
6.1.6.2.1. Is P orthogonal?#
print("inv(P) = \n",linalg.inv(P),"\n")
print("transpose(P) = \n",P.T,"\n")
print("P.T*P = \n",np.matmul(P.T,P),"\n")
inv(P) =
[[ 1. -0. -0. 0.]
[ 0. 1. -0. 0.]
[ 0. 0. -0. 1.]
[ 0. 0. 1. 0.]]
transpose(P) =
[[1. 0. 0. 0.]
[0. 1. 0. 0.]
[0. 0. 0. 1.]
[0. 0. 1. 0.]]
P.T*P =
[[1. 0. 0. 0.]
[0. 1. 0. 0.]
[0. 0. 1. 0.]
[0. 0. 0. 1.]]
Yes! We see that \(P^T = P^{-1}\) and \(P^T \cdot P = I\).
6.1.6.2.2. MATLAB#
Defines LU decomposition as follows:
Consider $\(A \cdot x=b\)$
Using MATLAB LU definition
Substitute into linear system:
Let \(y = U \cdot x\) and substitute.
Step 1. Solve \(L \cdot y = P \cdot b\) for \(y\)
Step 2. Solve \(U \cdot x = y\) for \(x\)
Pb = P.dot(b)
print("P*b = ",Pb,"\n")
yLU = linalg.solve(L,Pb)
print("yLU = ",yLU,"\n")
xLU = linalg.solve(U,yLU)
print("xLU = ",xLU,"\n")
P*b = [[0.00949348]
[0.82027783]
[0.37915841]
[0.90273651]]
yLU = [[ 0.00949348]
[ 0.81944583]
[-0.09758371]
[ 0.71492782]]
xLU = [[ 9.4407442 ]
[ 19.58501477]
[-12.08206663]
[-22.81406129]]
6.1.6.2.3. SciPy#
Defines LU decomposition as:
Consider $\(A \cdot x=b\)$
SciPy LU definition
Substitute into linear system:
Let \(y = U \cdot x\) and substitute.
Step 1. Solve \(L \cdot y = P^T \cdot b\) for \(y\)
Step 2. Solve \(U \cdot x = y\) for \(x\)
## LU decomposition algorithm.
Pb = (P.T).dot(b)
print("P.T*b = ",Pb,"\n")
yLU = linalg.solve(L,Pb)
print("yLU = ",yLU,"\n")
xLU = linalg.solve(U,yLU)
print("xLU = ",xLU,"\n")
P.T*b = [[0.00949348]
[0.82027783]
[0.37915841]
[0.90273651]]
yLU = [[ 0.00949348]
[ 0.81944583]
[-0.09758371]
[ 0.71492782]]
xLU = [[ 9.4407442 ]
[ 19.58501477]
[-12.08206663]
[-22.81406129]]
What is the difference? A transpose.
6.1.6.2.4. Verify our answer with linalg.solve#
Solve the linear system using linalg.solve
x = linalg.solve(Al,bl)
print("x = \n",x)
x =
[-0.5 1. ]
6.1.7. Invertable Matrix Theorem#
Screenshots from Linear Algebra and Its Applications by David C. Lay, 3rd edition
6.1.8. Eigenvectors and Eigenvalues#
Ae = np.array([(0, -4, -6), (-1, 0, -3), (1, 2, 5)])
print("Ae = \n",Ae)
Ae =
[[ 0 -4 -6]
[-1 0 -3]
[ 1 2 5]]
Activity
Do you expect this matrix to have negative eigenvalues? Choose: Yes / No / Cannot tell from inspection.
Calculate the eigenvalues.
Based on your calculations above, is this matrix i) positive definite, ii) positive semi definite, iii) negative definite, iv) negative semi definite, v) indefinite, or vi) cannot say without additional calculations? Write a sentence to justify your answer.
Calculate the determinate using only the eigenvalues.
# Expect negative eigenavalues?
# Yes / No / Cannot tell from inspect?
# Write a sentence to justify your answer.
### Calculate eigenvalues
# Matrix Al
print("Matrix = \n",Al,"\n")
l, v = linalg.eig(Al)
print("Eigenvalues = ",l,"\n")
print("Eigenvectors = \n",v,"\n")
Matrix =
[[4 3]
[6 3]]
Eigenvalues = [ 7.77200187+0.j -0.77200187+0.j]
Eigenvectors =
[[ 0.62246561 -0.53222953]
[ 0.78264715 0.8466001 ]]
6.1.9. Singular Value Decomposition#
Notes will be given in class.
Activity
Calculate singular value decomposition of \(A_l\) and \(C_d\)
What is the rank of each matrix?
# Matrix Al
print("Matrix = \n",Al,"\n")
U,s,Vh = linalg.svd(Al)
print("U = \n",U,"\n")
print("S = \n",s,"\n")
print("Vh = \n",Vh,"\n")
## Rank (from inspecting singular values)?
Matrix =
[[4 3]
[6 3]]
U =
[[-0.59581566 -0.80312122]
[-0.80312122 0.59581566]]
S =
[8.33557912 0.71980602]
Vh =
[[-0.86400595 -0.50348159]
[ 0.50348159 -0.86400595]]
6.1.10. Vector and Matrix Norms#
x = np.random.rand(3,1)
print("x = \n",x)
y = np.random.rand(3,1)
print("y = \n",y)
x =
[[0.02808233]
[0.20488865]
[0.77425602]]
y =
[[0.73113614]
[0.20737529]
[0.97754014]]
Activity
Verify the vector norm properties hold for the p-2 norm using \(x\) and \(y\). Hint: linalg.norm()
Calculate \(||A_l||_1\), \(||A_l||_{\infty}\) and \(||A_l||_2\) using both the rules given above and SciPy
6.1.11. Condition Number#
Activity
Estimate the condition number of \(A_l\) by inspecting the SVD results.
Calculate the condition number of \(A_l\) using np.linalg.cond()