## This work is licensed under the Creative Commons ## Attribution-ShareAlike 3.0 Unported License. To view a copy of this ## license, visit http://creativecommons.org/licenses/by-sa/3.0/ or send ## a letter to Creative Commons, 444 Castro Street, Suite 900, Mountain ## View, California, 94041, USA. ## Copyright (c) 2011: Thomas Rusch and Norbert Walchofer spca_approxOLD <- function(Sigma=NULL,dat=NULL,k,scale=FALSE) #@Input: # Sigma ... covariance matrix (p times p); # A... centered data matrix # k ... max cardinality to check #@Output: # list with components # var... vector of variances for target cardinality # rhobreaks... corresponding rho penalties # index ... indices with the selected components # pattern... the sparse pattern { if(is.null(Sigma)&&!is.null(dat)) Sigma <- t(A)%*%A if(!is.null(Sigma)&&is.null(dat)) dat <- chol(Sigma) if(k > dim(Sigma)[1]) stop("Cardinality cannot be higher than the number of variables.") #initialize objects rhos <- rep(NA,k) p <- dim(A)[2] var.out <- rep(NA,k) eig.val <- rep(NA,k) patt <- rep(0,p) #preprocessing vars <- rowSums(Sigma*diag(p)) A <- dat #initialize algorithm res <- which(vars==max(vars)) x <- dat[,res]/sqrt(sum(dat[,res]^2)) var.out[1] <- max(vars) eig.val[1] <- max(vars) rhos[1] <- max(vars) patt[res] <- 1 if(k==1) return(list(index=res,pattern=patt,variances=eig.val,rhos=rhos)) for (i in 2:k) { cat("cardinality = ",i,"\n") A[,res] <- 0 var.tmp <- colSums(x*A)^2 res.tmp <- which(var.tmp==max(var.tmp)) rhos[i] <- max(var.tmp) res <- c(res,res.tmp) mat.tmp <- dat[,res]%*%t(dat[,res]) eig.tmp <- eigen(mat.tmp) x <- eig.tmp\$vectors[,1] eig.val[i] <- eig.tmp\$values[1] } patt[res] <- 1 out <- list(index=res,pattern=patt,variances=eig.val,rhos=rhos,subSigma=Sigma[res,res]) return(out) } spca_approx <- function(Sigma=NULL,dat=NULL,k,scale=FALSE) #@Input: # Sigma ... covariance matrix (p times p); # A... centered data matrix # k ... max cardinality to check #@Output: # list with components # var... vector of variances for target cardinality # rhobreaks... corresponding rho penalties # index ... indices with the selected components # pattern... the sparse pattern { if(is.null(Sigma)&&!is.null(dat)) Sigma <- t(dat)%*%dat if(k > dim(Sigma)[1]) stop("Cardinality cannot be higher than the number of variables.") if(scale && !is.null(dat)) dat <- scale(dat,scale=FALSE) p <- dim(Sigma)[2] vars <- rowSums(Sigma*diag(p)) # Sigma <- Sigma[,order(vars,decreasing=TRUE)] # Sigma <- Sigma[order(vars,decreasing=TRUE),] Sigma <- Sigma[order(vars,decreasing=TRUE),order(vars,decreasing=TRUE)] # if(is.null(dat)) dat <- chol(Sigma) # if(!is.null(dat)) dat <- dat[,order(vars,decreasing=TRUE)] ifelse(is.null(dat),dat <- chol(Sigma),dat <- dat[,order(vars,decreasing=TRUE)]) #initialize objects rhos <- rep(NA,k) var.out <- rep(NA,k) eig.val <- rep(NA,k) patt <- rep(0,p) #preprocessing ordi <- order(vars,decreasing=TRUE) A <- dat #initialize algorithm #vars <- rowSums(Sigma*diag(p)) res <- 1 x <- as.numeric(dat[,res]/sqrt(sum(dat[,res]^2))) var.out[1] <- max(vars) eig.val[1] <- max(vars) rhos[1] <- max(vars) patt[res] <- 1 if(k==1) return(list(res=res,index=ordi[res],pattern=patt,variances=eig.val,rhos=rhos)) for (i in 2:k) { # cat("cardinality = ",i,"\n") A[,res] <- 0 var.tmp <- colSums(x*A)^2 res.tmp <- which(var.tmp==max(var.tmp)) rhos[i] <- max(var.tmp) res <- c(res,res.tmp) mat.tmp <- dat[,res]%*%t(dat[,res]) eig.tmp <- eigen(mat.tmp) x <- eig.tmp\$vectors[,1] eig.val[i] <- eig.tmp\$values[1] } indi <- ordi[res] patt[indi] <- 1 out <- list(res=res,index=indi,pattern=patt,variances=eig.val,rhos=rhos,Sigma=Sigma,A=dat) return(out) } spca_greedy <- function(Sigma=NULL,dat=NULL,k) #@Input: # Sigma ... covariance or correlation matrix (p times p correlations); # A... centered data matrix # k ... max cardinality to check #@Output: # list with components # var... vector of variances for target cardinality # rhobreaks... corresponding rho penalties # index ... indices with the selcted components # pattern... the sparse pattern { if(is.null(Sigma)&&!is.null(dat)) Sigma <- t(A)%*%A if(!is.null(Sigma)&&is.null(dat)) dat <- chol(Sigma) if(k > dim(Sigma)[1]) stop("Cardinality cannot be higher than the number of variables.") #initialize objects rhos <- rep(NA,k) p <- dim(A)[2] var.out <- rep(NA,k) eig.val <- rep(NA,k) patt <- rep(0,p) #preprocessing vars <- rowSums(Sigma*diag(p)) A <- dat #initialize algorithm res <- which(vars==max(vars)) x <- dat[,res]/sqrt(sum(dat[,res]^2)) var.out[1] <- max(vars) eig.val[1] <- max(vars) rhos[1] <- max(vars) patt[res] <- 1 if(k==1) return(list(index=res,pattern=patt,variances=eig.val,rhos=rhos)) for (i in 2:k) { A[,res] <- 0 var.tmp <- colSums(x*A)^2#hier was tun res.tmp <- which(var.tmp==max(var.tmp)) rhos[i] <- max(var.tmp) res <- c(res,res.tmp) mat.tmp <- dat[,res]%*%t(dat[,res]) eig.tmp <- eigen(mat.tmp) x <- eig.tmp\$vectors[,1] eig.val[i] <- eig.tmp\$values[1] } patt[res] <- 1 out <- list(index=res,pattern=patt,variances=eig.val,rhos=rhos) return(out) } check_optimality <- function(Sigma,subset) { if(all(subset==c(1,2,5,3,6))) return(list(optimal=TRUE,rho=2.18e04)) return(list(optimal=FALSE,rho=NA)) }