Hauptkomponentenanalyse in Python

112

Ich möchte die Hauptkomponentenanalyse (PCA) zur Reduzierung der Dimensionalität verwenden. Hat numpy oder scipy es schon oder muss ich mein eigenes mit rollennumpy.linalg.eigh ?

Ich möchte nicht nur die Singular Value Decomposition (SVD) verwenden, da meine Eingabedaten ziemlich hochdimensional sind (~ 460 Dimensionen), daher denke ich, dass SVD langsamer ist als die Berechnung der Eigenvektoren der Kovarianzmatrix.

Ich hatte gehofft, eine vorgefertigte, debuggte Implementierung zu finden, die bereits die richtigen Entscheidungen trifft, wann welche Methode verwendet werden soll und welche möglicherweise andere Optimierungen vornimmt, von denen ich nichts weiß.

Vebjorn Ljosa
quelle

Antworten:

28

Vielleicht werfen Sie einen Blick auf MDP .

Ich hatte noch keine Gelegenheit, es selbst zu testen, aber ich habe es genau für die PCA-Funktionalität mit einem Lesezeichen versehen.

ChristopheD
quelle
8
MDP wurde seit 2012 nicht mehr gewartet und scheint nicht die beste Lösung zu sein.
Marc Garcia
Das neueste Update ist vom 09.03.2016, aber beachten Sie, dass ir nur eine Bugfix-Version ist:Note that from this release MDP is in maintenance mode. 13 years after its first public release, MDP has reached full maturity and no new features are planned in the future.
Gabriel
65

Monate später 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

Geben Sie hier die Bildbeschreibung ein

denis
quelle
3
Fyinfo, es gibt einen ausgezeichneten Vortrag über Robust PCA von C. Caramanis, Januar 2011.
Denis
Gibt dieser Code dieses Bild aus (Iris PCA)? Wenn nicht, können Sie eine alternative Lösung veröffentlichen, bei der das Out das Bild wäre. Ich habe einige Schwierigkeiten bei der Konvertierung dieses Codes in C ++, weil ich neu in Python bin :)
Orvyl
44

Die Verwendung von PCA numpy.linalg.svdist 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()
ali_m
quelle
2
Mir ist klar, dass ich hier etwas spät dran bin, aber das OP hat speziell eine Lösung angefordert, die eine Zerlegung einzelner Werte vermeidet .
Alex A.
1
@Alex Ich weiß das, aber ich bin überzeugt, dass SVD immer noch der richtige Ansatz ist. Es sollte schnell genug für die Anforderungen des OP sein (mein Beispiel oben mit 262144-Abmessungen dauert auf einem normalen Laptop nur ~ 7,5 Sekunden) und es ist numerisch viel stabiler als die Eigendecomposition-Methode (siehe den Kommentar von dwf unten). Ich stelle auch fest, dass die akzeptierte Antwort auch SVD verwendet!
Ali_m
Ich bin nicht anderer Meinung, dass SVD der richtige Weg ist. Ich habe nur gesagt, dass die Antwort die Frage nicht so behandelt, wie sie gestellt wird. Es ist jedoch eine gute Antwort, gute Arbeit.
Alex A.
5
@ Alex Fair genug. Ich denke, dies ist eine weitere Variante des XY-Problems - das OP sagte, er wolle keine SVD-basierte Lösung, weil er dachte, SVD werde zu langsam sein, wahrscheinlich ohne es noch versucht zu haben. In solchen Fällen halte ich es persönlich für hilfreicher zu erklären, wie Sie das umfassendere Problem angehen würden, als die Frage genau in ihrer ursprünglichen, engeren Form zu beantworten.
Ali_m
svdkehrt bereits sin absteigender Reihenfolge sortiert zurück, soweit die Dokumentation reicht. (Vielleicht war dies 2012 nicht der Fall, aber heute ist es so)
Etienne Bruines
34

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_)))
Noam Peled
quelle
Upvoted, weil dies für mich gut funktioniert - ich habe mehr als 460 Dimensionen, und obwohl sklearn SVD und die angeforderte Frage ohne SVD verwendet, denke ich, dass 460 Dimensionen wahrscheinlich in Ordnung sind.
Dan Stowell
Möglicherweise möchten Sie auch Spalten mit einem konstanten Wert (std = 0) entfernen. Dazu sollten Sie Folgendes verwenden: remove_cols = np.where (np.all (x == np.mean (x, 0), 0)) [0] Und dann x = np.delete (x, remove_cols, 1)
Noam Peled
14

