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

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)









    



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



In [9]:

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



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

Trace



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

Determinant



In [17]:

    
la.det(C)









    Out[17]:





(-8-5j)

Rank



In [18]:

    
np.linalg.matrix_rank(C)









    Out[18]:





2

Norm



In [19]:

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









    Out[19]:





6.5574385243020004



In [20]:

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









    Out[20]:





6.3890280236012158



In [21]:

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









    Out[21]:





1.4765909770949921



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)









    



singular matrix



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)









    



expected square matrix



In [31]:

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









    Out[31]:





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.



In [33]:

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



In [34]:

    
least_squares(A2, b2)









    Out[34]:





array([ 0.4036226 ,  0.38604513,  0.40359296])

Matrix Decompositions



In [35]:

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









    Out[35]:





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

LU



In [36]:

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



In [37]:

    
p









    Out[37]:





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



In [38]:

    
l









    Out[38]:





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



In [39]:

    
u









    Out[39]:





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



In [40]:

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









    Out[40]:





True

Choleskey



In [41]:

    
U = la.cholesky(A)
U









    Out[41]:





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



In [42]:

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









    Out[42]:





True



In [43]:

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









    Out[43]:





True

QR



In [44]:

    
Q, R = la.qr(A)



In [45]:

    
Q









    Out[45]:





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



In [46]:

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









    Out[46]:





True



In [ ]:



In [47]:

    
np.allclose(Q@R, A)









    Out[47]:





True

Spectral



In [48]:

    
u, v = la.eig(A)



In [49]:

    
u









    Out[49]:





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



In [50]:

    
v









    Out[50]:





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



In [51]:

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









    Out[51]:





True



In [52]:

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









    Out[52]:





True

Inverting A



In [53]:

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









    Out[53]:





True

Powers of A



In [54]:

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









    Out[54]:





True

SVD



In [55]:

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



In [56]:

    
U









    Out[56]:





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



In [57]:

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









    Out[57]:





True



In [58]:

    
s









    Out[58]:





array([ 4.11554944,  0.88445056])



In [59]:

    
V









    Out[59]:





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



In [60]:

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









    Out[60]:





True



In [61]:

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









    Out[61]:





True



In [ ]: