Analisis komponen utama dengan Python

112

Saya ingin menggunakan analisis komponen utama (PCA) untuk reduksi dimensi. Apakah numpy atau scipy sudah memilikinya, atau apakah saya harus menggulung sendiri menggunakan numpy.linalg.eigh?

Saya tidak hanya ingin menggunakan dekomposisi nilai tunggal (SVD) karena data masukan saya cukup berdimensi tinggi (~ 460 dimensi), jadi menurut saya SVD akan lebih lambat daripada menghitung vektor eigen dari matriks kovarian.

Saya berharap untuk menemukan implementasi yang sudah dibuat sebelumnya dan di-debug yang sudah membuat keputusan yang tepat tentang kapan harus menggunakan metode mana, dan yang mungkin melakukan pengoptimalan lain yang tidak saya ketahui.

Vebjorn Ljosa
sumber

Jawaban:

28

Anda mungkin telah melihat MDP .

Saya belum sempat mengujinya sendiri, tetapi saya telah menandainya dengan tepat untuk fungsionalitas PCA.

ChristopheD
sumber
8
MDP belum dipertahankan sejak 2012, sepertinya bukan solusi terbaik.
Marc Garcia
Pembaruan terbaru adalah dari 09.03.2016, tetapi perhatikan bahwa ir hanya rilis perbaikan bug: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

Beberapa bulan kemudian, inilah PCA kelas kecil, dan gambarnya:

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

masukkan deskripsi gambar di sini

denis
sumber
3
Fyinfo, ada perbincangan yang sangat bagus tentang Robust PCA oleh C. Caramanis, Januari 2011.
denis
apakah kode ini akan mengeluarkan gambar itu (Iris PCA)? Jika tidak, dapatkah Anda memposting solusi alternatif di mana gambar itu akan keluar. Saya mengalami beberapa kesulitan dalam mengonversi kode ini ke c ++ karena saya baru dalam python :)
Orvyl
44

Penggunaan PCA numpy.linalg.svdsangat mudah. Berikut demo sederhana:

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
sumber
2
Saya menyadari saya agak terlambat di sini, tetapi OP secara khusus meminta solusi yang menghindari dekomposisi nilai tunggal.
Alex A.
1
@ Alex Saya menyadari hal itu, tetapi saya yakin bahwa SVD masih merupakan pendekatan yang tepat. Ini harus dengan mudah cukup cepat untuk kebutuhan OP (contoh saya di atas, dengan dimensi 262144 hanya membutuhkan ~ 7,5 detik pada laptop normal), dan itu jauh lebih stabil secara numerik daripada metode dekomposisi eigend (lihat komentar dwf di bawah). Saya juga mencatat bahwa jawaban yang diterima menggunakan SVD juga!
ali_m
Saya tidak setuju bahwa SVD adalah cara yang harus ditempuh, saya hanya mengatakan bahwa jawabannya tidak menjawab pertanyaan seperti yang dinyatakan oleh pertanyaan. Ini adalah jawaban yang bagus, kerja bagus.
Alex A.
5
@Alex Cukup adil. Saya pikir ini adalah varian lain dari masalah XY - OP mengatakan dia tidak menginginkan solusi berbasis SVD karena menurutnya SVD akan terlalu lambat, mungkin belum mencobanya. Dalam kasus seperti ini, menurut saya pribadi, akan lebih membantu untuk menjelaskan bagaimana Anda akan menangani masalah yang lebih luas, daripada menjawab pertanyaan dengan tepat dalam bentuk aslinya yang lebih sempit.
ali_m
svdsudah kembali sseperti yang diurutkan dalam urutan menurun, sejauh dokumentasi berjalan. (Mungkin ini tidak terjadi pada tahun 2012, tetapi hari ini adalah)
Etienne Bruines
34

Anda dapat menggunakan 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
sumber
Diberikan suara positif karena ini berfungsi dengan baik untuk saya - saya memiliki lebih dari 460 dimensi, dan meskipun sklearn menggunakan SVD dan pertanyaan yang diminta non-SVD, saya pikir 460 dimensi kemungkinan besar baik-baik saja.
Dan Stowell
Anda mungkin juga ingin menghapus kolom dengan nilai konstan (std = 0). Untuk itu Anda harus menggunakan: remove_cols = np.where (np.all (x == np.mean (x, 0), 0)) [0] Dan kemudian x = np.delete (x, remove_cols, 1)
Noam Peled
14

SVD harus bekerja dengan baik dengan 460 dimensi. Dibutuhkan sekitar 7 detik pada netbook Atom saya. Metode eig () membutuhkan lebih banyak waktu (sebagaimana mestinya, menggunakan lebih banyak operasi floating point) dan hampir selalu kurang akurat.

Jika Anda memiliki kurang dari 460 contoh maka yang ingin Anda lakukan adalah mendiagonalisasi matriks pencar (x - datamean) ^ T (x - mean), dengan asumsi titik data Anda adalah kolom, lalu mengalikan kiri dengan (x - datamean). Itu mungkin lebih cepat jika Anda memiliki lebih banyak dimensi daripada data.

dwf
sumber
Bisakah Anda menjelaskan lebih detail trik ini ketika Anda memiliki lebih banyak dimensi daripada data?
mrgloom
1
Pada dasarnya Anda mengasumsikan bahwa vektor eigen adalah kombinasi linier dari vektor data. Lihat Sirovich (1987). "Turbulensi dan dinamika struktur yang koheren."
dwf
11

Anda dapat dengan mudah "menggulung" sendiri menggunakan scipy.linalg(dengan asumsi kumpulan data yang dipusatkan sebelumnya data):

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