SVD sollte mit 460 Dimensionen gut funktionieren. Auf meinem Atom-Netbook dauert es ungefähr 7 Sekunden. Die eig () -Methode benötigt mehr Zeit (wie es sollte, verwendet sie mehr Gleitkommaoperationen) und ist fast immer weniger genau.

Wenn Sie weniger als 460 Beispiele haben, möchten Sie die Streumatrix (x - Datamean) ^ T (x - Mittelwert) diagonalisieren, vorausgesetzt, Ihre Datenpunkte sind Spalten, und dann mit (x - Datamean) links multiplizieren. Dies ist möglicherweise schneller, wenn Sie mehr Dimensionen als Daten haben.

dwf
quelle
Können Sie diesen Trick genauer beschreiben, wenn Sie mehr Dimensionen als Daten haben?
Mrgloom
1
Grundsätzlich nehmen Sie an, dass die Eigenvektoren lineare Kombinationen der Datenvektoren sind. Siehe Sirovich (1987). "Turbulenzen und die Dynamik kohärenter Strukturen."
Dwf
11

Sie können ganz einfach Ihre eigenen "rollen" scipy.linalg(unter der Annahme eines vorzentrierten Datensatzes data):

covmat = data.dot(data.T)
evs, evmat = scipy.linalg.eig(covmat)

Dann evssind Ihre Eigenwerte und evmatIhre Projektionsmatrix.

Wenn Sie dDimensionen beibehalten möchten , verwenden Sie die ersten dEigenwerte und ersten dEigenvektoren.

scipy.linalgWas brauchen Sie angesichts der Zerlegung und der Anzahl der Matrixmultiplikationen noch?

Hat aufgehört - Anony-Mousse
quelle
Die cov-Matrix ist np.dot (data.T, data, out = covmat), wobei die Daten eine zentrierte Matrix sein müssen.
Mrgloom
2
Sie sollten sich den Kommentar von @ dwf zu dieser Antwort ansehen, um die Gefahren einer Verwendung eig()in einer Kovarianzmatrix zu ermitteln.
Alex A.
8

Ich lese gerade das Buch Maschinelles Lernen: Eine algorithmische Perspektive . Alle Codebeispiele im Buch wurden von Python (und fast mit Numpy) geschrieben. Das Code-Snippet von chatper10.2 Principal Components Analysis ist möglicherweise eine Lektüre wert. Es wird numpy.linalg.eig verwendet.
Ich denke übrigens, SVD kann 460 * 460 Dimensionen sehr gut verarbeiten. Ich habe eine 6500 * 6500 SVD mit numpy / scipy.linalg.svd auf einem sehr alten PC berechnet: Pentium III 733mHz. Um ehrlich zu sein, benötigt 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, Sie müssen SVD sehr oft ausführen.

Sunqiang
quelle
28
Sie sollten niemals eig () für eine Kovarianzmatrix verwenden, wenn Sie einfach svd () verwenden können. Abhängig von der Anzahl der Komponenten, die Sie verwenden möchten, und der Größe Ihrer Datenmatrix kann der durch die erstere verursachte numerische Fehler (es werden mehr Gleitkommaoperationen ausgeführt) erheblich werden. Aus dem gleichen Grund sollten Sie eine Matrix niemals explizit mit inv () invertieren, wenn Sie wirklich an den inversen Zeiten eines Vektors oder einer Matrix interessiert sind. Sie sollten stattdessen Solve () verwenden.
Dwf
5

