[SciPy-User] struggling and fighting with linalg
josef.pktd at gmail.com
josef.pktd at gmail.com
Thu Oct 21 12:03:39 EDT 2010
There are a lot of recommendations in the mailing list and recipes in
various packages how to use linalg more efficiently.
However, yesterday I wasted several hours how to do the multivariate
normal distribution with cholesky factor and solve.
The other problem we have in some of the code is that svd or similar
are calculated several times when individual numpy/scipy functions are
called (e.g. pinv, logdet). Dag has oomatrix, and there are some
packages on the matlab fileexchange (I didn't look at those yet).
As a collection of recipes and reminder for myself, which is which and
what needs to be transposed in s v d, I started to write some linalg
array classes (in attachment). It's essentially just a holder for the
decomposition and operations defined with it, using arrays not
matrices and without any operator overloading.
Still incomplete and only partially tested, and I haven't yet gone
back to doing some statistics with it.
Initially, I was looking for a cholesky decomposition of
(block-)Toeplitz matrices. But those and sparse cholesky seem to be
unavailable in the BSD-compatible landscape, but are for GPL.
If anyone can provide code snippets or (short) references for missing
parts for the linalg, it would reduce my trial-and-error time. I
haven't worked my way through all the suggestions on the mailing list
yet. As for my initial intention, I would also like to get a
multivariate normal class with batteries included, so I don't have to
figure it out and test each time I need a piece.
Josef
-------------- next part --------------
'''Recipes for more efficient work with linalg using classes
intended for use for multivariate normal and linear regression
calculations
x is the data (nobs, nvars)
m is the moment matrix (x'x) or a covariance matrix Sigma
examples:
x'sigma^{-1}x
z = Px where P=Sigma^{-1/2} or P=Sigma^{1/2}
Initially assume positive definite, then add spectral cutoff and
regularization of moment matrix, and extend to PCA
Author: josef-pktd
Created on 2010-10-20
'''
import numpy as np
from scipy import linalg
#this has been copied from nitime a long time ago
#TODO: ceck whether class has changed in nitime
class OneTimeProperty(object):
"""A descriptor to make special properties that become normal attributes.
This is meant to be used mostly by the auto_attr decorator in this module.
Author: Fernando Perez, copied from nitime
"""
def __init__(self,func):
"""Create a OneTimeProperty instance.
Parameters
----------
func : method
The method that will be called the first time to compute a value.
Afterwards, the method's name will be a standard attribute holding
the value of this computation.
"""
self.getter = func
self.name = func.func_name
def __get__(self,obj,type=None):
"""This will be called on attribute access on the class or instance. """
if obj is None:
# Being called on the class, return the original function. This way,
# introspection works on the class.
#return func
print 'class access'
return self.getter
val = self.getter(obj)
#print "** auto_attr - loading '%s'" % self.name # dbg
setattr(obj, self.name, val)
return val
class PlainMatrixArray(object):
'''Class that defines linalg operation on an array
simplest version as benchmark
linear algebra recipes for multivariate normal and linear
regression calculations
'''
def __init__(self, data=None, sym=None):
if not data is None:
if sym is None:
self.x = np.asarray(data)
self.m = np.dot(self.x.T, self.x)
else:
raise ValueError('data and sym cannot be both given')
elif not sym is None:
self.m = np.asarray(sym)
self.x = np.eye(*self.m.shape) #default
else:
raise ValueError('either data or sym need to be given')
@OneTimeProperty
def minv(self):
return np.linalg.inv(self.m)
@OneTimeProperty
def m_y(self, y):
return np.dot(self.m, y)
def minv_y(self, y):
return np.dot(self.minv, y)
@OneTimeProperty
def mpinv(self):
return linalg.pinv(self.m)
@OneTimeProperty
def xpinv(self):
return linalg.pinv(self.x)
def yt_m_y(self, y):
return np.dot(y.T, np.dot(self.m, y))
def yt_minv_y(self, y):
return np.dot(y.T, np.dot(self.minv, y))
#next two are redundant
def y_m_yt(self, y):
return np.dot(y, np.dot(self.m, y.T))
def y_minv_yt(self, y):
return np.dot(y, np.dot(self.minv, y.T))
@OneTimeProperty
def mdet(self):
return linalg.det(self.m)
@OneTimeProperty
def mlogdet(self):
return np.log(linalg.det(self.m))
@OneTimeProperty
def meigh(self):
evals, evecs = linalg.eigh(self.m)
sortind = np.argsort(evals)[::-1]
return evals[sortind], evecs[:,sortind]
@OneTimeProperty
def mhalf(self):
evals, evecs = self.meigh
return np.dot(np.diag(evals**0.5), evecs.T)
#return np.dot(evecs, np.dot(np.diag(evals**0.5), evecs.T))
#return np.dot(evecs, 1./np.sqrt(evals) * evecs.T))
@OneTimeProperty
def minvhalf(self):
evals, evecs = self.meigh
return np.dot(evecs, 1./np.sqrt(evals) * evecs.T)
class SvdArray(PlainMatrixArray):
'''Class that defines linalg operation on an array
svd version, where svd is taken on original data array, if
or when it matters
no spectral cutoff in first version
'''
def __init__(self, data=None, sym=None):
super(SvdArray, self).__init__(data=data, sym=sym)
u, s, v = np.linalg.svd(self.x, full_matrices=1)
self.u, self.s, self.v = u, s, v
self.sdiag = linalg.diagsvd(s, *x.shape)
self.sinvdiag = linalg.diagsvd(1./s, *x.shape)
def _sdiagpow(self, p):
return linalg.diagsvd(np.power(self.s, p), *x.shape)
@OneTimeProperty
def minv(self):
sinvv = np.dot(self.sinvdiag, self.v)
return np.dot(sinvv.T, sinvv)
@OneTimeProperty
def meigh(self):
evecs = self.v.T
evals = self.s**2
return evals, evecs
@OneTimeProperty
def mdet(self):
return self.meigh[0].prod()
@OneTimeProperty
def mlogdet(self):
return np.log(self.meigh[0]).sum()
@OneTimeProperty
def mhalf(self):
return np.dot(np.diag(self.s), self.v)
@OneTimeProperty
def xxthalf(self):
return np.dot(self.u, self.sdiag)
@OneTimeProperty
def xxtinvhalf(self):
return np.dot(self.u, self.sinvdiag)
class CholArray(PlainMatrixArray):
'''Class that defines linalg operation on an array
cholesky version, where svd is taken on original data array, if
or when it matters
plan: use cholesky factor and cholesky solve
nothing implemented yet
'''
def __init__(self, data=None, sym=None):
super(SvdArray, self).__init__(data=data, sym=sym)
def testcompare(m1, m2):
from numpy.testing import assert_almost_equal, assert_approx_equal
decimal = 12
#inv
assert_almost_equal(m1.minv, m2.minv, decimal=decimal)
#matrix half and invhalf
#fix sign in test, should this be standardized
s1 = np.sign(m1.mhalf.sum(1))[:,None]
s2 = np.sign(m2.mhalf.sum(1))[:,None]
scorr = s1/s2
assert_almost_equal(m1.mhalf, m2.mhalf * scorr, decimal=decimal)
assert_almost_equal(m1.minvhalf, m2.minvhalf, decimal=decimal)
#eigenvalues, eigenvectors
evals1, evecs1 = m1.meigh
evals2, evecs2 = m2.meigh
assert_almost_equal(evals1, evals2, decimal=decimal)
#normalization can be different: evecs in columns
s1 = np.sign(evecs1.sum(0))
s2 = np.sign(evecs2.sum(0))
scorr = s1/s2
assert_almost_equal(evecs1, evecs2 * scorr, decimal=decimal)
#determinant
assert_approx_equal(m1.mdet, m2.mdet, significant=13)
assert_approx_equal(m1.mlogdet, m2.mlogdet, significant=13)
####### helper function for interactive work
def tiny2zero(x, eps = 1e-15):
'''replace abs values smaller than eps by zero, makes copy
'''
mask = np.abs(x.copy()) < eps
x[mask] = 0
return x
def maxabs(x):
return np.max(np.abs(x))
#if __name__ == '__main__':
n = 5
y = np.arange(n)
x = np.random.randn(100,n)
autocov = 2*0.8**np.arange(n) +0.01 * np.random.randn(n)
sigma = linalg.toeplitz(autocov)
mat = PlainMatrixArray(sym=sigma)
print tiny2zero(mat.mhalf)
mih = mat.minvhalf
print tiny2zero(mih) #for nicer printing
mat2 = PlainMatrixArray(data=x)
print maxabs(mat2.yt_minv_y(np.dot(x.T, x)) - mat2.m)
print tiny2zero(mat2.minv_y(mat2.m))
mat3 = SvdArray(data=x)
print mat3.meigh[0]
print mat2.meigh[0]
testcompare(mat2, mat3)
'''
m = np.dot(x.T, x)
u,s,v = np.linalg.svd(x, full_matrices=1)
Sig = linalg.diagsvd(s,*x.shape)
>>> np.max(np.abs(np.dot(u, np.dot(Sig, v)) - x))
3.1086244689504383e-015
>>> np.max(np.abs(np.dot(u.T, u) - np.eye(100)))
3.3306690738754696e-016
>>> np.max(np.abs(np.dot(v.T, v) - np.eye(5)))
6.6613381477509392e-016
>>> np.max(np.abs(np.dot(Sig.T, Sig) - np.diag(s**2)))
5.6843418860808015e-014
>>> evals,evecs = linalg.eigh(np.dot(x.T, x))
>>> evals[::-1]
array([ 123.36404464, 112.17036442, 102.04198468, 76.60832278,
74.70484487])
>>> s**2
array([ 123.36404464, 112.17036442, 102.04198468, 76.60832278,
74.70484487])
>>> np.max(np.abs(np.dot(v.T, np.dot(np.diag(s**2), v)) - m))
1.1368683772161603e-013
>>> us = np.dot(u, Sig)
>>> np.max(np.abs(np.dot(us, us.T) - np.dot(x, x.T)))
1.0658141036401503e-014
>>> sv = np.dot(Sig, v)
>>> np.max(np.abs(np.dot(sv.T, sv) - np.dot(x.T, x)))
1.1368683772161603e-013
'''
More information about the SciPy-User
mailing list