Hauptkomponentenanalyse in Python
https://stackoverflow.com/questions/1730600
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
RussianFrage
Ich möchte hauptsächlich verwenden Komponentenanalyse (PCA) zur Dimensionsreduktion. Hat numpy oder SciPy bereits haben, oder muss ich meine eigene mit numpy.linalg.eigh ?
Ich will nicht nur Einzelwertzerlegung (SVD) verwenden, da meine Eingangsdaten sind sehr hochdimensionalen (~ 460 Dimensionen), so denke ich, SVD langsamer sein wird als die Eigenvektoren der Kovarianzmatrix berechnet wird.
Ich habe gehofft, eine premade zu finden, die Umsetzung debuggt, dass bereits die richtigen Entscheidungen für wann macht die Methode zu verwenden, und was vielleicht andere Optimierungen tut, dass ich weiß nicht, über.
Lösung
Sie können einen Blick auf MDP .
Ich habe nicht hatte die Chance, es selbst zu testen, aber ich habe es vorgemerkt genau für die PCA-Funktionalität.
Andere Tipps
Einige Monate später, ist hier eine kleine Klasse PCA, und ein Bild:
#!/usr/bin/env python
""" a small class for Principal Component Analysis
Usage:
p = PCA( A, fraction=0.90 )
In:
A: an array of e.g. 1000 observations x 20 variables, 1000 rows x 20 columns
fraction: use principal components that account for e.g.
90 % of the total variance
Out:
p.U, p.d, p.Vt: from numpy.linalg.svd, A = U . d . Vt
p.dinv: 1/d or 0, see NR
p.eigen: the eigenvalues of A*A, in decreasing order (p.d**2).
eigen[j] / eigen.sum() is variable j's fraction of the total variance;
look at the first few eigen[] to see how many PCs get to 90 %, 95 % ...
p.npc: number of principal components,
e.g. 2 if the top 2 eigenvalues are >= `fraction` of the total.
It's ok to change this; methods use the current value.
Methods:
The methods of class PCA transform vectors or arrays of e.g.
20 variables, 2 principal components and 1000 observations,
using partial matrices U' d' Vt', parts of the full U d Vt:
A ~ U' . d' . Vt' where e.g.
U' is 1000 x 2
d' is diag([ d0, d1 ]), the 2 largest singular values
Vt' is 2 x 20. Dropping the primes,
d . Vt 2 principal vars = p.vars_pc( 20 vars )
U 1000 obs = p.pc_obs( 2 principal vars )
U . d . Vt 1000 obs, p.obs( 20 vars ) = pc_obs( vars_pc( vars ))
fast approximate A . vars, using the `npc` principal components
Ut 2 pcs = p.obs_pc( 1000 obs )
V . dinv 20 vars = p.pc_vars( 2 principal vars )
V . dinv . Ut 20 vars, p.vars( 1000 obs ) = pc_vars( obs_pc( obs )),
fast approximate Ainverse . obs: vars that give ~ those obs.
Notes:
PCA does not center or scale A; you usually want to first
A -= A.mean(A, axis=0)
A /= A.std(A, axis=0)
with the little class Center or the like, below.
See also:
http://en.wikipedia.org/wiki/Principal_component_analysis
http://en.wikipedia.org/wiki/Singular_value_decomposition
Press et al., Numerical Recipes (2 or 3 ed), SVD
PCA micro-tutorial
iris-pca .py .png
"""
from __future__ import division
import numpy as np
dot = np.dot
# import bz.numpyutil as nu
# dot = nu.pdot
__version__ = "2010-04-14 apr"
__author_email__ = "denis-bz-py at t-online dot de"
#...............................................................................
class PCA:
def __init__( self, A, fraction=0.90 ):
assert 0 <= fraction <= 1
# A = U . diag(d) . Vt, O( m n^2 ), lapack_lite --
self.U, self.d, self.Vt = np.linalg.svd( A, full_matrices=False )
assert np.all( self.d[:-1] >= self.d[1:] ) # sorted
self.eigen = self.d**2
self.sumvariance = np.cumsum(self.eigen)
self.sumvariance /= self.sumvariance[-1]
self.npc = np.searchsorted( self.sumvariance, fraction ) + 1
self.dinv = np.array([ 1/d if d > self.d[0] * 1e-6 else 0
for d in self.d ])
def pc( self ):
""" e.g. 1000 x 2 U[:, :npc] * d[:npc], to plot etc. """
n = self.npc
return self.U[:, :n] * self.d[:n]
# These 1-line methods may not be worth the bother;
# then use U d Vt directly --
def vars_pc( self, x ):
n = self.npc
return self.d[:n] * dot( self.Vt[:n], x.T ).T # 20 vars -> 2 principal
def pc_vars( self, p ):
n = self.npc
return dot( self.Vt[:n].T, (self.dinv[:n] * p).T ) .T # 2 PC -> 20 vars
def pc_obs( self, p ):
n = self.npc
return dot( self.U[:, :n], p.T ) # 2 principal -> 1000 obs
def obs_pc( self, obs ):
n = self.npc
return dot( self.U[:, :n].T, obs ) .T # 1000 obs -> 2 principal
def obs( self, x ):
return self.pc_obs( self.vars_pc(x) ) # 20 vars -> 2 principal -> 1000 obs
def vars( self, obs ):
return self.pc_vars( self.obs_pc(obs) ) # 1000 obs -> 2 principal -> 20 vars
class Center:
""" A -= A.mean() /= A.std(), inplace -- use A.copy() if need be
uncenter(x) == original A . x
"""
# mttiw
def __init__( self, A, axis=0, scale=True, verbose=1 ):
self.mean = A.mean(axis=axis)
if verbose:
print "Center -= A.mean:", self.mean
A -= self.mean
if scale:
std = A.std(axis=axis)
self.std = np.where( std, std, 1. )
if verbose:
print "Center /= A.std:", self.std
A /= self.std
else:
self.std = np.ones( A.shape[-1] )
self.A = A
def uncenter( self, x ):
return np.dot( self.A, x * self.std ) + np.dot( x, self.mean )
#...............................................................................
if __name__ == "__main__":
import sys
csv = "iris4.csv" # wikipedia Iris_flower_data_set
# 5.1,3.5,1.4,0.2 # ,Iris-setosa ...
N = 1000
K = 20
fraction = .90
seed = 1
exec "\n".join( sys.argv[1:] ) # N= ...
np.random.seed(seed)
np.set_printoptions( 1, threshold=100, suppress=True ) # .1f
try:
A = np.genfromtxt( csv, delimiter="," )
N, K = A.shape
except IOError:
A = np.random.normal( size=(N, K) ) # gen correlated ?
print "csv: %s N: %d K: %d fraction: %.2g" % (csv, N, K, fraction)
Center(A)
print "A:", A
print "PCA ..." ,
p = PCA( A, fraction=fraction )
print "npc:", p.npc
print "% variance:", p.sumvariance * 100
print "Vt[0], weights that give PC 0:", p.Vt[0]
print "A . Vt[0]:", dot( A, p.Vt[0] )
print "pc:", p.pc()
print "\nobs <-> pc <-> x: with fraction=1, diffs should be ~ 0"
x = np.ones(K)
# x = np.ones(( 3, K ))
print "x:", x
pc = p.vars_pc(x) # d' Vt' x
print "vars_pc(x):", pc
print "back to ~ x:", p.pc_vars(pc)
Ax = dot( A, x.T )
pcx = p.obs(x) # U' d' Vt' x
print "Ax:", Ax
print "A'x:", pcx
print "max |Ax - A'x|: %.2g" % np.linalg.norm( Ax - pcx, np.inf )
b = Ax # ~ back to original x, Ainv A x
back = p.vars(b)
print "~ back again:", back
print "max |back - x|: %.2g" % np.linalg.norm( back - x, np.inf )
# end pca.py
PCA numpy.linalg.svd Verwendung ist super einfach. Hier ist eine einfache Demo:
import numpy as np
import matplotlib.pyplot as plt
from scipy.misc import lena
# the underlying signal is a sinusoidally modulated image
img = lena()
t = np.arange(100)
time = np.sin(0.1*t)
real = time[:,np.newaxis,np.newaxis] * img[np.newaxis,...]
# we add some noise
noisy = real + np.random.randn(*real.shape)*255
# (observations, features) matrix
M = noisy.reshape(noisy.shape[0],-1)
# singular value decomposition factorises your data matrix such that:
#
# M = U*S*V.T (where '*' is matrix multiplication)
#
# * U and V are the singular matrices, containing orthogonal vectors of
# unit length in their rows and columns respectively.
#
# * S is a diagonal matrix containing the singular values of M - these
# values squared divided by the number of observations will give the
# variance explained by each PC.
#
# * if M is considered to be an (observations, features) matrix, the PCs
# themselves would correspond to the rows of S^(1/2)*V.T. if M is
# (features, observations) then the PCs would be the columns of
# U*S^(1/2).
#
# * since U and V both contain orthonormal vectors, U*V.T is equivalent
# to a whitened version of M.
U, s, Vt = np.linalg.svd(M, full_matrices=False)
V = Vt.T
# PCs are already sorted by descending order
# of the singular values (i.e. by the
# proportion of total variance they explain)
# if we use all of the PCs we can reconstruct the noisy signal perfectly
S = np.diag(s)
Mhat = np.dot(U, np.dot(S, V.T))
print "Using all PCs, MSE = %.6G" %(np.mean((M - Mhat)**2))
# if we use only the first 20 PCs the reconstruction is less accurate
Mhat2 = np.dot(U[:, :20], np.dot(S[:20, :20], V[:,:20].T))
print "Using first 20 PCs, MSE = %.6G" %(np.mean((M - Mhat2)**2))
fig, [ax1, ax2, ax3] = plt.subplots(1, 3)
ax1.imshow(img)
ax1.set_title('true image')
ax2.imshow(noisy.mean(0))
ax2.set_title('mean of noisy images')
ax3.imshow((s[0]**(1./2) * V[:,0]).reshape(img.shape))
ax3.set_title('first spatial PC')
plt.show()
Sie können sklearn verwenden:
import sklearn.decomposition as deco
import numpy as np
x = (x - np.mean(x, 0)) / np.std(x, 0) # You need to normalize your data first
pca = deco.PCA(n_components) # n_components is the components number after reduction
x_r = pca.fit(x).transform(x)
print ('explained variance (first %d components): %.2f'%(n_components, sum(pca.explained_variance_ratio_)))
matplotlib.mlab hat eine PCA Implementierung .
SVD sollte mit 460 Dimensionen funktionieren. Es dauert etwa 7 Sekunden auf meinem Atom Netbook. Die eig () -Methode nimmt mehr Zeit (wie es sollte, es wird mehr Gleitkommaoperationen) und wird fast immer ungenauer.
Wenn Sie weniger als 460 Beispiele dann, was Sie tun möchten, ist die Streumatrix diagonalisieren (x - datamean) ^ T (x - Mittelwert), vorausgesetzt, Ihre Datenpunkte Spalten sind, und dann durch (x links Multiplikation - datamean). Die Macht schneller in dem Fall, dass Sie mehr Dimensionen als Daten haben.
Sie können ganz einfach "roll" Ihre eigene mit scipy.linalg (eine Annahme vorzentrierte dataSet data):
covmat = data.dot(data.T)
evs, evmat = scipy.linalg.eig(covmat)
Dann evs sind Ihre Eigenwerte und evmat ist Ihre Projektionsmatrix.
Wenn Sie d Dimensionen behalten möchten, verwenden Sie die erste d Eigenwerte und erste d Eigenvektoren.
Da scipy.linalg die Zersetzung hat und numpy die Matrix-Multiplikationen, was sonst benötigen Sie?
Ich beende gerade das Buch Maschinelles Lernen zu lesen: eine algorithmische Perspektive . Alle Code-Beispiele in dem Buch wurden von Python geschrieben (und fast mit Numpy). Der Code-Schnipsel chatper10.2 Hauptkomponentenanalyse vielleicht ein lesenswert. Es verwendet numpy.linalg.eig.
By the way, ich denke, SVD 460 * 460 Abmessungen sehr gut umgehen kann. Ich habe einen 6500 * 6500 SVD berechnen mit numpy / scipy.linalg.svd auf einen sehr alten PC: Pentium III 733 MHz. Um ehrlich zu sein, muss das Skript viel Speicher (ca. 1.xG) und viel Zeit (ca. 30 Minuten), um das SVD Ergebnis zu erhalten.
Aber ich denke, 460 * 460 auf einem modernen PC wird kein großes Problem sein, es sei denn, u SVD muss tun, um eine große Anzahl von Zeiten.
Sie haben keine volle Einzelwertzerlegung (SVD) müssen es alle Eigenwerte und Eigenvektoren berechnet und können für große Matrizen unerschwinglich sein. scipy und seine spärlichen Modul bieten generische Funktionen lineare algrebra auf beide spärlich und dichte Matrizen arbeiten, unter denen sich die eig * Familie von Funktionen:
http://docs.scipy.org/ doc / scipy / reference / sparse.linalg.html # Matrix-Faktorisierung
Scikit-Learn ein Python PCA Implementierung die nur dichte Matrizen für jetzt unterstützen.
Timings:
In [1]: A = np.random.randn(1000, 1000)
In [2]: %timeit scipy.sparse.linalg.eigsh(A)
1 loops, best of 3: 802 ms per loop
In [3]: %timeit np.linalg.svd(A)
1 loops, best of 3: 5.91 s per loop
Hier wird eine weitere Implementierung eines PCA-Modul für Python mit numpy, scipy und C -extensions. Das Modul führt PCA entweder einen oder der SVD NIPALS (Nonlinear Iterative Partial Least Squares) Algorithmus verwendet wird, der in C umgesetzt wird.
Wenn Sie mit 3D-Vektoren arbeiten, können Sie sich bewerben SVD kurz mit der toolbelt vg . Es ist eine dünne Schicht auf der Oberseite der numpy.
import numpy as np
import vg
vg.principal_components(data)
Es gibt auch eine bequeme Alias, wenn Sie nur die erste Hauptkomponente wollen:
vg.major_axis(data)
Ich habe die Bibliothek an meinem letzten Start, wo sie von Anwendungen wie diese motiviert war. Einfache Ideen, die ausführliche oder undurchsichtige in NumPy
