Hacking Math I
Matrix Inverses and Numpy
This topic:¶
- Intro to Numpy and Scipy
- Inverse varieties and approaches
- Inverses and QR
- Polynomial interpolation
- Pseudo-inverses
Reading:
- I2ALA Chapter 11 (Matrix Inverses)
- https://docs.scipy.org/doc/numpy/user/quickstart.html
I. Intro to Numpy and SciPy¶
Leading Python tools:¶
Jupyter - "notebooks" for inline code + LaTex math + markup, etc.
NumPy - low-level array & matrix handling and algorithms
SciPy - higher level numerical algorithms (still fairly basic)
Matplotlib - matlab-style plotting & image display
SkLearn - (Scikit-learn) Machine Learning library
NumPy¶
- Numerical algorithm toolbox
- Similar to Matlab but with key differences, such as zero-based indexing (like C) instead of one-based (like math texts). Also Julia.
- Uses optimized linear algebra libraries called BLAS/LAPACK typically
- NumPy Arrays - special data structure which allows direct and efficient linear algebra manipulations (basically a list with added functionality).
- Used/extended in Tensorflow, Keras, pyTorch, Dask, ...
import numpy as np
x = np.array([[2,0,1,8],[1,2,3,4]])
print(x)
NumPy for Matlab Users¶
help(np.array)
"Vectorization"¶
Converting loops into linear algebra functions where loop is performed at lower level
E.g. instead of looping over a list to compute squares of elements, make into array and "square" array
v = [1,2,3,4,5]
v2 = []
for v_k in v:
v2.append(v_k**2)
v2
np.array(v)*2
np.array(v)**2
v+v2
np.array(v)+np.array(v2)
Matrix multiplication¶
multiplication of numpy arrays is elementwise by default
a = np.array([1,2,3])
b = np.array([3,3,3])
a*b
A = np.array([[1,2,3],[4,5,6]])
b = np.array([3,3,3])
A*b
A@b
np.matmul(A,b)
np.dot(A,b)
np.dot(A,A.T) # transpose also A.transpose()
Fast Numerical Mathematics¶
l = range(1234)
[i ** 2 for i in l]
l = range(1234)
%timeit [i ** 2 for i in l]
a = np.arange(1234)
a ** 2
a = np.arange(1234)
%timeit a ** 2
SciPy¶
Implements higher-level scientific algorithms using NumPy. Examples:
- Integration (scipy.integrate)
- Optimization (scipy.optimize)
- Interpolation (scipy.interpolate)
- Signal Processing (scipy.signal)
- Linear Algebra (scipy.linalg)
- Statistics (scipy.stats)
- File IO (scipy.io)
Also uses BLAS/LAPACK
import scipy
dir(scipy)
from scipy import stats
dir(scipy.stats)
dir(np.linalg)
dir(np.random)
II. Left and Right inverses¶
A = np.array([[[-3,-4],[4,6],[1,1]]])
B = 1/9*np.array([[-11,-10,16],[7,8,-11]])
B@A
C = 1/2*np.array([[0,-1,6],[0,1,-4]])
#np.matmul(C,A)
C@A
A = np.array([[[-3,-4],[4,6],[1,1]]])
B = 1/9*np.array([[-11,-10,16],[7,8,-11]])
C = 1/2*np.array([[0,-1,6],[0,1,-4]])
print(B@A)
print(C@A)
b = np.array([1,-2,0])
print(B@b)
print(C@b)
x = B@b
A@x
np.linalg.inv(A)
A = np.array([[-3,-4],[4,6],[1,1]])
B = 1/9*np.array([[-11,-10,16],[7,8,-11]])
C = 1/2*np.array([[0,-1,6],[0,1,-4]])
At = A.transpose()
Bt = B.transpose()
Ct = C.transpose()
print(At@Bt)
print(At@Ct)
b = np.array([1,2])
print(Bt@b)
print(Ct@b)
xB = Bt@b
xC = Ct@b
print(At@xB)
print(At@xC)
print(At.shape)
print(Bt.shape)
III. "Inverses"¶
A = np.array([[1,-2,3],[0,2,2],[-3,-4,-4]])
Ainv = 1/30*np.array([[0,-20,-10],[-6,5,-2],[6,10,2]])
print(A@Ainv)
print(Ainv@A)
print(np.linalg.inv(A))
print(Ainv)
II.A. Inverse Examples¶
Diagonal matrices¶
Identity Matrix
$I_3 = \begin{bmatrix} 1 & 0 & 0\\ 0 & 1 & 0 \\ 0 & 0 & 1\\ \end{bmatrix}$
Diagonal matrix with nonzero diagonal elements $D = \begin{bmatrix} 2 & 0 & 0\\ 0 & 3 & 0 \\ 0 & 0 & 4\\ \end{bmatrix}$
Diagonal matrix with zeros on the diagonal $D = \begin{bmatrix} 2 & 0 & 0\\ 0 & 0 & 0 \\ 0 & 0 & 4\\ \end{bmatrix}$
Use simplicity of matrix multiplication for diagonal matrices to get inverse
"Orthogonal matrices"¶
- square matrix with othonormal columns $Q^TQ = I$ so $Q^{-1} = Q^T$
- should be called "orthonormal matrix"
- Generalized rotation and/or reflection matrix
Apply definition of orthonormal to columns to prove inverse
III. Inverses and QR¶
import numpy as np
A = np.array([[1,2,3],[3,1,5],[1,7,8]])
b = np.array([5,6,7])
# solve Ax = b
x1 = np.linalg.inv(A)@b
print(x1)
# use QR
(Q,R) = np.linalg.qr(A)
print(Q)
print(R)
Q@Q.transpose()
x2 = np.linalg.inv(R)@(Q.transpose()@b)
print(x2)
IV. Polynomial interpolation¶
V. Pseudoinverse¶