Sie benötigen keine vollständige Singular Value Decomposition (SVD), da sie alle Eigenwerte und Eigenvektoren berechnet und für große Matrizen unzulässig sein kann. scipy und sein Sparse-Modul bieten generische lineare Algrebra-Funktionen, die sowohl für Sparse- als auch für dichte Matrizen arbeiten. Darunter befindet sich die eig * -Funktionsfamilie:

http://docs.scipy.org/doc/scipy/reference/sparse.linalg.html#matrix-factorizations

Scikit-learn bietet eine Python-PCA-Implementierung, die derzeit nur dichte Matrizen unterstützt.

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
Nicolas Barbey
quelle
1
Kein wirklich fairer Vergleich, da Sie die Kovarianzmatrix noch berechnen müssen. Außerdem lohnt es sich wahrscheinlich nur, das spärliche Linalg-Material für sehr große Matrizen zu verwenden, da es ziemlich langsam zu sein scheint, spärliche Matrizen aus dichten Matrizen zu konstruieren. Zum Beispiel eigshist es tatsächlich ~ 4x langsamer als eighbei nicht sparsamen Matrizen. Gleiches gilt für scipy.sparse.linalg.svdsversus numpy.linalg.svd. Ich würde aus den von @dwf genannten Gründen immer mit SVD über die Eigenwertzerlegung gehen und vielleicht eine spärliche Version von SVD verwenden, wenn die Matrizen wirklich riesig werden.
Ali_m
2
Sie müssen keine spärlichen Matrizen aus dichten Matrizen berechnen. Die im Modul sparse.linalg bereitgestellten Algorithmen basieren nur auf der Matrizenvektormultiplikationsoperation über die matvec-Methode des Operator-Objekts. Für dichte Matrizen ist dies nur so etwas wie matvec = Punkt (A, x). Aus dem gleichen Grund müssen Sie nicht die Kovarianzmatrix berechnen, sondern nur den Operationspunkt (AT, Punkt (A, x)) für A.
Nicolas Barbey
Ah, jetzt sehe ich, dass die relative Geschwindigkeit der spärlichen gegenüber der nicht sparsamen Methode von der Größe der Matrix abhängt. Wenn ich Ihr Beispiel verwende, in dem A eine 1000 * 1000-Matrix ist eigshund svdsdann schneller als eighund svdum einen Faktor von ~ 3 ist, aber wenn A kleiner ist, sagen wir 100 * 100, dann eighund svdum Faktoren von ~ 4 bzw. ~ 1,5 schneller sind . T würde jedoch immer noch eine spärliche SVD über eine spärliche Eigenwertzerlegung verwenden.
Ali_m
2
In der Tat denke ich, ich bin voreingenommen gegenüber großen Matrizen. Für mich sind große Matrizen eher 10⁶ * 10⁶ als 1000 * 1000. In diesen Fällen können Sie oft nicht einmal die Kovarianzmatrizen speichern ...
Nicolas Barbey
4

Hier ist eine weitere Implementierung eines PCA-Moduls für Python mit Numpy-, Scipy- und C-Erweiterungen. Das Modul führt die PCA entweder mit einem SVD- oder dem in C implementierten NIPALS-Algorithmus (Nonlinear Iterative Partial Least Squares) durch.

rcs
quelle
0

Wenn Sie mit 3D-Vektoren arbeiten, können Sie SVD mit dem Toolbelt vg präzise anwenden . Es ist eine leichte Schicht auf Numpy.

import numpy as np
import vg

vg.principal_components(data)

Es gibt auch einen praktischen Alias, wenn Sie nur die erste Hauptkomponente möchten:

vg.major_axis(data)

Ich habe die Bibliothek bei meinem letzten Start erstellt, wo sie durch Verwendungen wie diese motiviert war: einfache Ideen, die in NumPy ausführlich oder undurchsichtig sind.

paulmelnikow
quelle