Analiza głównych komponentów w Pythonie

112

Chciałbym użyć analizy głównych składowych (PCA) do redukcji wymiarowości. Czy Numpy lub Scipy już to ma, czy muszę używać własnego numpy.linalg.eigh?

Nie chcę tylko używać rozkładu według wartości osobliwych (SVD), ponieważ moje dane wejściowe są dość wielowymiarowe (~ 460 wymiarów), więc myślę, że SVD będzie wolniejsze niż obliczanie wektorów własnych macierzy kowariancji.

Miałem nadzieję, że uda mi się znaleźć gotową, zdebugowaną implementację, która już podejmuje właściwe decyzje dotyczące tego, kiedy użyć której metody i która może zapewniać inne optymalizacje, o których nie wiem.

Vebjorn Ljosa
źródło

Odpowiedzi:

28

Możesz rzucić okiem na MDP .

Nie miałem okazji tego osobiście przetestować, ale dodałem zakładkę dokładnie pod kątem funkcjonalności PCA.

ChristopheD
źródło
8
MDP nie jest utrzymywany od 2012 roku, nie wygląda na najlepsze rozwiązanie.
Marc Garcia,
Najnowsza aktualizacja pochodzi z 09.03.2016, ale pamiętaj, że ir jest tylko wydaniem poprawiającym błędy: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

Kilka miesięcy później, oto mały PCA klasy i zdjęcie:

#!/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

wprowadź opis obrazu tutaj

denis
źródło
3
Fyinfo, jest doskonała rozmowa C. Caramanisa o Robust PCA , styczeń 2011.
denis
czy ten kod wyświetli ten obraz (Iris PCA)? Jeśli nie, czy możesz opublikować alternatywne rozwiązanie, w którym wyjście byłoby z tego obrazu. Komunikatory mają problemy z konwersją tego kodu do c ++, ponieważ jestem nowy w Pythonie :)
Orvyl
44

Korzystanie z PCA numpy.linalg.svdjest bardzo łatwe. Oto proste 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
źródło
2
Zdaję sobie sprawę, że jestem trochę spóźniony, ale OP specjalnie zażądał rozwiązania, które pozwala uniknąć rozkładu według wartości osobliwych.
Alex A.
1
@Alex Zdaję sobie z tego sprawę, ale jestem przekonany, że SVD to nadal właściwe podejście. Powinien być wystarczająco szybki dla potrzeb OP (mój przykład powyżej, z wymiarami 262144 zajmuje tylko ~ 7,5 sekundy na normalnym laptopie) i jest znacznie bardziej stabilny numerycznie niż metoda eigendecomposition (patrz komentarz dwf poniżej). Zwracam również uwagę, że zaakceptowana odpowiedź również wykorzystuje SVD!
ali_m
Nie zaprzeczam, że SVD jest drogą do zrobienia, powiedziałem tylko, że odpowiedź nie odnosi się do pytania tak, jak zostało zadane. Ale to dobra odpowiedź, niezła robota.
Alex A.
5
@Alex Fair wystarczająco. Myślę, że jest to kolejny wariant problemu XY - OP powiedział, że nie chce rozwiązania opartego na SVD, ponieważ uważał , że SVD będzie zbyt wolny, prawdopodobnie jeszcze go nie wypróbował. W takich przypadkach osobiście uważam, że bardziej pomocne jest wyjaśnienie, w jaki sposób można rozwiązać szerszy problem, niż odpowiadanie na pytanie dokładnie w jego pierwotnej, węższej formie.
ali_m
svdjuż zwraca sposortowane w porządku malejącym, o ile dokumentacja idzie. (Być może tak nie było w 2012 r., Ale dziś jest)
Etienne Bruines
34

Możesz użyć sklearn:

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
źródło
Uznano za pozytywnymi, ponieważ to dobrze działa dla mnie - mam ponad 460 wymiarów i chociaż sklearn używa SVD, a pytanie żądane nie-SVD, myślę, że 460 wymiarów prawdopodobnie będzie OK.
Dan Stowell
Możesz również chcieć usunąć kolumny o stałej wartości (std = 0). W tym celu należy użyć: remove_cols = np.where (np.all (x == np.mean (x, 0), 0)) [0] A następnie x = np.delete (x, remove_cols, 1)
Noam Peled
14

SVD powinno działać dobrze z 460 wymiarami. Na moim netbooku Atom zajmuje to około 7 sekund. Metoda eig () przyjmuje więcej czasu (tak jak powinna, wykorzystuje więcej operacji zmiennoprzecinkowych) i prawie zawsze będzie mniej dokładna.

Jeśli masz mniej niż 460 przykładów, to co chcesz zrobić, to przekątna macierz rozproszenia (x - datamean) ^ T (x - średnia), zakładając, że twoje punkty danych są kolumnami, a następnie mnożąc w lewo przez (x - datamean). Że może być szybciej w przypadku, gdy masz więcej niż wymiary danych.

dwf
źródło
czy możesz opisać bardziej szczegółowo tę sztuczkę, gdy masz więcej wymiarów niż danych?
mrgloom
1
Zasadniczo zakładasz, że wektory własne są liniowymi kombinacjami wektorów danych. Zobacz Sirovich (1987). „Turbulencja i dynamika spójnych struktur”.
dwf
11

Możesz dość łatwo „przetoczyć” własne, używając scipy.linalg(zakładając wstępnie wyśrodkowany zbiór danych data):

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

Następnie evssą twoje wartości własne ievmat jest twoją macierzą projekcji.

Jeśli chcesz zachować dwymiary, użyj pierwszych dwartości własnych i pierwszejd wektorów własnych.

