In [3]:
import numpy as np
import numpy.random as npr
import pylab as pl
%matplotlib inline

In [2]:
a = 1 + 5j

In [3]:
a


Out[3]:
(1+5j)

In [4]:
x = 1

In [25]:
a = np.complex(x,x+2)

In [26]:
a


Out[26]:
(1+3j)

In [74]:
np.complex(*npr.rand(2))


Out[74]:
(0.6874750903219116+0.49300799260849226j)

In [75]:
rand_complex = lambda n=1 : np.array([np.complex(*npr.randn(2)) for __ in xrange(n)])

In [87]:
n = 100 # samples
p = 5 # dimension
X = np.array([rand_complex(p) for __ in xrange(n)])

In [88]:
X.shape


Out[88]:
(100, 5)

In [89]:
a.conjugate()*(a)


Out[89]:
(10+0j)

In [90]:
a*a.conjugate()


Out[90]:
(10+0j)

In [91]:
X[0]*(X[0].conjugate())/n


Out[91]:
array([ 0.00622855+0.j,  0.00153904+0.j,  0.00072866+0.j,  0.00143073+0.j,
        0.08361050+0.j])

In [92]:
X.T.shape


Out[92]:
(5, 100)

In [93]:
x = X.T[0]

In [98]:
(x.dot((x.T).conjugate()))


Out[98]:
(168.37844244984677+0j)

In [95]:
# Joel Tropp: Introduction to Matrix Concentration inequalities says on p.6 that
# this is an unbiased estimator for the covariance matrix of X, but I don't understand:
Y = sum([x*x.conjugate() for x in X.T]) / n

In [ ]:
Y = np.zeros((p,p))
for i in range(p):
    Y[i] =

In [58]:
Y.shape


Out[58]:
(100,)

In [59]:
Y


Out[59]:
array([ 0.18102339+0.j,  0.44827120+0.j,  0.17766590+0.j,  0.28715848+0.j,
        0.11166032+0.j,  0.35993556+0.j,  0.17921374+0.j,  0.18321532+0.j,
        0.18539421+0.j,  0.20290436+0.j,  0.15008532+0.j,  0.24967354+0.j,
        0.25000338+0.j,  0.13448715+0.j,  0.15879320+0.j,  0.17102191+0.j,
        0.16692174+0.j,  0.15227429+0.j,  0.14850775+0.j,  0.22891312+0.j,
        0.15522525+0.j,  0.14109556+0.j,  0.12357581+0.j,  0.19688546+0.j,
        0.25709227+0.j,  0.16899224+0.j,  0.16745321+0.j,  0.14542843+0.j,
        0.15685478+0.j,  0.08720036+0.j,  0.19738370+0.j,  0.18705173+0.j,
        0.20032441+0.j,  0.20302085+0.j,  0.25405256+0.j,  0.15419715+0.j,
        0.23919532+0.j,  0.23659836+0.j,  0.17278417+0.j,  0.28279633+0.j,
        0.25462865+0.j,  0.19720029+0.j,  0.09808417+0.j,  0.19592428+0.j,
        0.24245900+0.j,  0.19234139+0.j,  0.29313953+0.j,  0.33177463+0.j,
        0.26364197+0.j,  0.13769819+0.j,  0.09364633+0.j,  0.15144355+0.j,
        0.14509015+0.j,  0.12930946+0.j,  0.24461192+0.j,  0.11364361+0.j,
        0.10633440+0.j,  0.26226123+0.j,  0.21427967+0.j,  0.14398721+0.j,
        0.21803783+0.j,  0.12282650+0.j,  0.20879291+0.j,  0.20772165+0.j,
        0.24455720+0.j,  0.17584388+0.j,  0.10433116+0.j,  0.23035885+0.j,
        0.19758060+0.j,  0.12926350+0.j,  0.11166044+0.j,  0.23388821+0.j,
        0.18382383+0.j,  0.19361171+0.j,  0.21246544+0.j,  0.18450932+0.j,
        0.20414653+0.j,  0.26029765+0.j,  0.28176651+0.j,  0.24083597+0.j,
        0.26683940+0.j,  0.12463813+0.j,  0.35934650+0.j,  0.13469039+0.j,
        0.24455330+0.j,  0.15415118+0.j,  0.11207082+0.j,  0.14338405+0.j,
        0.14644119+0.j,  0.18799786+0.j,  0.09356931+0.j,  0.23407897+0.j,
        0.19657983+0.j,  0.19050683+0.j,  0.22846786+0.j,  0.14180992+0.j,
        0.09258563+0.j,  0.13536979+0.j,  0.25461074+0.j,  0.15827486+0.j])

In [56]:
np.cov(X.T).shape


Out[56]:
(10, 10)

In [61]:
np.cov(X.T)


