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

• Tutorial for ``scipy.linalg` <http://docs.scipy.org/doc/scipy/reference/tutorial/linalg.html>`__

Exact solution of linear system of equations

\begin{align} x + 2y &= 3 \\ 3x + 4y &= 17 \end{align}

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)

806 ms ± 9.5 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

%timeit la.inv(A1) @ b1

870 ms ± 9.07 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

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.

• How to interpret lapack function names

• Summary of BLAS functions

• Sumary of Lapack functions

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.557438524302

la.norm(C, 2) # largest sinular value

6.389028023601216

la.norm(C, -2) # smallest singular value

1.476590977094992

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)

Matrix is singular.

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)

Input a needs to be a square matrix.

x, resid, rank, s = la.lstsq(A2, b2)
x

array([ 0.41443305, -0.53024457,  0.61557362])

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.41443305, -0.53024457,  0.61557362])

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