Kemudian evsadalah nilai eigen Anda, dan evmatmatriks proyeksi Anda.

Jika Anda ingin mempertahankan ddimensi, gunakan nilai deigen pertama dan dvektor eigen pertama .

Mengingat scipy.linalgperkalian matriks yang memiliki dekomposisi dan numpy, apa lagi yang Anda butuhkan?

Memiliki QUIT - Anony-Mousse
sumber
cov matrix adalah np.dot (data.T, data, out = covmat), dimana data harus berpusat pada matriks.
mrgloom
2
Anda harus melihat komentar @ dwf tentang jawaban ini untuk bahaya penggunaan eig()matriks kovarian.
Alex A.
8

Saya baru saja selesai membaca buku Machine Learning: An Algorithmic Perspective . Semua contoh kode di buku itu ditulis oleh Python (dan hampir dengan Numpy). Potongan kode chatper10.2 Analisis Komponen Utama mungkin layak dibaca. Ini menggunakan numpy.linalg.eig.
Omong-omong, saya rasa SVD dapat menangani dimensi 460 * 460 dengan sangat baik. Saya telah menghitung 6500 * 6500 SVD dengan numpy / scipy.linalg.svd pada PC yang sangat tua: Pentium III 733mHz. Sejujurnya, skrip membutuhkan banyak memori (sekitar 1.xG) dan banyak waktu (sekitar 30 menit) untuk mendapatkan hasil SVD. Tapi saya pikir 460 * 460 pada PC modern tidak akan menjadi masalah besar kecuali jika Anda perlu melakukan SVD berkali-kali.

sunqiang
sumber
28
Anda tidak boleh menggunakan eig () pada matriks kovarians jika Anda cukup menggunakan svd (). Bergantung pada berapa banyak komponen yang Anda rencanakan untuk digunakan dan ukuran matriks data Anda, kesalahan numerik yang diperkenalkan oleh yang pertama (ini melakukan lebih banyak operasi floating point) bisa menjadi signifikan. Untuk alasan yang sama, Anda tidak boleh secara eksplisit membalik matriks dengan inv () jika yang Anda minati adalah waktu inversi vektor atau matriks; Anda harus menggunakan Solving () sebagai gantinya.
dwf
5

Anda tidak perlu Dekomposisi Nilai Singular (SVD) penuh karena ia menghitung semua nilai eigen dan vektor eigen dan dapat menjadi penghalang untuk matriks besar. scipy dan modul renggangnya menyediakan fungsi algrebra linier generik yang bekerja pada matriks renggang dan padat, di antaranya terdapat keluarga fungsi eig *:

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

Scikit-learn menyediakan implementasi Python PCA yang hanya mendukung matriks padat untuk saat ini.

Waktu:

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
sumber
1
Bukan perbandingan yang adil, karena Anda masih perlu menghitung matriks kovarians. Juga mungkin hanya layak menggunakan barang linalg renggang untuk matriks yang sangat besar, karena tampaknya cukup lambat untuk membuat matriks renggang dari matriks padat. misalnya, eigshsebenarnya ~ 4x lebih lambat daripada eighmatriks nonsparse. Hal yang sama berlaku untuk scipy.sparse.linalg.svdsversus numpy.linalg.svd. Saya akan selalu menggunakan SVD daripada dekomposisi eigenvalue karena alasan yang disebutkan @dwf, dan mungkin menggunakan versi SVD yang jarang jika matriks menjadi sangat besar.
ali_m
2
Anda tidak perlu menghitung matriks renggang dari matriks padat. Algoritme yang disediakan dalam modul sparse.linalg hanya mengandalkan operasi perkalian vektor matriks melalui metode matvec pada objek Operator. Untuk matriks padat, ini hanya seperti matvec = dot (A, x). Untuk alasan yang sama, Anda tidak perlu menghitung matriks kovarian tetapi hanya memberikan titik operasi (AT, titik (A, x)) untuk A.
Nicolas Barbey
Ah, sekarang saya melihat bahwa kecepatan relatif metode sparse vs nonsparse bergantung pada ukuran matriks. Jika saya menggunakan contoh Anda di mana A adalah matriks 1000 * 1000 lalu eigshdan svdslebih cepat daripada eighdan svddengan faktor ~ 3, tetapi jika A lebih kecil, katakan 100 * 100, maka eighdan svdlebih cepat dengan faktor ~ 4 dan ~ 1,5 masing-masing . T masih akan menggunakan SVD jarang pada dekomposisi nilai eigen yang jarang.
ali_m
2
Memang, saya pikir saya bias terhadap matriks besar. Bagi saya, matriks besar lebih seperti 10⁶ * 10⁶ daripada 1000 * 1000. Dalam kasus tersebut, Anda bahkan sering tidak dapat menyimpan matriks kovarian ...
Nicolas Barbey
4

Berikut adalah implementasi lain dari modul PCA untuk python menggunakan numpy, scipy dan C-extension. Modul ini menjalankan PCA menggunakan algoritma SVD atau NIPALS (Nonlinear Iterative Partial Least Squares) yang diimplementasikan di C.

rcs
sumber
0

Jika Anda bekerja dengan vektor 3D, Anda dapat menerapkan SVD secara ringkas menggunakan toolbelt vg . Ini adalah lapisan tipis di atas numpy.

import numpy as np
import vg

vg.principal_components(data)

Ada juga alias yang nyaman jika Anda hanya menginginkan komponen utama pertama:

vg.major_axis(data)

Saya membuat perpustakaan ini pada startup terakhir saya, yang dimotivasi oleh penggunaan seperti ini: ide sederhana yang bertele-tele atau tidak jelas di NumPy.

paulmelnikow.dll
sumber