Out[61]:
array([[  1.69921994e+00+0.j        ,   8.92856574e-02+0.01686832j,
          1.34784793e-01+0.158611j  ,  -7.42028633e-02-0.19772521j,
         -1.96657398e-01-0.05822664j,   3.19507912e-02+0.07432521j,
          9.27351924e-02-0.04605331j,   1.63933659e-01-0.14195614j,
         -6.60988819e-04+0.16464858j,  -1.83180013e-01+0.02376909j],
       [  8.92856574e-02-0.01686832j,   2.25555608e+00+0.j        ,
          8.86498761e-02+0.00172851j,  -2.61418922e-02-0.13402592j,
          8.14839818e-02+0.25809225j,   1.04613341e-01+0.08052215j,
         -6.35370435e-02+0.3208301j ,   4.87101237e-03+0.28323388j,
          3.36270444e-02+0.06370114j,   1.11268402e-01+0.09221383j],
       [  1.34784793e-01-0.158611j  ,   8.86498761e-02-0.00172851j,
          1.69321628e+00+0.j        ,  -1.43539098e-01+0.28361066j,
          2.44214620e-01-0.21974763j,   1.86281050e-01-0.0354298j ,
         -2.49927279e-01-0.0344726j ,  -1.80497851e-01+0.16198084j,
          8.73404300e-03+0.07984407j,   1.11319606e-01-0.01421182j],
       [ -7.42028633e-02+0.19772521j,  -2.61418922e-02+0.13402592j,
         -1.43539098e-01-0.28361066j,   1.90239623e+00+0.j        ,
         -1.36673552e-02-0.17075837j,  -2.67402276e-02-0.04381924j,
         -7.71042074e-02-0.01403706j,   1.08328786e-01+0.06787448j,
         -2.83654076e-02+0.21684511j,   2.80496583e-02+0.08129267j],
       [ -1.96657398e-01+0.05822664j,   8.14839818e-02-0.25809225j,
          2.44214620e-01+0.21974763j,  -1.36673552e-02+0.17075837j,
          2.41594657e+00+0.j        ,  -6.84852807e-02-0.04017138j,
          1.32816749e-01+0.02558448j,  -2.69296322e-01-0.04357293j,
          1.12562803e-01+0.02908231j,   2.07439846e-02-0.01222748j],
       [  3.19507912e-02-0.07432521j,   1.04613341e-01-0.08052215j,
          1.86281050e-01+0.0354298j ,  -2.67402276e-02+0.04381924j,
         -6.84852807e-02+0.04017138j,   1.82254980e+00+0.j        ,
         -3.41275973e-02+0.00388601j,   1.56677769e-02-0.11289085j,
          1.28775659e-01+0.00936737j,   1.81752030e-01-0.06974281j],
       [  9.27351924e-02+0.04605331j,  -6.35370435e-02-0.3208301j ,
         -2.49927279e-01+0.0344726j ,  -7.71042074e-02+0.01403706j,
          1.32816749e-01-0.02558448j,  -3.41275973e-02-0.00388601j,
          2.05579563e+00+0.j        ,   3.35091176e-01-0.06164157j,
         -1.08200265e-01-0.23249233j,   6.24350565e-02-0.28277935j],
       [  1.63933659e-01+0.14195614j,   4.87101237e-03-0.28323388j,
         -1.80497851e-01-0.16198084j,   1.08328786e-01-0.06787448j,
         -2.69296322e-01+0.04357293j,   1.56677769e-02+0.11289085j,
          3.35091176e-01+0.06164157j,   1.92869961e+00+0.j        ,
          1.67929241e-01-0.07928676j,   3.56160070e-02-0.13114128j],
       [ -6.60988819e-04-0.16464858j,   3.36270444e-02-0.06370114j,
          8.73404300e-03-0.07984407j,  -2.83654076e-02-0.21684511j,
          1.12562803e-01-0.02908231j,   1.28775659e-01-0.00936737j,
         -1.08200265e-01+0.23249233j,   1.67929241e-01+0.07928676j,
          1.77929773e+00+0.j        ,  -3.33767487e-02+0.32536837j],
       [ -1.83180013e-01-0.02376909j,   1.11268402e-01-0.09221383j,
          1.11319606e-01+0.01421182j,   2.80496583e-02-0.08129267j,
          2.07439846e-02+0.01222748j,   1.81752030e-01+0.06974281j,
          6.24350565e-02+0.28277935j,   3.56160070e-02+0.13114128j,
         -3.33767487e-02-0.32536837j,   1.73280727e+00+0.j        ]])

In [2]:
from scipy import linalg

Quick review of scipy.linalg and numpy.linalg


In [4]:
a = npr.rand(100,100)
b = npr.rand(1000,1000)

Inverting a matrix

Find the inverse of $\mathbf{A}$, $\mathbf{A}^{-1}$, such that $$\mathbf{A} \mathbf{A}^{-1} = \mathbf{I}$$

  • numpy.linalg seems significantly faster than scipy.linalg for large and small matrices, but scipy.linalg is faster for matrices aroud $100 \times 100$

In [5]:
tiny = npr.rand(10,10)

In [118]:
%timeit linalg.inv(tiny)


10000 loops, best of 3: 29.4 µs per loop

In [119]:
%timeit np.linalg.inv(tiny)


100000 loops, best of 3: 12.5 µs per loop

In [115]:
%timeit ans = linalg.inv(a)


1000 loops, best of 3: 270 µs per loop

In [107]:
%timeit ans = linalg.inv(b)


10 loops, best of 3: 59.1 ms per loop

In [27]:
c = npr.rand(5000,5000)

In [108]:
%timeit ans = linalg.inv(c)


1 loops, best of 3: 6.96 s per loop

In [113]:
%timeit np.linalg.inv(a)


1000 loops, best of 3: 347 µs per loop

In [110]:
%timeit ans = np.linalg.inv(b)


10 loops, best of 3: 37.7 ms per loop

In [111]:
%timeit ans = np.linalg.inv(c)


1 loops, best of 3: 3.92 s per loop

Solving a system of equations

Solve for $\mathbf{x}$ in $ \mathbf{A} \mathbf{x} = \mathbf{b}$.

  • numpy again seems faster for most cases (except around $100\times 100$)

In [6]:
ab = npr.rand(100)
bb = npr.rand(1000)
cb = npr.rand(5000)

In [124]:
%timeit linalg.solve(a,ab)


10000 loops, best of 3: 174 µs per loop

In [125]:
%timeit linalg.solve(b,bb)


10 loops, best of 3: 22 ms per loop

In [126]:
%timeit linalg.solve(c,cb)


1 loops, best of 3: 2.58 s per loop

In [127]:
%timeit np.linalg.solve(a,ab)


1000 loops, best of 3: 233 µs per loop

In [128]:
%timeit np.linalg.solve(b,bb)


100 loops, best of 3: 12 ms per loop

In [129]:
%timeit np.linalg.solve(c,cb)


1 loops, best of 3: 1.25 s per loop

In [131]:
tinyb = npr.rand(len(tiny))

In [133]:
%timeit linalg.solve(tiny,tinyb)


10000 loops, best of 3: 34.3 µs per loop

In [132]:
%timeit np.linalg.solve(tiny,tinyb)


100000 loops, best of 3: 14.4 µs per loop

Determinants

  • A matrix of coefficients in a system of linear equations has a unique solution exactly when the determinant is nonzero (zero means either no or many solutions)
  • A matrix corresponding to a linear transformation of a vector space has an inverse operation exactly when the determinant is nonzero
    • Absolute value of determinant (when the matrix has real entries) gives the "scale factor" by which area/volume/measure is multiplied under the transformation
    • Sign of the determinant indicates whether the transformation preserves orientation

In [134]:
%timeit linalg.det(c)


1 loops, best of 3: 2.19 s per loop

In [136]:
%timeit np.linalg.det(c)


1 loops, best of 3: 1.16 s per loop

Norm

Least squares solution to $\mathbf{Ax} = \mathbf{b}$


In [7]:
%timeit linalg.lstsq(a,ab)


100 loops, best of 3: 2.3 ms per loop

In [8]:
%timeit linalg.lstsq(b,bb)


1 loops, best of 3: 464 ms per loop

Moore-Penrose pseudo-inverse

  • Of a matrix or Hermitian matrix

Kronecker product

Eigenvalue problems

  • Solve (or compute eigenvalues from) ordinary or generalized eigenvalue problems for square, Hermitian, symmetric, or banded matrices

In [10]:
%timeit linalg.eig(a)


100 loops, best of 3: 8.11 ms per loop

In [11]:
%timeit linalg.eig(b)


1 loops, best of 3: 1.32 s per loop

In [12]:
%timeit np.linalg.eig(b)


1 loops, best of 3: 1.48 s per loop

In [33]:
%timeit linalg.eigvals(b)


1 loops, best of 3: 806 ms per loop

Decompositions

LU decomposition

  • Decompose a matrix into a lower triangular matrix $L$ and an upper triangular matrix $U$ $$ A = LU$$
    • Useful in computing matrix inverses, determinants, and solutions to systems of linear equations

Singular value decomposition


In [19]:
%timeit linalg.svd(a)


100 loops, best of 3: 2.51 ms per loop

In [20]:
%timeit linalg.svd(b)


1 loops, best of 3: 453 ms per loop

Cholesky decomposition

Polar decomposition

QR decomposition

QZ decomposition

  • For generalized eigenvalues of a pair of matrices

Schur decomposition

Interpolative decompositions

  • Cheaper and more efficient to construct than SVD
  • Preserves structure of $\mathbf{A}$

In [14]:
from scipy.linalg import interpolative

In [17]:
%timeit interpolative.svd(a,0.1)


100 loops, best of 3: 6.5 ms per loop

In [23]:
%timeit interpolative.svd(a,5)


1000 loops, best of 3: 340 µs per loop

In [25]:
%timeit interpolative.svd(b,50)


10 loops, best of 3: 40.1 ms per loop

In [28]:
%timeit interpolative.svd(c,50)


1 loops, best of 3: 1.1 s per loop

In [30]:
%timeit interpolative.estimate_spectral_norm(a)


1000 loops, best of 3: 214 µs per loop

In [31]:
%timeit interpolative.estimate_spectral_norm(b)


100 loops, best of 3: 7.04 ms per loop

In [32]:
%timeit interpolative.estimate_spectral_norm(c)


1 loops, best of 3: 315 ms per loop

Matrix

Special matrices