Biorąc pod uwagę, że scipy.linalgma rozkład i zdrętwienie mnożenia macierzy, czego jeszcze potrzebujesz?

Has QUIT - Anony-Mousse
źródło
macierz cov to np.dot (data.T, data, out = covmat), gdzie dane muszą być wyśrodkowane.
mrgloom
2
Powinieneś przyjrzeć się komentarzowi @ dwf na temat tej odpowiedzi, aby sprawdzić niebezpieczeństwa związane z używaniem eig()na macierzy kowariancji.
Alex A.
8

Właśnie kończę czytać książkę Machine Learning: An Algorithmic Perspective . Wszystkie przykłady kodu w książce zostały napisane przez Pythona (i prawie z Numpy). Warto przeczytać fragment kodu programu chatper10.2 Analiza głównych komponentów . Używa numpy.linalg.eig.
Nawiasem mówiąc, myślę, że SVD bardzo dobrze radzi sobie z wymiarami 460 * 460. Obliczyłem 6500 * 6500 SVD z numpy / scipy.linalg.svd na bardzo starym komputerze: Pentium III 733 MHz. Szczerze mówiąc, skrypt potrzebuje dużo pamięci (około 1.xG) i dużo czasu (około 30 minut), aby uzyskać wynik SVD. Ale myślę, że 460 * 460 na nowoczesnym komputerze nie będzie dużym problemem, chyba że będziesz potrzebować SVD wiele razy.

sunqiang
źródło
28
Nigdy nie powinieneś używać eig () na macierzy kowariancji, kiedy możesz po prostu użyć svd (). W zależności od tego, ile komponentów planujesz użyć i od rozmiaru macierzy danych, błąd numeryczny wprowadzony przez ten pierwszy (wykonuje więcej operacji zmiennoprzecinkowych) może stać się znaczący. Z tego samego powodu nie powinieneś nigdy jawnie odwracać macierzy za pomocą inv (), jeśli to, co naprawdę cię interesuje, to odwrotność razy wektor lub macierz; zamiast tego powinieneś użyćolve ().
dwf
5

Nie potrzebujesz pełnej dekompozycji na wartości osobliwe (SVD), ponieważ oblicza ona wszystkie wartości własne i wektory własne i może być przeszkodą w przypadku dużych macierzy. scipy i jego rzadki moduł zapewniają ogólne liniowe funkcje algrebry działające zarówno na rzadkich, jak i gęstych macierzach, wśród których znajduje się rodzina funkcji eig *:

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

Scikit-learn zapewnia implementację Python PCA, która na razie obsługuje tylko gęste macierze.

Czasy:

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
źródło
1
Niezbyt uczciwe porównanie, ponieważ nadal musisz obliczyć macierz kowariancji. Również prawdopodobnie warto używać rzadkiego materiału linalg tylko dla bardzo dużych macierzy, ponieważ wydaje się, że tworzenie rzadkich macierzy z gęstych macierzy wydaje się być dość powolne. na przykład eigshjest faktycznie ~ 4x wolniejszy niż eighdla nierzadkich macierzy. To samo dotyczy scipy.sparse.linalg.svdsversus numpy.linalg.svd. Zawsze brałbym SVD zamiast rozkładu wartości własnych z powodów, o których wspomniał @dwf, i być może użyłbym rzadkiej wersji SVD, jeśli macierze stają się naprawdę duże.
ali_m
2
Nie musisz obliczać rzadkich macierzy z gęstych macierzy. Algorytmy zawarte w module sparse.linalg opierają się jedynie na operacji mnożenia wektorów macierzy za pomocą metody matvec obiektu Operator. W przypadku gęstych macierzy jest to coś w rodzaju matvec = dot (A, x). Z tego samego powodu nie musisz obliczać macierzy kowariancji, ale tylko podać kropkę operacyjną (AT, kropka (A, x)) dla A.
Nicolas Barbey
Ach, teraz widzę, że względna prędkość metod rzadkich i nierzadkich zależy od rozmiaru macierzy. Jeśli użyję twojego przykładu, w którym A jest matrycą 1000 * 1000, to eigshi svdssą szybsze niż eighi svdo współczynnik ~ 3, ale jeśli A jest mniejsze, powiedzmy 100 * 100, to eighi svdsą szybsze odpowiednio o współczynniki ~ 4 i ~ 1,5 . Jednak T nadal używałby rzadkiego SVD zamiast rzadkiego rozkładu wartości własnej.
ali_m
2
Rzeczywiście, myślę, że jestem nastawiony na duże matryce. Dla mnie duże matryce są bardziej jak 10⁶ * 10⁶ niż 1000 * 1000. W takim przypadku często nie można nawet przechowywać macierzy kowariancji ...
Nicolas Barbey
4

Oto kolejna implementacja modułu PCA dla Pythona przy użyciu rozszerzeń numpy, scipy i C. Moduł przeprowadza PCA przy użyciu algorytmu SVD lub NIPALS (Nonlinear Iterative Partial Least Squares), który jest zaimplementowany w C.

rcs
źródło
0

Jeśli pracujesz z wektorami 3D, możesz zwięźle zastosować SVD za pomocą paska narzędzi vg . To lekka warstwa na wierzchu numpy.

import numpy as np
import vg

vg.principal_components(data)

Istnieje również wygodny alias, jeśli chcesz tylko pierwszy główny składnik:

vg.major_axis(data)

Bibliotekę stworzyłem podczas mojego ostatniego uruchomienia, gdzie była motywowana takimi zastosowaniami: prostymi pomysłami, które są pełne lub nieprzejrzyste w NumPy.

paulmelnikow
źródło