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')

Resources

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