R/kernelEVD.R

kernelEVD <-
function(x){
	if(!is.matrix(x)) {
		stop("The function kernelEVD requires input of type 'matrix'.")
	}
	n <- nrow(x)
	p <- ncol(x)
	if(n > p) {
		return(classSVD(x))
	}
	else {
		centerofX <- apply(x, 2, mean)
		Xcentered <- scale(x, center=TRUE, scale=FALSE)
		if(n == 1) {
			stop("The sample size is 1. No singular value decomposition can be performed.")
		}
		eigen <- eigen(Xcentered %*% t(Xcentered)/(n-1))
		eigen.descending <- greatsort(eigen$values)
		loadings <- eigen$vectors[,eigen.descending$index]
		tolerance <- n * max(eigen$values) * .Machine$double.eps
		rank <- sum(eigen.descending$sortedvector > tolerance)
		eigenvalues <- eigen.descending$sortedvector[1:rank]
		loadings <- t((Xcentered/sqrt(n-1))) %*% loadings[,1:rank] %*% diag(1/sqrt(eigenvalues), nrow=length(eigenvalues), ncol=length(eigenvalues))
		scores <- Xcentered %*% loadings
		return(list(loadings=as.matrix(loadings), scores=as.matrix(scores), eigenvalues=as.vector(eigenvalues), rank=rank,
						Xcentered=as.matrix(Xcentered), centerofX=as.vector(centerofX)))
	}
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.