From aba9321920c935565650ca45e588c9bcdfb4c514 Mon Sep 17 00:00:00 2001 From: picost Date: Wed, 3 Aug 2016 14:23:19 +0200 Subject: [PATCH 1/4] Update supervised_pca.py # In both classes Added docstrings and comments. Changed a few names for readability. #In SupervisedPCA - removing unused alpha argument - X and Y tables shapes changed according to sklearn conventional shape (n_samples, n_features) - maximum number of component enforced to be less than minimum axis size of X (ie rank(X) assuming X has maximum rank) - gamma parameter changed according to the docstring - Removed X_scale variable that was unused (scaling was not inplace. should better be done by the user using a Pipeline) # In Kernel SPCA: - correction of call to _get_kernel method in _fit method. Correction of parameters passed to pairwise_kernels in _get_kernel method (issue with kwargs metric and kernel). - correction of the call to eigsh (use of kwargs) - TO BE CHECKED: correction of transform method. Next steps: finish debugging and modify X and Y shapes in KernelSPCA Separate both classes in distinct modules. --- MachineLearningScikitLearn/supervised_pca.py | 279 +++++++++++-------- 1 file changed, 162 insertions(+), 117 deletions(-) diff --git a/MachineLearningScikitLearn/supervised_pca.py b/MachineLearningScikitLearn/supervised_pca.py index d8cd972..8ca0b7e 100644 --- a/MachineLearningScikitLearn/supervised_pca.py +++ b/MachineLearningScikitLearn/supervised_pca.py @@ -15,137 +15,171 @@ class SupervisedPCA(BaseEstimator, TransformerMixin): """Supervised Principal component analysis (SPCA) - Non-linear dimensionality reduction through the use of kernels. - Parameters ---------- n_components: int or None Number of components. If None, all non-zero components are kept. - kernel: "linear" | "poly" | "rbf" | "sigmoid" | "precomputed" Kernel. Default: "linear" - degree : int, optional Degree for poly, rbf and sigmoid kernels. Default: 3. - gamma : float, optional Kernel coefficient for rbf and poly kernels. Default: 1/n_features. - coef0 : float, optional Independent term in poly and sigmoid kernels. - - eigen_solver: string ['auto'|'dense'|'arpack'] Select eigensolver to use. If n_components is much less than the number of training samples, arpack may be more efficient than the dense eigensolver. - tol: float convergence tolerance for arpack. Default: 0 (optimal value will be chosen by arpack) - max_iter : int maximum number of iterations for arpack Default: None (optimal value will be chosen by arpack) - Attributes ---------- - `lambdas_`, `alphas_`: Eigenvalues and eigenvectors of the centered kernel matrix - - """ - - def __init__(self, n_components=None, kernel="linear", gamma=0, degree=3, - coef0=1, alpha=1.0, fit_inverse_transform=False, + + def __init__(self, n_components=None, kernel="linear", gamma=None, degree=3, + coef0=1, fit_inverse_transform=False, eigen_solver='auto', tol=0, max_iter=None): - - + + self.n_components = n_components self.kernel = kernel.lower() - self.gamma = gamma + self.gamma = None self.degree = degree self.coef0 = coef0 - self.alpha = alpha self.fit_inverse_transform = fit_inverse_transform self.eigen_solver = eigen_solver self.tol = tol self.max_iter = max_iter self.centerer = KernelCenterer() - - + def transform(self, X): """ Returns a new X, X_trans, based on previous self.fit() estimates """ - return X.dot( self.alphas_ ) - - + return(np.dot(X,self.alphas_)) + + def fit(self,X,Y): + """Returns (gen. eigenvalues, gen. eigenvectors) of matrix Q=KxKyKx. + The snapshots in X and Y are stores by rows. + + Parameters + ---------- + X: (n_samples,n_inputs) np.ndarray + Matrix of normalized input snapshots + Y: (n_samples,n_outputs) np.ndarray + Matrix of normalized output snapshots + + Returns: + -------- + lambdas: (n_components,) np.ndarray + Array of eigenvalues of XKyX^t matrix stored by decreasing + order. Only eigenvalues >= 0 are conserved. + The number of components is at most equal to the rank of X. + alphas: (n_inputs,n_components) np.ndarray + Array containing the eigenvectors associated with the eigenva + -lues in lambdas stored by columns. It is the matrix of new + basis vectors coordinates in the original basis. + Info + ---- + Solves the generalized eigenvalues problem to build the spca coeffs as + described in algo. 1 from Ghodsi 2010 paper. + """ + if Y.ndim == 1: + Y.reshape((-1,1)) + assert X.shape[0] == Y.shape[0], "X and Y don't contain the same\ + numbers of samples. Please check the matrices." self._fit(X,Y) - return - + return(self.lambdas_,self.alphas_) + def fit_transform( self, X, Y): - - + """Fits the SPCA model (solves eigenvalues problem) and returns trans + -formed input variables. + + Parameters + ---------- + X: (n_samples,n_inputs) np.ndarray + Matrix of normalized input snapshots + Y: (n_samples,n_outputs) np.ndarray + Matrix of normalized output snapshots + + Returns: + -------- + X_fit_trans: (nsamples,n_components) np.ndarray + Matrix of coordinates of input samples used for fit projected + on the space spanned by the eigenvectors computed during the + fit. The coordinates are expressed in the basis defined by the + eigenvectors. + """ self.fit( X,Y) - return self._transform() - + return(self._transform()) + def _transform(self): - - return self.X_fit.dot(self.alphas_) + """Returns the transformed version of input variables used for fit. + """ + return(np.dot(self.X_fit,self.alphas_)) + - def _fit(self, X, Y): - #find kenerl matrix of Y - K = self.centerer.fit_transform(self._get_kernel(Y)) - #scale X - X_scale = scale(X) - - + #store the input X and Y used for fit + self.X_fit = X + self._nfeature = X.shape[1] + #Compute kernel matrix of Y + Ky = self.centerer.fit_transform(self._get_kernel(Y)) + #maximum number of components that can be kept is the min between the + #number of features and samples. + n_comp_max = min(X.shape) if self.n_components is None: - n_components = K.shape[0] + n_components = n_comp_max else: - n_components = min(K.shape[0], self.n_components) - - #compute eigenvalues of X^TKX - - M = (X.T).dot(K).dot(X) - print "here" + n_components = min(n_comp_max, self.n_components) + #compute eigenvalues and eigenvetors of XTKX^T + #--------------------------------------------- + M = (X.T).dot(Ky).dot(X) + # Choosing eigensolver if self.eigen_solver == 'auto': if M.shape[0] > 200 and n_components < 10: eigen_solver = 'arpack' else: eigen_solver = 'dense' - else: - eigen_solver = self.eigen_solver - - if eigen_solver == 'dense': - self.lambdas_, self.alphas_ = linalg.eigh( - M, eigvals=(M.shape[0] - n_components, M.shape[0] - 1)) - elif eigen_solver == 'arpack': + else: + eigen_solver = self.eigen_solver + #Eigen-problem resolution + if eigen_solver == 'dense': + index_eigv_min = M.shape[0] - n_components + index_eigv_max = M.shape[0] - 1 + self.lambdas_, self.alphas_ = linalg.eigh(M, + eigvals=(index_eigv_min, + index_eigv_max)) + elif eigen_solver == 'arpack': self.lambdas_, self.alphas_ = eigsh(M, n_components, which="LA", tol=self.tol) + #----------------------------------------------------------- + #Sorting eigenvalues and eigenvectors by decreasing order of + #eigenvalues indices = self.lambdas_.argsort()[::-1] self.lambdas_ = self.lambdas_[indices] self.alphas_ = self.alphas_[:, indices] - #remove the zero/negative eigenvalues self.alphas_ = self.alphas_[:, self.lambdas_ > 0 ] self.lambdas_ = self.lambdas_[ self.lambdas_ > 0 ] - print self.alphas_.shape - - self.X_fit = X; - - + return() + + def _get_kernel(self, X, Y=None): - params = {"gamma": self.gamma, + params = {"gamma": self.gamma or 1/self._nfeature, "degree": self.degree, "coef0": self.coef0} try: @@ -155,65 +189,48 @@ def _get_kernel(self, X, Y=None): raise ValueError("%s is not a valid kernel. Valid kernels are: " "rbf, poly, sigmoid, linear and precomputed." % self.kernel) - - + -class KernelSupervisedPCA( BaseEstimator, TransformerMixin): - +class KernelSupervisedPCA(BaseEstimator, TransformerMixin): """Kernel Supervised Principal component analysis (SPCA) - Non-linear dimensionality reduction through the use of kernels. Parameters ---------- n_components: int or None Number of components. If None, all non-zero components are kept. - x||ykernel: "linear" | "poly" | "rbf" | "sigmoid" | "precomputed" Kernel. Default: "linear" - degree : int, optional Degree for poly, rbf and sigmoid kernels. Default: 3. - gamma : float, optional Kernel coefficient for rbf and poly kernels. Default: 1/n_features. - coef0 : float, optional Independent term in poly and sigmoid kernels. - - eigen_solver: string ['auto'|'dense'|'arpack'] Select eigensolver to use. If n_components is much less than the number of training samples, arpack may be more efficient than the dense eigensolver. - tol: float convergence tolerance for arpack. Default: 0 (optimal value will be chosen by arpack) - max_iter : int maximum number of iterations for arpack Default: None (optimal value will be chosen by arpack) - Attributes ---------- - `lambdas_`, `alphas_`: Eigenvalues and eigenvectors of the centered kernel matrix - - """ - + def __init__(self, n_components=None, xkernel={'kernel': "linear", 'gamma':0, 'degree':3, 'coef0':1}, ykernel = {'kernel': "linear", 'gamma':0, 'degree':3, 'coef0':1}, fit_inverse_transform=False, eigen_solver='auto', tol=0, max_iter=None): - - self.n_components = n_components self.xkernel = xkernel self.ykernel = ykernel @@ -222,76 +239,104 @@ def __init__(self, n_components=None, xkernel={'kernel': "linear", 'gamma':0, 'd self.tol = tol self.max_iter = max_iter self.centerer = KernelCenterer() - - + + def transform(self, X): """ - Returns a new X, X_trans, based on previous self.fit() estimates + Returns a new X, X_trans, based on previous self.fit() estimates. """ - K = self._get_kernel(self.X_fit, self.xkernel, X ) - return K.T.dot( self.alphas_ ) - - + K = self.centerer.fit_transform(self._get_kernel(self.X_fit, self.xkernel, X)) + X_trans = np.dot(self.alphas_.T,K) + return(X_trans) + + def fit(self,X,Y): + """Returns (gen. eigenvalues, gen. eigenvectors) of matrix Q=KxKyKx. + The snapshots in X and Y are stores by rows. + + Parameters + ---------- + X: (n_samples,n_inputs) np.ndarray + Matrix of normalized input snapshots + Y: (n_samples,n_outputs) np.ndarray + Matrix of normalized output snapshots + Info + ---- + Solves the generalized eigenvalues problem to build the kspca coeffs as + described in algo. 3 from Ghodsi 2010 paper. + """ self._fit(X,Y) - return - + return(self.lambdas_,self.alphas_) + def fit_transform( self, X, Y): - - + """Fits the KSPCA model (solves eigenvalues problem) and returns trans + -formed input variables. + """ self.fit( X,Y) - return self._transform() - + return(self._transform()) + def _transform(self): - - return self.Kx_fit.dot(self.alphas_) + """Returns the transformed version of input variables used for fit. + """ + return(self.Kx_fit.dot(self.alphas_)) - - def _fit(self, X, Y): - #find kenerl matrix of Y - Ky = self.centerer.fit_transform(self._get_kernel(Y), self.ykernel) - Kx = self.centerer.fit_transform( self._get_kernel(X), self.xkernel) - - + def _fit(self, X, Y): + #find kenerl matrices of X and Y + Ky = self.centerer.fit_transform(self._get_kernel(Y, self.ykernel)) + Kx = self.centerer.fit_transform( self._get_kernel(X, self.xkernel)) + #stores the X, Kx and Ky matrices used for fit + self.X_fit = X + self.Kx_fit = Kx + self.Ky_fit = Ky + #n_components is set as the min between the specified number of compo- + #nents and the number of samples if self.n_components is None: n_components = Ky.shape[0] else: n_components = min(Ky.shape[0], self.n_components) - - #compute eigenvalues of X^TKX - + #-------------------------------------------------------------- + #Compute generalized eigenvalues and eigenvectors of Kx^T.Ky.Kx M = (Kx).dot(Ky).dot(Kx) - if self.eigen_solver == 'auto': + + #Chose the eigensolver to be used + if self.eigen_solver == 'auto': if M.shape[0] > 200 and n_components < 10: eigen_solver = 'arpack' else: eigen_solver = 'dense' else: eigen_solver = self.eigen_solver - + #Solve the generalized eigenvalues problem if eigen_solver == 'dense': self.lambdas_, self.alphas_ = linalg.eigh( M, Kx, eigvals=(M.shape[0] - n_components, M.shape[0] - 1)) elif eigen_solver == 'arpack': - self.lambdas_, self.alphas_ = eigsh(M, Kx, n_components, + self.lambdas_, self.alphas_ = eigsh(A=M, M=Kx, k=n_components, which="LA", tol=self.tol) + else: + #useless for now + self.lambdas_, self.alphas_ = self.eigen_solver(M, Kx, n_components) + + #Sort the eigenvalues in increasing order indices = self.lambdas_.argsort()[::-1] self.lambdas_ = self.lambdas_[indices] self.alphas_ = self.alphas_[:, indices] - + #remove the zero/negative eigenvalues self.alphas_ = self.alphas_[:, self.lambdas_ > 0 ] self.lambdas_ = self.lambdas_[ self.lambdas_ > 0 ] - - self.X_fit = X; - self.Kx_fit = Kx; - + return() + def _get_kernel(self, X, params, Y=None): + """Computes and returns kernel of X with given params. + If Y is specified, then returns the pairwise kernel between X and Y. + """ try: - return pairwise_kernels(X, Y, metric=params['kernel'], - n_jobs = -1, **params) + coparams = copy.copy(params) + return pairwise_kernels(X, Y, metric=coparams.pop('kernel'), + n_jobs = 1, **coparams) except AttributeError: raise ValueError("%s is not a valid kernel. Valid kernels are: " "rbf, poly, sigmoid, linear and precomputed." From 95c93d0ee38c55e0fba2193d56253e1dc5e2c6ad Mon Sep 17 00:00:00 2001 From: picost Date: Fri, 5 Aug 2016 14:59:53 +0200 Subject: [PATCH 2/4] Update supervised_pca.py Added a parameter to enable the user to chose whether or not the negative or dropped eigenvalues should be dropped. When performing trials, it can be usefull to keep all eigenvalues such that the span of the initial and transformed vectors is the same. Also added a little more documentation about alphas_ and lambdas_ --- MachineLearningScikitLearn/supervised_pca.py | 53 ++++++++++++++------ 1 file changed, 37 insertions(+), 16 deletions(-) diff --git a/MachineLearningScikitLearn/supervised_pca.py b/MachineLearningScikitLearn/supervised_pca.py index 8ca0b7e..d31aaea 100644 --- a/MachineLearningScikitLearn/supervised_pca.py +++ b/MachineLearningScikitLearn/supervised_pca.py @@ -41,16 +41,24 @@ class SupervisedPCA(BaseEstimator, TransformerMixin): max_iter : int maximum number of iterations for arpack Default: None (optimal value will be chosen by arpack) + drop_negative_eigenvalues: bool, optional, default: False + If set to True, negative or null eigenvalues are dropped as well as the + associated eigenvectors. Attributes ---------- - `lambdas_`, `alphas_`: - Eigenvalues and eigenvectors of the centered kernel matrix + `lambdas_`: (n_components,) np.ndarray + Array of eigenvalues of XKyX^t matrix stored by decreasing + order. Only eigenvalues >= 0 are conserved. + The number of components is at most equal to the rank of X. + `alphas_`: (n_inputs,n_components) np.ndarray + Array containing the eigenvectors associated with the eigenvalues + in lambdas stored by columns. It is the matrix of new basis vectors + coordinates in the original basis. """ def __init__(self, n_components=None, kernel="linear", gamma=None, degree=3, - coef0=1, fit_inverse_transform=False, - eigen_solver='auto', tol=0, max_iter=None): - + coef0=1, fit_inverse_transform=False, eigen_solver='auto', + tol=0, max_iter=None,drop_negative_eigenvalues=False): self.n_components = n_components self.kernel = kernel.lower() @@ -62,7 +70,8 @@ def __init__(self, n_components=None, kernel="linear", gamma=None, degree=3, self.tol = tol self.max_iter = max_iter self.centerer = KernelCenterer() - + self._drop_negative_eigenvalues = drop_negative_eigenvalues + def transform(self, X): """ Returns a new X, X_trans, based on previous self.fit() estimates @@ -172,9 +181,10 @@ def _fit(self, X, Y): indices = self.lambdas_.argsort()[::-1] self.lambdas_ = self.lambdas_[indices] self.alphas_ = self.alphas_[:, indices] - #remove the zero/negative eigenvalues - self.alphas_ = self.alphas_[:, self.lambdas_ > 0 ] - self.lambdas_ = self.lambdas_[ self.lambdas_ > 0 ] + if self._drop_negative_eigenvalues: + #remove the zero/negative eigenvalues + self.alphas_ = self.alphas_[:, self.lambdas_ > 0 ] + self.lambdas_ = self.lambdas_[ self.lambdas_ > 0 ] return() @@ -221,16 +231,26 @@ class KernelSupervisedPCA(BaseEstimator, TransformerMixin): max_iter : int maximum number of iterations for arpack Default: None (optimal value will be chosen by arpack) + drop_negative_eigenvalues: bool, optional, default: False + If set to True, negative or null eigenvalues are dropped as well as the + associated eigenvectors. Attributes ---------- - `lambdas_`, `alphas_`: - Eigenvalues and eigenvectors of the centered kernel matrix + `lambdas_`: (n_components,) np.ndarray + Array of eigenvalues of XKyX^t matrix stored by decreasing + order. Only eigenvalues >= 0 are conserved. + The number of components is at most equal to the rank of X. + `alphas_`: (n_inputs,n_components) np.ndarray + Array containing the eigenvectors associated with the eigenvalues + in lambdas stored by columns. It is the matrix of new basis vectors + coordinates in the original basis. """ def __init__(self, n_components=None, xkernel={'kernel': "linear", 'gamma':0, 'degree':3, 'coef0':1}, ykernel = {'kernel': "linear", 'gamma':0, 'degree':3, 'coef0':1}, fit_inverse_transform=False, - eigen_solver='auto', tol=0, max_iter=None): + eigen_solver='auto', tol=0, max_iter=None, + drop_negative_eigenvalues=False): self.n_components = n_components self.xkernel = xkernel self.ykernel = ykernel @@ -239,7 +259,7 @@ def __init__(self, n_components=None, xkernel={'kernel': "linear", 'gamma':0, 'd self.tol = tol self.max_iter = max_iter self.centerer = KernelCenterer() - + self._drop_negative_eigenvalues = drop_negative_eigenvalues def transform(self, X): """ @@ -324,9 +344,10 @@ def _fit(self, X, Y): self.lambdas_ = self.lambdas_[indices] self.alphas_ = self.alphas_[:, indices] - #remove the zero/negative eigenvalues - self.alphas_ = self.alphas_[:, self.lambdas_ > 0 ] - self.lambdas_ = self.lambdas_[ self.lambdas_ > 0 ] + if self._drop_negative_eigenvalues: + #remove the zero/negative eigenvalues + self.alphas_ = self.alphas_[:, self.lambdas_ > 0 ] + self.lambdas_ = self.lambdas_[ self.lambdas_ > 0 ] return() def _get_kernel(self, X, params, Y=None): From 0908cf80204e05031e130d618ca0f773170eaf5d Mon Sep 17 00:00:00 2001 From: picost Date: Fri, 5 Aug 2016 15:10:43 +0200 Subject: [PATCH 3/4] Create kernel_supervised_pca Create a separate module to define KernelSupervisedPCA. --- .../kernel_supervised_pca | 177 ++++++++++++++++++ 1 file changed, 177 insertions(+) create mode 100644 MachineLearningScikitLearn/kernel_supervised_pca diff --git a/MachineLearningScikitLearn/kernel_supervised_pca b/MachineLearningScikitLearn/kernel_supervised_pca new file mode 100644 index 0000000..54e5005 --- /dev/null +++ b/MachineLearningScikitLearn/kernel_supervised_pca @@ -0,0 +1,177 @@ +# -*- coding: utf-8 -*- +""" +Supervised PCA according to Supervised Principal Compontent Anaysis by Ghodsi et al. 2010 +doi:10.1016/j.patcog.2010.12.015 +Algorithm 3 of the paper is implemented here. + +@autor: CamDavidsonPilon +@contributor: picost +""" +import numpy as np +from scipy import linalg + +from ..utils.arpack import eigsh +from ..base import BaseEstimator, TransformerMixin +from ..preprocessing import KernelCenterer, scale +from ..metrics.pairwise import pairwise_kernels + +class KernelSupervisedPCA(BaseEstimator, TransformerMixin): + """Kernel Supervised Principal component analysis (SPCA) + Non-linear dimensionality reduction through the use of kernels. + + Parameters + ---------- + n_components: int or None + Number of components. If None, all non-zero components are kept. + x||ykernel: "linear" | "poly" | "rbf" | "sigmoid" | "precomputed" + Kernel. + Default: "linear" + degree : int, optional + Degree for poly, rbf and sigmoid kernels. + Default: 3. + gamma : float, optional + Kernel coefficient for rbf and poly kernels. + Default: 1/n_features. + coef0 : float, optional + Independent term in poly and sigmoid kernels. + eigen_solver: string ['auto'|'dense'|'arpack'] + Select eigensolver to use. If n_components is much less than + the number of training samples, arpack may be more efficient + than the dense eigensolver. + tol: float + convergence tolerance for arpack. + Default: 0 (optimal value will be chosen by arpack) + max_iter : int + maximum number of iterations for arpack + Default: None (optimal value will be chosen by arpack) + drop_negative_eigenvalues: bool, optional, default: False + If set to True, negative or null eigenvalues are dropped as well as the + associated eigenvectors. + Attributes + ---------- + `lambdas_`: (n_components,) np.ndarray + Array of eigenvalues of XKyX^t matrix stored by decreasing + order. Only eigenvalues >= 0 are conserved. + The number of components is at most equal to the rank of X. + `alphas_`: (n_inputs,n_components) np.ndarray + Array containing the eigenvectors associated with the eigenvalues + in lambdas stored by columns. It is the matrix of new basis vectors + coordinates in the original basis. + """ + + def __init__(self, n_components=None, xkernel={'kernel': "linear", 'gamma':0, 'degree':3, + 'coef0':1}, ykernel = {'kernel': "linear", 'gamma':0, 'degree':3, + 'coef0':1}, fit_inverse_transform=False, + eigen_solver='auto', tol=0, max_iter=None, + drop_negative_eigenvalues=False): + self.n_components = n_components + self.xkernel = xkernel + self.ykernel = ykernel + self.fit_inverse_transform = fit_inverse_transform + self.eigen_solver = eigen_solver + self.tol = tol + self.max_iter = max_iter + self.centerer = KernelCenterer() + self._drop_negative_eigenvalues = drop_negative_eigenvalues + + def transform(self, X): + """ + Returns a new X, X_trans, based on previous self.fit() estimates. + """ + K = self.centerer.fit_transform(self._get_kernel(self.X_fit, self.xkernel, X)) + X_trans = np.dot(self.alphas_.T,K) + return(X_trans) + + + def fit(self,X,Y): + """Returns (gen. eigenvalues, gen. eigenvectors) of matrix Q=KxKyKx. + The snapshots in X and Y are stores by rows. + + Parameters + ---------- + X: (n_samples,n_inputs) np.ndarray + Matrix of normalized input snapshots + Y: (n_samples,n_outputs) np.ndarray + Matrix of normalized output snapshots + Info + ---- + Solves the generalized eigenvalues problem to build the kspca coeffs as + described in algo. 3 from Ghodsi 2010 paper. + """ + self._fit(X,Y) + return(self.lambdas_,self.alphas_) + + def fit_transform( self, X, Y): + """Fits the KSPCA model (solves eigenvalues problem) and returns trans + -formed input variables. + """ + self.fit( X,Y) + return(self._transform()) + + def _transform(self): + """Returns the transformed version of input variables used for fit. + """ + return(self.Kx_fit.dot(self.alphas_)) + + + def _fit(self, X, Y): + #find kenerl matrices of X and Y + Ky = self.centerer.fit_transform(self._get_kernel(Y, self.ykernel)) + Kx = self.centerer.fit_transform( self._get_kernel(X, self.xkernel)) + #stores the X, Kx and Ky matrices used for fit + self.X_fit = X + self.Kx_fit = Kx + self.Ky_fit = Ky + #n_components is set as the min between the specified number of compo- + #nents and the number of samples + if self.n_components is None: + n_components = Ky.shape[0] + else: + n_components = min(Ky.shape[0], self.n_components) + #-------------------------------------------------------------- + #Compute generalized eigenvalues and eigenvectors of Kx^T.Ky.Kx + M = (Kx).dot(Ky).dot(Kx) + + #Chose the eigensolver to be used + if self.eigen_solver == 'auto': + if M.shape[0] > 200 and n_components < 10: + eigen_solver = 'arpack' + else: + eigen_solver = 'dense' + else: + eigen_solver = self.eigen_solver + #Solve the generalized eigenvalues problem + if eigen_solver == 'dense': + self.lambdas_, self.alphas_ = linalg.eigh( + M, Kx, eigvals=(M.shape[0] - n_components, M.shape[0] - 1)) + elif eigen_solver == 'arpack': + self.lambdas_, self.alphas_ = eigsh(A=M, M=Kx, k=n_components, + which="LA", + tol=self.tol) + else: + #useless for now + self.lambdas_, self.alphas_ = self.eigen_solver(M, Kx, n_components) + + #Sort the eigenvalues in increasing order + indices = self.lambdas_.argsort()[::-1] + self.lambdas_ = self.lambdas_[indices] + self.alphas_ = self.alphas_[:, indices] + + if self._drop_negative_eigenvalues: + #remove the zero/negative eigenvalues + self.alphas_ = self.alphas_[:, self.lambdas_ > 0 ] + self.lambdas_ = self.lambdas_[ self.lambdas_ > 0 ] + return() + + def _get_kernel(self, X, params, Y=None): + """Computes and returns kernel of X with given params. + If Y is specified, then returns the pairwise kernel between X and Y. + """ + try: + coparams = copy.copy(params) + return pairwise_kernels(X, Y, metric=coparams.pop('kernel'), + n_jobs = 1, **coparams) + except AttributeError: + raise ValueError("%s is not a valid kernel. Valid kernels are: " + "rbf, poly, sigmoid, linear and precomputed." + % params['kernel']) From 0e1ef80d058036f544c86524a9a909938153c3a7 Mon Sep 17 00:00:00 2001 From: picost Date: Fri, 5 Aug 2016 15:12:15 +0200 Subject: [PATCH 4/4] Update supervised_pca.py Removes KernelSupervisedPCA from this module as it is implemented in a separated module now. (+ added the doi of the paper in the header) --- MachineLearningScikitLearn/supervised_pca.py | 175 +------------------ 1 file changed, 9 insertions(+), 166 deletions(-) diff --git a/MachineLearningScikitLearn/supervised_pca.py b/MachineLearningScikitLearn/supervised_pca.py index d31aaea..a4a67ee 100644 --- a/MachineLearningScikitLearn/supervised_pca.py +++ b/MachineLearningScikitLearn/supervised_pca.py @@ -1,4 +1,12 @@ -#supervised PCA according to Supervised Principal Compontent Anaysis by Ghodsi et al. 2010 +# -*- coding: utf-8 -*- +""" +Supervised PCA according to Supervised Principal Compontent Anaysis by Ghodsi et al. 2010 +doi:10.1016/j.patcog.2010.12.015 +Algorithm 1 of the paper is implemented here. + +@autor: CamDavidsonPilon +@contributor: picost +""" import numpy as np from scipy import linalg @@ -11,8 +19,6 @@ from time import clock - - class SupervisedPCA(BaseEstimator, TransformerMixin): """Supervised Principal component analysis (SPCA) Non-linear dimensionality reduction through the use of kernels. @@ -199,168 +205,5 @@ def _get_kernel(self, X, Y=None): raise ValueError("%s is not a valid kernel. Valid kernels are: " "rbf, poly, sigmoid, linear and precomputed." % self.kernel) - - - -class KernelSupervisedPCA(BaseEstimator, TransformerMixin): - """Kernel Supervised Principal component analysis (SPCA) - Non-linear dimensionality reduction through the use of kernels. - - Parameters - ---------- - n_components: int or None - Number of components. If None, all non-zero components are kept. - x||ykernel: "linear" | "poly" | "rbf" | "sigmoid" | "precomputed" - Kernel. - Default: "linear" - degree : int, optional - Degree for poly, rbf and sigmoid kernels. - Default: 3. - gamma : float, optional - Kernel coefficient for rbf and poly kernels. - Default: 1/n_features. - coef0 : float, optional - Independent term in poly and sigmoid kernels. - eigen_solver: string ['auto'|'dense'|'arpack'] - Select eigensolver to use. If n_components is much less than - the number of training samples, arpack may be more efficient - than the dense eigensolver. - tol: float - convergence tolerance for arpack. - Default: 0 (optimal value will be chosen by arpack) - max_iter : int - maximum number of iterations for arpack - Default: None (optimal value will be chosen by arpack) - drop_negative_eigenvalues: bool, optional, default: False - If set to True, negative or null eigenvalues are dropped as well as the - associated eigenvectors. - Attributes - ---------- - `lambdas_`: (n_components,) np.ndarray - Array of eigenvalues of XKyX^t matrix stored by decreasing - order. Only eigenvalues >= 0 are conserved. - The number of components is at most equal to the rank of X. - `alphas_`: (n_inputs,n_components) np.ndarray - Array containing the eigenvectors associated with the eigenvalues - in lambdas stored by columns. It is the matrix of new basis vectors - coordinates in the original basis. - """ - - def __init__(self, n_components=None, xkernel={'kernel': "linear", 'gamma':0, 'degree':3, - 'coef0':1}, ykernel = {'kernel': "linear", 'gamma':0, 'degree':3, - 'coef0':1}, fit_inverse_transform=False, - eigen_solver='auto', tol=0, max_iter=None, - drop_negative_eigenvalues=False): - self.n_components = n_components - self.xkernel = xkernel - self.ykernel = ykernel - self.fit_inverse_transform = fit_inverse_transform - self.eigen_solver = eigen_solver - self.tol = tol - self.max_iter = max_iter - self.centerer = KernelCenterer() - self._drop_negative_eigenvalues = drop_negative_eigenvalues - - def transform(self, X): - """ - Returns a new X, X_trans, based on previous self.fit() estimates. - """ - K = self.centerer.fit_transform(self._get_kernel(self.X_fit, self.xkernel, X)) - X_trans = np.dot(self.alphas_.T,K) - return(X_trans) - - - def fit(self,X,Y): - """Returns (gen. eigenvalues, gen. eigenvectors) of matrix Q=KxKyKx. - The snapshots in X and Y are stores by rows. - - Parameters - ---------- - X: (n_samples,n_inputs) np.ndarray - Matrix of normalized input snapshots - Y: (n_samples,n_outputs) np.ndarray - Matrix of normalized output snapshots - Info - ---- - Solves the generalized eigenvalues problem to build the kspca coeffs as - described in algo. 3 from Ghodsi 2010 paper. - """ - self._fit(X,Y) - return(self.lambdas_,self.alphas_) - - def fit_transform( self, X, Y): - """Fits the KSPCA model (solves eigenvalues problem) and returns trans - -formed input variables. - """ - self.fit( X,Y) - return(self._transform()) - - def _transform(self): - """Returns the transformed version of input variables used for fit. - """ - return(self.Kx_fit.dot(self.alphas_)) - - - def _fit(self, X, Y): - #find kenerl matrices of X and Y - Ky = self.centerer.fit_transform(self._get_kernel(Y, self.ykernel)) - Kx = self.centerer.fit_transform( self._get_kernel(X, self.xkernel)) - #stores the X, Kx and Ky matrices used for fit - self.X_fit = X - self.Kx_fit = Kx - self.Ky_fit = Ky - #n_components is set as the min between the specified number of compo- - #nents and the number of samples - if self.n_components is None: - n_components = Ky.shape[0] - else: - n_components = min(Ky.shape[0], self.n_components) - #-------------------------------------------------------------- - #Compute generalized eigenvalues and eigenvectors of Kx^T.Ky.Kx - M = (Kx).dot(Ky).dot(Kx) - - #Chose the eigensolver to be used - if self.eigen_solver == 'auto': - if M.shape[0] > 200 and n_components < 10: - eigen_solver = 'arpack' - else: - eigen_solver = 'dense' - else: - eigen_solver = self.eigen_solver - #Solve the generalized eigenvalues problem - if eigen_solver == 'dense': - self.lambdas_, self.alphas_ = linalg.eigh( - M, Kx, eigvals=(M.shape[0] - n_components, M.shape[0] - 1)) - elif eigen_solver == 'arpack': - self.lambdas_, self.alphas_ = eigsh(A=M, M=Kx, k=n_components, - which="LA", - tol=self.tol) - else: - #useless for now - self.lambdas_, self.alphas_ = self.eigen_solver(M, Kx, n_components) - - #Sort the eigenvalues in increasing order - indices = self.lambdas_.argsort()[::-1] - self.lambdas_ = self.lambdas_[indices] - self.alphas_ = self.alphas_[:, indices] - - if self._drop_negative_eigenvalues: - #remove the zero/negative eigenvalues - self.alphas_ = self.alphas_[:, self.lambdas_ > 0 ] - self.lambdas_ = self.lambdas_[ self.lambdas_ > 0 ] - return() - - def _get_kernel(self, X, params, Y=None): - """Computes and returns kernel of X with given params. - If Y is specified, then returns the pairwise kernel between X and Y. - """ - try: - coparams = copy.copy(params) - return pairwise_kernels(X, Y, metric=coparams.pop('kernel'), - n_jobs = 1, **coparams) - except AttributeError: - raise ValueError("%s is not a valid kernel. Valid kernels are: " - "rbf, poly, sigmoid, linear and precomputed." - % params['kernel'])