Linear Algebra Examples¶

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

In [1]:

import numpy as np
import scipy.linalg as la
import matplotlib.pyplot as plt
%matplotlib inline

In [2]:

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}

In [3]:

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

Out[3]:

array([[1, 2],
       [3, 4]])

In [4]:

b = np.array([3,17])
b

Out[4]:

array([ 3, 17])

In [5]:

x = la.solve(A, b)
x

Out[5]:

array([11., -4.])

In [6]:

np.allclose(A @ x, b)

Out[6]:

True

In [7]:

A1 = np.random.random((1000,1000))
b1 = np.random.random(1000)

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

In [8]:

%timeit la.solve(A1, b1)

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

In [9]:

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

In [10]:

import scipy.linalg.lapack as lapack

In [11]:

lu, piv, x, info = lapack.dgesv(A, b)
x

Out[11]:

array([11., -4.])

Basic information about a matrix¶

In [12]:

C = np.array([[1, 2+3j], [3-2j, 4]])
C

Out[12]:

array([[1.+0.j, 2.+3.j],
       [3.-2.j, 4.+0.j]])

In [13]:

C.conjugate()

Out[13]:

array([[1.-0.j, 2.-3.j],
       [3.+2.j, 4.-0.j]])

In [14]:

def trace(M):
    return np.diag(M).sum()

In [15]:

trace(C)

Out[15]:

(5+0j)

In [16]:

np.allclose(trace(C), la.eigvals(C).sum())

Out[16]:

True

In [17]:

la.det(C)

Out[17]:

(-8-5j)

In [18]:

np.linalg.matrix_rank(C)

Out[18]:

In [19]:

la.norm(C, None) # Frobenius (default)

Out[19]:

6.557438524302

In [20]:

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

Out[20]:

6.389028023601216

In [21]:

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

Out[21]:

1.476590977094992

In [22]:

la.svdvals(C)

Out[22]:

array([6.38902802, 1.47659098])

Least-squares solution¶

In [23]:

la.solve(A, b)

Out[23]:

array([11., -4.])

In [24]:

x, resid, rank, s = la.lstsq(A, b)
x

Out[24]:

array([11., -4.])

In [25]:

A1 = np.array([[1,2],[2,4]])
A1

Out[25]:

array([[1, 2],
       [2, 4]])

In [26]:

b1 = np.array([3, 17])
b1

Out[26]:

array([ 3, 17])

In [27]:

try:
    la.solve(A1, b1)
except la.LinAlgError as e:
    print(e)

Matrix is singular.

In [28]:

x, resid, rank, s = la.lstsq(A1, b1)
x

Out[28]:

array([1.48, 2.96])

In [29]:

A2 = np.random.random((10,3))
b2 = np.random.random(10)

In [30]:

try:
    la.solve(A2, b2)
except ValueError as e:
    print(e)

Input a needs to be a square matrix.

In [31]:

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

Out[31]:

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.

In [32]:

def least_squares(X, y):
    return la.solve(X.T @ X, X.T @ y)

In [33]:

least_squares(A2, b2)

Out[33]:

array([ 0.41443305, -0.53024457,  0.61557362])

Matrix Decompositions¶

In [34]:

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

Out[34]:

array([[1. , 0.6],
       [0.6, 4. ]])

LU¶

In [35]:

p, l, u = la.lu(A)

In [36]:

Out[36]:

array([[1., 0.],
       [0., 1.]])

In [37]:

Out[37]:

array([[1. , 0. ],
       [0.6, 1. ]])

In [38]:

Out[38]:

array([[1.  , 0.6 ],
       [0.  , 3.64]])

In [39]:

np.allclose(p@l@u, A)

Out[39]:

True

Choleskey¶

In [40]:

U = la.cholesky(A)
U

Out[40]:

array([[1.       , 0.6      ],
       [0.       , 1.9078784]])

In [41]:

np.allclose(U.T @ U, A)

Out[41]:

True

In [42]:

# If workiing wiht complex matrices
np.allclose(U.T.conj() @ U, A)

Out[42]:

True

QR¶

In [43]:

Q, R = la.qr(A)

In [44]:

Out[44]:

array([[-0.85749293, -0.51449576],
       [-0.51449576,  0.85749293]])

In [45]:

np.allclose((la.norm(Q[:,0]), la.norm(Q[:,1])), (1,1))

Out[45]:

True

In [ ]:

In [46]:

np.allclose(Q@R, A)

Out[46]:

True

Spectral¶

In [47]:

u, v = la.eig(A)

In [48]:

Out[48]:

array([0.88445056+0.j, 4.11554944+0.j])

In [49]:

Out[49]:

array([[-0.98195639, -0.18910752],
       [ 0.18910752, -0.98195639]])

In [50]:

np.allclose((la.norm(v[:,0]), la.norm(v[:,1])), (1,1))

Out[50]:

True

In [51]:

np.allclose(v @ np.diag(u) @ v.T, A)

Out[51]:

True

Inverting A¶

In [52]:

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

Out[52]:

True

Powers of A¶

In [53]:

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

Out[53]:

True

SVD¶

In [54]:

U, s, V = la.svd(A)

In [55]:

Out[55]:

array([[ 0.18910752,  0.98195639],
       [ 0.98195639, -0.18910752]])

In [56]:

np.allclose((la.norm(U[:,0]), la.norm(U[:,1])), (1,1))

Out[56]:

True

In [57]:

Out[57]:

array([4.11554944, 0.88445056])

In [58]:

Out[58]:

array([[ 0.18910752,  0.98195639],
       [ 0.98195639, -0.18910752]])

In [59]:

np.allclose((la.norm(V[:,0]), la.norm(V[:,1])), (1,1))

Out[59]:

True

In [60]:

np.allclose(U @ np.diag(s) @ V, A)

Out[60]:

True

In [ ]: