numpy matrix math

! Linear Algebra in Numpy

This is a collection of linear algebra routines in NumPy

I find myself frequently reinventing the wheel, so these are minimal working examples for different problems I have

cholesky

import numpy as np

b=np.random.randn(4,4)
c=np.random.randn(4)
c=abs(c)
d=b.dot(np.diag(c)).dot(b.T)
# must be positive semidefinite and symmetric - numpy may require positive
# definite, but algorithms exist for approximate cholesky orthogonalization
c=np.linalg.cholesky(d)
print d-c.dot(c.T)

eigendecomposition

import numpy as np

b=np.random.randn(5,5)
b=b+b.T
evals,evecs=np.linalg.eig(b)
print b-evecs.dot(np.diag(evals)).dot(evecs.transpose())

outer products

import numpy as np

a=np.arange(1,4,1)
b=np.outer(a,a)
evals,evecs=np.linalg.eig(b)
d=evecs[:,0]
print np.outer(d,d)*evals[0]-b

svd versus eigenvalue decomp

C=np.random.randn(2,2)
#C=.5*(C+C.T) # make symmetric
#right.dot(right.T) #only diagonal for symmetric matrices
evals,left,right=LA.eig(C,left=True)
left/=np.linalg.norm(left,axis=0)
right/=np.linalg.norm(right,axis=0)
print(np.conjugate(right.T).dot(C).dot(right).diagonal())
print((C.dot(right)/(right))[0,:])
print(np.conjugate(left.T).dot(C).dot(left).diagonal())
print((np.conjugate(left.T).dot(C)/np.conjugate(left.T))[:,0])
print(evals)

u,s,v=LA.svd(C)
print(np.linalg.norm(u.dot(np.diag(s)).dot(v)-C))
print(C.dot(v.T.dot(np.diag(s**-1)).dot(u.T)))
print(v.T.dot(np.diag(s**-1)).dot(u.T).dot(C))
print(s)

symmetric orthogonalization

import numpy as np

b=np.random.randn(4,4)
c=np.random.randn(4)
c=abs(c)
d=b.dot(np.diag(c)).dot(b.T)
evals,evecs=np.linalg.eig(d)
s=evecs.dot(np.diag(evals**-.5)).dot(evecs.T)
print s.dot(d).dot(s)-np.identity(4)