Linear Algebra Examples¶

This just shows the machanics of linear algebra calculations with python. See Lecture 5 for motivation and understanding.

import numpy as np
import scipy.linalg as la
import matplotlib.pyplot as plt
%matplotlib inline
plt.style.use('ggplot')

Exact solution of linear system of equations¶

A = np.array([[1,2],[3,4]])
A
array([[1, 2],
[3, 4]])
b = np.array([3,17])
b
array([ 3, 17])
x = la.solve(A, b)
x
array([ 11.,  -4.])
np.allclose(A @ x, b)
True
A1 = np.random.random((1000,1000))
b1 = np.random.random(1000)

Using solve is faster and more stable numerically than using matrix inversion¶

%timeit la.solve(A1, b1)
The slowest run took 5.24 times longer than the fastest. This could mean that an intermediate result is being cached
1 loops, best of 3: 90.3 ms per loop
%timeit la.inv(A1) @ b1
1 loops, best of 3: 140 ms per loop

Under the hood (Optional)¶

The solve function uses the dgesv fortran function to do the actual work. Here is an example of how to do this directly with the lapack function. There is rarely any reason to use blas or lapack functions directly becuase the linalg package provides more convenient functions that also perfrom error checking, but you can use Python to experiment with lapack or blas before using them in a language like C or Fortran.

import scipy.linalg.lapack as lapack
lu, piv, x, info = lapack.dgesv(A, b)
x
array([ 11.,  -4.])

Basic information about a matrix¶

C = np.array([[1, 2+3j], [3-2j, 4]])
C
array([[ 1.+0.j,  2.+3.j],
[ 3.-2.j,  4.+0.j]])
C.conjugate()
array([[ 1.-0.j,  2.-3.j],
[ 3.+2.j,  4.-0.j]])
def trace(M):
return np.diag(M).sum()
trace(C)
(5+0j)
np.allclose(trace(C), la.eigvals(C).sum())
True
la.det(C)
(-8-5j)
np.linalg.matrix_rank(C)
2
la.norm(C, None) # Frobenius (default)
6.5574385243020004
la.norm(C, 2) # largest sinular value
6.3890280236012158
la.norm(C, -2) # smallest singular value
1.4765909770949921
la.svdvals(C)
array([ 6.38902802,  1.47659098])

Least-squares solution¶

la.solve(A, b)
array([ 11.,  -4.])
x, resid, rank, s = la.lstsq(A, b)
x
array([ 11.,  -4.])
A1 = np.array([[1,2],[2,4]])
A1
array([[1, 2],
[2, 4]])
b1 = np.array([3, 17])
b1
array([ 3, 17])
try:
la.solve(A1, b1)
except la.LinAlgError as e:
print(e)
singular matrix
x, resid, rank, s = la.lstsq(A1, b1)
x
array([ 1.48,  2.96])
A2 = np.random.random((10,3))
b2 = np.random.random(10)
try:
la.solve(A2, b2)
except ValueError as e:
print(e)
expected square matrix
x, resid, rank, s = la.lstsq(A2, b2)
x
array([ 0.4036226 ,  0.38604513,  0.40359296])

Normal equations¶

One way to solve least squares equations $$X\beta = y$$ for $$\beta$$ is by using the formula $$\beta = (X^TX)^{-1}X^Ty$$ as you may have learnt in statistical theory classes (or can derive yourself with a bit of calculus). This is implemented below.

Note: This is not how the la.lstsq function solves least square problems as it can be inefficent for large matrices.

def least_squares(X, y):
return la.solve(X.T @ X, X.T @ y)
least_squares(A2, b2)
array([ 0.4036226 ,  0.38604513,  0.40359296])

Matrix Decompositions¶

A = np.array([[1,0.6],[0.6,4]])
A
array([[ 1. ,  0.6],
[ 0.6,  4. ]])

LU¶

p, l, u = la.lu(A)
p
array([[ 1.,  0.],
[ 0.,  1.]])
l
array([[ 1. ,  0. ],
[ 0.6,  1. ]])
u
array([[ 1.  ,  0.6 ],
[ 0.  ,  3.64]])
np.allclose(p@l@u, A)
True

Choleskey¶

U = la.cholesky(A)
U
array([[ 1.       ,  0.6      ],
[ 0.       ,  1.9078784]])
np.allclose(U.T @ U, A)
True
# If workiing wiht complex matrices
np.allclose(U.T.conj() @ U, A)
True

QR¶

Q, R = la.qr(A)
Q
array([[-0.85749293, -0.51449576],
[-0.51449576,  0.85749293]])
np.allclose((la.norm(Q[:,0]), la.norm(Q[:,1])), (1,1))
True
np.allclose(Q@R, A)
True

Spectral¶

u, v = la.eig(A)
u
array([ 0.88445056+0.j,  4.11554944+0.j])
v
array([[-0.98195639, -0.18910752],
[ 0.18910752, -0.98195639]])
np.allclose((la.norm(v[:,0]), la.norm(v[:,1])), (1,1))
True
np.allclose(v @ np.diag(u) @ v.T, A)
True

Inverting A¶

np.allclose(v @ np.diag(1/u) @ v.T, la.inv(A))
True

Powers of A¶

np.allclose(v @ np.diag(u**5) @ v.T, np.linalg.matrix_power(A, 5))
True

SVD¶

U, s, V = la.svd(A)
U
array([[ 0.18910752,  0.98195639],
[ 0.98195639, -0.18910752]])
np.allclose((la.norm(U[:,0]), la.norm(U[:,1])), (1,1))
True
s
array([ 4.11554944,  0.88445056])
V
array([[ 0.18910752,  0.98195639],
[ 0.98195639, -0.18910752]])
np.allclose((la.norm(V[:,0]), la.norm(V[:,1])), (1,1))
True
np.allclose(U @ np.diag(s) @ V, A)
True