771 lines
21 KiB
Python
771 lines
21 KiB
Python
"""Module contain algorithms for low-rank L-shaped model.
|
|
|
|
"""
|
|
|
|
__all__ = ['pca', 'pcr', 'pls', 'nipals_lpls']
|
|
__docformat__ = "restructuredtext en"
|
|
|
|
from math import sqrt as msqrt
|
|
|
|
from numpy import dot,empty,zeros,apply_along_axis,newaxis,finfo,sqrt,r_,expand_dims,\
|
|
minimum
|
|
from numpy.linalg import inv, svd
|
|
from scipy.sandbox import arpack
|
|
|
|
def pca(X, aopt, scale='scores', mode='normal', center_axis=0):
|
|
""" Principal Component Analysis.
|
|
|
|
PCA is a low rank bilinear aprroximation to a data matrix that sequentially
|
|
extracts orthogonal components of maximum variance.
|
|
|
|
Parameters:
|
|
|
|
X : {array}
|
|
Data measurement matrix, (samples x variables)
|
|
aopt : {integer}
|
|
Number of components to use, aopt<=min(samples, variables)
|
|
center_axis : {integer}
|
|
Center along given axis. If neg.: no centering (-inf,..., matrix modes)
|
|
|
|
Returns:
|
|
|
|
T : {array}
|
|
Scores, (samples, components)
|
|
P : {array}
|
|
Loadings, (variables, components)
|
|
E : {array}
|
|
Residuals, (samples, variables)
|
|
evx : {array}
|
|
X-explained variance, (components,)
|
|
mnx : {array}
|
|
X location, (variables,)
|
|
aopt : {integer}
|
|
The number of components used
|
|
ssqx : {list}
|
|
Sum of squared residuals in X along each dimesion
|
|
[(samples, ), (variables,)]
|
|
leverage : {array}
|
|
Leverages, (samples,)
|
|
|
|
OtherParameters:
|
|
|
|
scale : {string}, optional
|
|
Where to put the weights [['scores'], 'loadings']
|
|
mode : {string}, optional
|
|
Amount of info retained, [['normal'], 'fast', 'detailed']
|
|
|
|
|
|
:SeeAlso:
|
|
|
|
`center` : Data centering
|
|
|
|
|
|
*Notes*
|
|
|
|
Uses kernel speed-up if m>>n or m<<n.
|
|
|
|
If residuals turn rank deficient, a lower number of component than given
|
|
in input will be used.
|
|
|
|
*Examples*:
|
|
|
|
>>> import scipy,engines
|
|
>>> a=scipy.asarray([[1,2,3],[2,4,5]])
|
|
>>> dat=engines.pca(a, 2)
|
|
>>> dat['evx']
|
|
array([0.,99.8561562, 100.])
|
|
|
|
"""
|
|
|
|
m, n = X.shape
|
|
assert(aopt<=min(m,n))
|
|
if center_axis>=0:
|
|
X = X - expand_dims(X.mean(center_axis), center_axis)
|
|
if m>(n+100) or n>(m+100):
|
|
u, s, v = esvd(X, aopt)
|
|
else:
|
|
u, s, vt = svd(X, 0)
|
|
v = vt.T
|
|
u = u[:,:aopt]
|
|
s = s[:aopt]
|
|
v = v[:,:aopt]
|
|
|
|
# ranktest
|
|
tol = 1e-10
|
|
eff_rank = sum(s>s[0]*tol)
|
|
aopt = minimum(aopt, eff_rank)
|
|
T = u*s
|
|
s = s[:aopt]
|
|
T = T[:,:aopt]
|
|
P = v[:,:aopt]
|
|
e = s**2
|
|
|
|
if scale=='loads':
|
|
T = T/s
|
|
P = P*s
|
|
|
|
if mode == 'fast':
|
|
return {'T':T, 'P':P, 'aopt':aopt}
|
|
|
|
if mode=='detailed':
|
|
E = empty((aopt, m, n))
|
|
ssq = []
|
|
lev = []
|
|
for ai in range(aopt):
|
|
E[ai,:,:] = X - dot(T[:,:ai+1], P[:,:ai+1].T)
|
|
ssq.append([(E[ai,:,:]**2).mean(0), (E[ai,:,:]**2).mean(1)])
|
|
if scale=='loads':
|
|
lev.append([((s*T)**2).sum(1), (P**2).sum(1)])
|
|
else:
|
|
lev.append([(T**2).sum(1), ((s*P)**2).sum(1)])
|
|
else:
|
|
# residuals
|
|
E = X - dot(T, P.T)
|
|
sep = E**2
|
|
ssq = [sep.sum(0), sep.sum(1)]
|
|
# leverages
|
|
if scale=='loads':
|
|
lev = [(1./m)+(T**2).sum(1), (1./n)+((P/s)**2).sum(1)]
|
|
else:
|
|
lev = [(1./m)+((T/s)**2).sum(1), (1./n)+(P**2).sum(1)]
|
|
# variances
|
|
expvarx = r_[0, 100*e.cumsum()/(X*X).sum()]
|
|
|
|
return {'T': T, 'P': P, 'E': E, 'evx': expvarx, 'leverage': lev, 'ssqx': ssq,
|
|
'aopt': aopt, 'eigvals': e}
|
|
|
|
def pcr(a, b, aopt, scale='scores',mode='normal',center_axis=0):
|
|
""" Principal Component Regression.
|
|
|
|
Performs PCR on given matrix and returns results in a dictionary.
|
|
|
|
Parameters:
|
|
|
|
a : array
|
|
Data measurement matrix, (samples x variables)
|
|
b : array
|
|
Data response matrix, (samples x responses)
|
|
aopt : int
|
|
Number of components to use, aopt<=min(samples, variables)
|
|
|
|
Returns:
|
|
|
|
results : dict
|
|
keys -- values, T -- scores, P -- loadings, E -- residuals,
|
|
levx -- leverages, ssqx -- sum of squares, expvarx -- cumulative
|
|
explained variance, aopt -- number of components used
|
|
|
|
OtherParameters:
|
|
|
|
mode : str
|
|
Amount of info retained, ('fast', 'normal', 'detailed')
|
|
center_axis : int
|
|
Center along given axis. If neg.: no centering (-inf,..., matrix modes)
|
|
|
|
SeeAlso:
|
|
|
|
- pca : other blm
|
|
- pls : other blm
|
|
- lpls : other blm
|
|
|
|
*Notes*
|
|
|
|
-----
|
|
|
|
Uses kernel speed-up if m>>n or m<<n.
|
|
|
|
If residuals turn rank deficient, a lower number of component than given
|
|
in input will be used. The number of components used is given in results-dict.
|
|
|
|
|
|
Examples
|
|
--------
|
|
|
|
>>> import scipy,engines
|
|
>>> a=scipy.asarray([[1,2,3],[2,4,5]])
|
|
>>> b=scipy.asarray([[1,1],[2,3]])
|
|
>>> dat=engines.pcr(a, 2)
|
|
>>> dat['evx']
|
|
array([0.,99.8561562, 100.])
|
|
|
|
"""
|
|
try:
|
|
k, l = b.shape
|
|
except:
|
|
b = atleast_2d(b).T
|
|
k, l = b.shape
|
|
if center_axis>=0:
|
|
b = b - expand_dims(b.mean(center_axis), center_axis)
|
|
dat = pca(a, aopt=aopt, scale=scale, mode=mode, center_axis=center_axis)
|
|
T = dat['T']
|
|
weights = apply_along_axis(vnorm, 0, T)**2
|
|
if scale=='loads':
|
|
Q = dot(b.T, T*weights)
|
|
else:
|
|
Q = dot(b.T, T/weights)
|
|
|
|
if mode=='fast':
|
|
dat.update({'Q':Q})
|
|
return dat
|
|
|
|
if mode=='detailed':
|
|
F = empty((aopt, k, l))
|
|
ssqy = []
|
|
for i in range(aopt):
|
|
F[i,:,:] = b - dot(T[:,:i+1], Q[:,:i+1].T)
|
|
ssqy.append([(F[i,:,:]**2).mean(0), (F[i,:,:]**2).mean(1)])
|
|
else:
|
|
F = b - dot(T, Q.T)
|
|
sepy = F**2
|
|
ssqy = [sepy.sum(0), sepy.sum(1)]
|
|
|
|
expvary = r_[0, 100*((T**2).sum(0)*(Q**2).sum(0)/(b**2).sum()).cumsum()[:aopt]]
|
|
|
|
dat.update({'Q': Q, 'F': F, 'evy': expvary, 'ssqy': ssqy})
|
|
return dat
|
|
|
|
def pls(X, Y, aopt=2, scale='scores', mode='normal', center_axis=-1):
|
|
"""Partial Least Squares Regression.
|
|
|
|
Performs PLS on given matrix and returns results in a dictionary.
|
|
|
|
*Parameters*:
|
|
|
|
X : {array}
|
|
Data measurement matrix, (samples x variables)
|
|
Y : {array}
|
|
Data response matrix, (samples x responses)
|
|
aopt : {integer}, optional
|
|
Number of components to use, aopt<=min(samples, variables)
|
|
scale : ['scores', 'loadings'], optional
|
|
Which component should get the scale
|
|
center_axis : {-1, integer}
|
|
Perform centering across given axis, (-1 is no centering)
|
|
|
|
*Returns*:
|
|
|
|
T : {array}
|
|
X-scores
|
|
W : {array}
|
|
X-loading-weights
|
|
P : {array}
|
|
X-loadings
|
|
R : {array}
|
|
X-loadings-basis
|
|
Q : {array}
|
|
Y-loadings
|
|
B : {array}
|
|
Regression coefficients
|
|
E : {array}
|
|
X-block residuals
|
|
F : {array}
|
|
Y-block residuals
|
|
evx : {array}
|
|
X-explained variance
|
|
evy : {array}
|
|
Y-explained variance
|
|
mnx : {array}
|
|
X location
|
|
mny : {array}
|
|
Y location
|
|
aopt : {array}
|
|
The number of components used
|
|
ssqx : {list}, optional
|
|
Sum of squared residuals in X along each dimesion
|
|
ssqy : {list}
|
|
Sum of squared residuals in Y along each dimesion
|
|
leverage : {array}
|
|
Sample leverages
|
|
|
|
*OtherParameters*:
|
|
|
|
mode : ['normal', 'fast', 'detailed'], optional
|
|
How much details to compute
|
|
|
|
*SeeAlso*:
|
|
|
|
`center` : data centering
|
|
|
|
*Notes*
|
|
|
|
- The output with mode='fast' will only return T and W
|
|
|
|
- If residuals turn rank deficient, a lower number of component than given in input will be used. The number of components used is given in results.
|
|
|
|
|
|
*Examples*
|
|
|
|
>>> import numpy, engines
|
|
>>> a = numpy.asarray([[1,2,3],[2,4,5]])
|
|
>>> b = numpy.asarray([[1,1],[2,3]])
|
|
>>> dat =engines.pls(a, b, 2)
|
|
>>> dat['evx']
|
|
array([0.,99.8561562, 100.])
|
|
|
|
"""
|
|
|
|
m, n = X.shape
|
|
try:
|
|
k, l = Y.shape
|
|
except:
|
|
Y = atleast_2d(Y).T
|
|
k, l = Y.shape
|
|
assert(m==k)
|
|
assert(aopt<min(m, n))
|
|
mnx, mny = 0,0
|
|
if center_axis>=0:
|
|
mnx = expand_dims(X.mean(center_axis), center_axis)
|
|
X = X - mnx
|
|
mny = expand_dims(Y.mean(center_axis), center_axis)
|
|
Y = Y - mny
|
|
|
|
W = empty((n, aopt))
|
|
P = empty((n, aopt))
|
|
R = empty((n, aopt))
|
|
Q = empty((l, aopt))
|
|
T = empty((m, aopt))
|
|
B = empty((aopt, n, l))
|
|
tt = empty((aopt,))
|
|
|
|
XY = dot(X.T, Y)
|
|
for i in range(aopt):
|
|
if XY.shape[1]==1: #pls 1
|
|
w = XY.reshape(n, l)
|
|
w = w/vnorm(w)
|
|
elif n<l: # more yvars than xvars
|
|
s, w = arpack.eigen_symmetric(dot(XY, XY.T),k=1, tol=1e-10, maxiter=100)
|
|
#w, s, vh = svd(dot(XY, XY.T))
|
|
#w = w[:,:1]
|
|
else: # more xvars than yvars
|
|
s, q = arpack.eigen_symmetric(dot(XY.T, XY), k=1, tol=1e-10, maxiter=100)
|
|
#q, s, vh = svd(dot(XY.T, XY))
|
|
#q = q[:,:1]
|
|
|
|
w = dot(XY, q)
|
|
w = w/vnorm(w)
|
|
r = w.copy()
|
|
if i>0:
|
|
for j in range(0, i, 1):
|
|
r = r - dot(P[:,j].T, w)*R[:,j][:,newaxis]
|
|
|
|
t = dot(X, r)
|
|
tt[i] = tti = dot(t.T, t).ravel()
|
|
p = dot(X.T, t)/tti
|
|
q = dot(r.T, XY).T/tti
|
|
XY = XY - dot(p, q.T)*tti
|
|
T[:,i] = t.ravel()
|
|
W[:,i] = w.ravel()
|
|
|
|
if mode=='fast' and i==aopt-1:
|
|
if scale=='loads':
|
|
tnorm = sqrt(tt)
|
|
T = T/tnorm
|
|
W = W*tnorm
|
|
return {'T':T, 'W':W}
|
|
|
|
P[:,i] = p.ravel()
|
|
R[:,i] = r.ravel()
|
|
Q[:,i] = q.ravel()
|
|
B[i] = dot(R[:,:i+1], Q[:,:i+1].T)
|
|
|
|
qnorm = apply_along_axis(vnorm, 0, Q)
|
|
tnorm = sqrt(tt)
|
|
pp = (P**2).sum(0)
|
|
if mode=='detailed':
|
|
E = empty((aopt, m, n))
|
|
F = empty((aopt, k, l))
|
|
ssqx, ssqy = [], []
|
|
leverage = empty((aopt, m))
|
|
#h2x = [] #hotellings T^2
|
|
#h2y = []
|
|
for ai in range(aopt):
|
|
E[ai,:,:] = X - dot(T[:,:ai+1], P[:,:ai+1].T)
|
|
F[ai,:,:] = Y - dot(T[:,:i], Q[:,:i].T)
|
|
ssqx.append([(E[ai,:,:]**2).mean(0), (E[ai,:,:]**2).mean(1)])
|
|
ssqy.append([(F[ai,:,:]**2).mean(0), (F[ai,:,:]**2).mean(1)])
|
|
leverage[ai,:] = 1./m + ((T[:,:ai+1]/tnorm[:ai+1])**2).sum(1)
|
|
#h2y.append(1./k + ((Q[:,:ai+1]/qnorm[:ai+1])**2).sum(1))
|
|
else:
|
|
# residuals
|
|
E = X - dot(T, P.T)
|
|
F = Y - dot(T, Q.T)
|
|
sepx = E**2
|
|
ssqx = [sepx.sum(0), sepx.sum(1)]
|
|
sepy = F**2
|
|
ssqy = [sepy.sum(0), sepy.sum(1)]
|
|
leverage = 1./m + ((T/tnorm)**2).sum(1)
|
|
|
|
# variances
|
|
tp= tt*pp
|
|
tq = tt*qnorm*qnorm
|
|
expvarx = r_[0, 100*tp/(X*X).sum()]
|
|
expvary = r_[0, 100*tq/(Y*Y).sum()]
|
|
|
|
if scale=='loads':
|
|
T = T/tnorm
|
|
W = W*tnorm
|
|
Q = Q*tnorm
|
|
P = P*tnorm
|
|
|
|
return {'Q': Q, 'P': P, 'T': T, 'W': W, 'R': R, 'E': E, 'F': F, 'B': B,
|
|
'evx': expvarx, 'evy': expvary, 'ssqx': ssqx, 'ssqy': ssqy,
|
|
'leverage': leverage, 'mnx': mnx, 'mny': mny}
|
|
|
|
def nipals_lpls(X, Y, Z, a_max, alpha=.7, mean_ctr=[2, 0, 2], scale='scores', zorth = False, verbose=False):
|
|
""" L-shaped Partial Least Sqaures Regression by the nipals algorithm.
|
|
|
|
An L-shaped low rank model aproximates three matrices in a hyploid
|
|
structure. That means that the main matrix (X), has one matrix asociated
|
|
with its row space and one to its column space. A simultanously low rank estiamte
|
|
of these three matrices tries to discover common directions/subspaces.
|
|
|
|
*Parameters*:
|
|
|
|
X : {array}
|
|
Main data matrix (m, n)
|
|
Y : {array}
|
|
External row data (m, l)
|
|
Z : {array}
|
|
External column data (n, o)
|
|
a_max : {integer}
|
|
Maximum number of components to calculate (0, min(m,n))
|
|
alpha : {float}, optional
|
|
Parameter to control the amount of influence from Z-matrix.
|
|
0 is none, which returns a pls-solution, 1 is max
|
|
mean_center : {array-like}, optional
|
|
A three element array-like structure with elements in [-1,0,1,2],
|
|
that decides the type of centering used.
|
|
-1 : nothing
|
|
0 : row center
|
|
1 : column center
|
|
2 : double center
|
|
|
|
*Returns*:
|
|
|
|
T : {array}
|
|
X-scores
|
|
W : {array}
|
|
X-weights/Z-weights
|
|
P : {array}
|
|
X-loadings
|
|
Q : {array}
|
|
Y-loadings
|
|
U : {array}
|
|
X-Y relation
|
|
L : {array}
|
|
Z-scores
|
|
K : {array}
|
|
Z-loads
|
|
B : {array}
|
|
Regression coefficients X->Y
|
|
evx : {array}
|
|
X-explained variance
|
|
evy : {array}
|
|
Y-explained variance
|
|
evz : {array}
|
|
Z-explained variance
|
|
mnx : {array}
|
|
X location
|
|
mny : {array}
|
|
Y location
|
|
mnz : {array}
|
|
Z location
|
|
|
|
*OtherParameters*:
|
|
|
|
scale : {'scores', 'loads'}, optional
|
|
Option to decide on where the scale goes.
|
|
zorth : {False, boolean}, optional
|
|
Option to force orthogonality between latent components
|
|
in Z
|
|
verbose : {boolean}, optional
|
|
Verbosity of console output. For use in debugging.
|
|
|
|
*References*
|
|
|
|
Saeboe et al., LPLS-regression: a method for improved prediction and
|
|
classification through inclusion of background information on
|
|
predictor variables, J. of chemometrics and intell. laboratory syst.
|
|
|
|
Martens et.al, Regression of a data matrix on descriptors of
|
|
both its rows and of its columns via latent variables: L-PLSR,
|
|
Computational statistics & data analysis, 2005
|
|
|
|
"""
|
|
m, n = X.shape
|
|
k, l = Y.shape
|
|
u, o = Z.shape
|
|
max_rank = min(m, n)
|
|
assert (a_max>0 and a_max<max_rank), "Number of comp error:\
|
|
tried:%d, max_rank:%d" %(a_max,max_rank)
|
|
|
|
if mean_ctr!=None:
|
|
xctr, yctr, zctr = mean_ctr
|
|
X, mnX = center(X, xctr)
|
|
Y, mnY = center(Y, yctr)
|
|
Z, mnZ = center(Z, zctr)
|
|
|
|
varX = (X**2).sum()
|
|
varY = (Y**2).sum()
|
|
varZ = (Z**2).sum()
|
|
|
|
# initialize
|
|
U = empty((k, a_max))
|
|
Q = empty((l, a_max))
|
|
T = zeros((m, a_max))
|
|
W = empty((n, a_max))
|
|
P = empty((n, a_max))
|
|
K = empty((o, a_max))
|
|
L = empty((u, a_max))
|
|
B = empty((a_max, n, l))
|
|
E = X.copy()
|
|
F = Y.copy()
|
|
G = Z.copy()
|
|
#b0 = empty((a_max, 1, l))
|
|
var_x = empty((a_max,))
|
|
var_y = empty((a_max,))
|
|
var_z = empty((a_max,))
|
|
|
|
MAX_ITER = 4500
|
|
LIM = finfo(X.dtype).resolution
|
|
is_rd = False
|
|
for a in range(a_max):
|
|
if verbose:
|
|
print "\nWorking on comp. %s" %a
|
|
u = F[:,:1]
|
|
w = E[:1,:].T
|
|
l = G[:,:1]
|
|
diff = 1
|
|
niter = 0
|
|
while (diff>LIM and niter<MAX_ITER):
|
|
niter += 1
|
|
u1 = u.copy()
|
|
w1 = w.copy()
|
|
l1 = l.copy()
|
|
w = dot(E.T, u)
|
|
wn = msqrt(dot(w.T, w))
|
|
if wn < LIM:
|
|
print "Rank exhausted in X! Comp: %d " %a
|
|
is_rd = True
|
|
break
|
|
w = w/wn
|
|
#w = w/dot(w.T, w)
|
|
l = dot(G, w)
|
|
k = dot(G.T, l)
|
|
k = k/msqrt(dot(k.T, k))
|
|
#k = k/dot(k.T, k)
|
|
w = alpha*k + (1-alpha)*w
|
|
w = w/msqrt(dot(w.T, w))
|
|
t = dot(E, w)
|
|
c = dot(F.T, t)
|
|
c = c/msqrt(dot(c.T, c))
|
|
u = dot(F, c)
|
|
diff = dot((u - u1).T, (u - u1))
|
|
if verbose:
|
|
if niter==MAX_ITER:
|
|
print "Maximum nunber of iterations reached!"
|
|
print "Iterations: %d " %niter
|
|
print "Error: %.2E" %diff
|
|
|
|
if is_rd:
|
|
print "Hei og haa ... rank deficient, this should really not happen"
|
|
break
|
|
|
|
tt = dot(t.T, t)
|
|
p = dot(E.T, t)/tt
|
|
q = dot(F.T, t)/tt
|
|
if zorth:
|
|
k = dot(G.T, l)/dot(l.T, l)
|
|
else:
|
|
k = w
|
|
l = dot(G, w)
|
|
|
|
U[:,a] = u.ravel()
|
|
W[:,a] = w.ravel()
|
|
P[:,a] = p.ravel()
|
|
T[:,a] = t.ravel()
|
|
Q[:,a] = q.ravel()
|
|
L[:,a] = l.ravel()
|
|
K[:,a] = k.ravel()
|
|
|
|
# rank-one deflations
|
|
E = E - dot(t, p.T)
|
|
F = F - dot(t, q.T)
|
|
G = G - dot(l, k.T)
|
|
|
|
var_x[a] = pow(E, 2).sum()
|
|
var_y[a] = pow(F, 2).sum()
|
|
var_z[a] = pow(G, 2).sum()
|
|
|
|
B[a] = dot(dot(W[:,:a+1], inv(dot(P[:,:a+1].T, W[:,:a+1]))), Q[:,:a+1].T)
|
|
#b0[a] = mnY - dot(mnX, B[a])
|
|
|
|
|
|
# variance explained
|
|
evx = 100.*(1 - var_x/varX)
|
|
evy = 100.*(1 - var_y/varY)
|
|
evz = 100.*(1 - var_z/varZ)
|
|
|
|
if scale=='loads':
|
|
tnorm = apply_along_axis(vnorm, 0, T)
|
|
T = T/tnorm
|
|
W = W*tnorm
|
|
Q = Q*tnorm
|
|
knorm = apply_along_axis(vnorm, 0, K)
|
|
L = L*knorm
|
|
K = K/knorm
|
|
|
|
return {'T':T, 'W':W, 'P':P, 'Q':Q, 'U':U, 'L':L, 'K':K, 'B':B, 'E': E, 'F': F, 'G': G, 'evx':evx, 'evy':evy, 'evz':evz,'mnx': mnX, 'mny': mnY, 'mnz': mnZ}
|
|
|
|
def vnorm(a):
|
|
"""Returns the norm of a vector.
|
|
|
|
*Parameters*:
|
|
|
|
a : {array}
|
|
Input data, 1-dim, or column vector (m, 1)
|
|
|
|
*Returns*:
|
|
|
|
a_norm : {array}
|
|
Norm of input vector
|
|
|
|
"""
|
|
return msqrt(dot(a.T,a))
|
|
|
|
def center(a, axis):
|
|
""" Matrix centering.
|
|
|
|
*Parameters*:
|
|
|
|
a : {array}
|
|
Input data
|
|
axis : {integer}
|
|
Which centering to perform.
|
|
0 = col center, 1 = row center, 2 = double center
|
|
-1 = nothing
|
|
|
|
*Returns*:
|
|
|
|
a_centered : {array}
|
|
Centered data matrix
|
|
mn : {array}
|
|
Location vector/matrix
|
|
"""
|
|
|
|
# check if we have a vector
|
|
is_vec = len(a.shape)==1
|
|
if not is_vec:
|
|
is_vec = a.shape[0]==1 or a.shape[1]==1
|
|
if is_vec:
|
|
if axis==2:
|
|
warnings.warn("Double centering of vecor ignored, using ordinary centering")
|
|
if axis==-1:
|
|
mn = 0
|
|
else:
|
|
mn = a.mean()
|
|
return a - mn, mn
|
|
# !!!fixme: use broadcasting
|
|
if axis==-1:
|
|
mn = zeros((1,a.shape[1],))
|
|
#mn = tile(mn, (a.shape[0], 1))
|
|
elif axis==0:
|
|
mn = a.mean(0)[newaxis]
|
|
#mn = tile(mn, (a.shape[0], 1))
|
|
elif axis==1:
|
|
mn = a.mean(1)[:,newaxis]
|
|
#mn = tile(mn, (1, a.shape[1]))
|
|
elif axis==2:
|
|
mn = a.mean(0)[newaxis] + a.mean(1)[:,newaxis] - a.mean()
|
|
return a - mn , a.mean(0)[newaxis]
|
|
else:
|
|
raise IOError("input error: axis must be in [-1,0,1,2]")
|
|
|
|
return a - mn, mn
|
|
|
|
def _scale(a, axis):
|
|
""" Matrix scaling to unit variance.
|
|
|
|
*Parameters*:
|
|
|
|
a : {array}
|
|
Input data
|
|
axis : {integer}
|
|
Which scaling to perform.
|
|
0 = column, 1 = row, -1 = nothing
|
|
|
|
*Returns*:
|
|
|
|
a_scaled : {array}
|
|
Scaled data matrix
|
|
mn : {array}
|
|
Scaling vector/matrix
|
|
"""
|
|
|
|
if axis==-1:
|
|
sc = zeros((a.shape[1],))
|
|
elif axis==0:
|
|
sc = a.std(0)
|
|
elif axis==1:
|
|
sc = a.std(1)[:,newaxis]
|
|
else:
|
|
raise IOError("input error: axis must be in [-1,0,1]")
|
|
|
|
return a - sc, sc
|
|
|
|
def esvd(data, a_max=None):
|
|
""" SVD with kernel calculation
|
|
|
|
Calculate subspaces of X'X or XX' depending on the shape
|
|
of the matrix.
|
|
|
|
Parameters:
|
|
|
|
data : {array}
|
|
Data matrix
|
|
a_max : {integer}
|
|
Number of components to extract
|
|
|
|
Returns:
|
|
|
|
u : {array}
|
|
Right hand eigenvectors
|
|
s : {array}
|
|
Singular values
|
|
v : {array}
|
|
Left hand eigenvectors
|
|
|
|
notes:
|
|
|
|
Uses Anoldi iterations (ARPACK)
|
|
|
|
"""
|
|
|
|
m, n = data.shape
|
|
if m>=n:
|
|
kernel = dot(data.T, data)
|
|
|
|
if a_max==None:
|
|
a_max = n - 1
|
|
s, v = arpack.eigen_symmetric(kernel,k=a_max, which='LM',
|
|
maxiter=200,tol=1e-5)
|
|
s = s[::-1]
|
|
v = v[:,::-1]
|
|
#u, s, vt = svd(kernel)
|
|
#v = vt.T
|
|
s = sqrt(s)
|
|
u = dot(data, v)/s
|
|
else:
|
|
kernel = dot(data, data.T)
|
|
if a_max==None:
|
|
a_max = m -1
|
|
s, u = arpack.eigen_symmetric(kernel,k=a_max, which='LM',
|
|
maxiter=200,tol=1e-5)
|
|
s = s[::-1]
|
|
u = u[:,::-1]
|
|
#u, s, vt = svd(kernel)
|
|
s = sqrt(s)
|
|
v = dot(data.T, u)/s
|
|
|
|
return u, s, v
|