initial commit

6 files changed, 408 insertions, 0 deletions

diff --git a/R/cd.sparsestep.R b/R/cd.sparsestep.R
new file mode 100644
index 0000000..29dbec8
--- /dev/null
+++ b/R/cd.sparsestep.R
@@ -0,0 +1,73 @@
+#' @title Coordinate descent algorithm for SparseStep
+#'
+#'
+#' @export
+#'
+sparsestep.cd <- function(x, y, lambdas=NULL, epsilon=1e-5, intercept=TRUE,
+			  ...)
+{
+	nm <- dim(X)
+	n <- nm[1]
+	m <- nm[2]
+	one <- rep(1, n)
+
+	if (intercept) {
+		meanx <- drop(one %*% x)/n
+		x <- scale(x, meanx, FALSE)
+		mu <- mean(y)
+		y <- drop(y - mu)
+	} else {
+		meanx <- rep(0, m)
+		mu <- 0
+		y <- drop(y)
+	}
+
+	XX <- t(x) %*% x
+	Xy <- t(x) %*% y
+
+	#gammas <- 2^(seq(log(1e6)/log(2), log(1e-8)/log(2)))
+
+	num.lambdas <- length(lambdas)
+	#num.gammas <- length(gammas)
+	
+	#betas <- array(0, dim=c(num.gammas, num.lambdas, m))
+	betas <- array(0, dim=c(num.lambdas, m))
+	for (l in num.lambdas:1) {
+		lambda <- lambdas[l]
+		# initialize beta
+		if (l == num.lambdas) {
+			beta <- as.vector(matrix(0, 1, m))
+		} else {
+			#beta <- betas[num.gammas, l+1, ]
+			beta <- betas[l+1, ]
+		}
+
+		j <- 1
+		last.beta <- as.vector(matrix(0, 1, m))
+		while (TRUE) {
+			# code
+			b <- -2 * x[, j] %*% (y - x[, -j] %*% last.beta[-j])
+			a <- x[, j] %*% x[, j]
+			if (abs(last.beta[j]) > epsilon) {
+				beta[j] <- b/a
+			} else {
+				beta[j] <- (2*b - lambda*sign(last.beta[j]))/
+					(2*a)
+			}
+			# check convergence
+			if (sum(abs(beta - last.beta)) < 1e-10) {
+				break
+			} else {
+				last.beta <- beta
+			}
+			# continue
+			j <- j %% m + 1
+		}
+		betas[l, ] <- beta
+	}
+
+	return(betas)
+}
+
+
+
diff --git a/R/fit.sparsestep.R b/R/fit.sparsestep.R
new file mode 100644
index 0000000..ca40f43
--- /dev/null
+++ b/R/fit.sparsestep.R
@@ -0,0 +1,139 @@
+#' @title Fits the SparseStep model
+#'
+#' @description Fits the SparseStep model for a single value of the 
+#' regularization parameter.
+#'
+#' @param x matrix of predictors
+#' @param y response
+#' @param lambda regularization parameter
+#' @param gamma0 starting value of the gamma parameter
+#' @param gammastop stopping value of the gamma parameter
+#' @param IMsteps number of steps of the majorization algorithm to perform for 
+#' each value of gamma
+#' @param gammastep factor to decrease gamma with at each step
+#' @param normalize if TRUE, each variable is standardized to have unit L2 
+#' norm, otherwise it is left alone.
+#' @param intercept if TRUE, an intercept is included in the model (and not 
+#' penalized), otherwise no intercept is included
+#' @param force.zero if TRUE, absolute coefficients smaller than the provided 
+#' threshold value are set to absolute zero as a post-processing step, 
+#' otherwise no thresholding is performed
+#' @param threshold threshold value to use for setting coefficients to 
+#' absolute zero
+#' @param XX The X'X matrix; useful for repeated runs where X'X stays the same
+#' @param Xy The X'y matrix; useful for repeated runs where X'y stays the same
+#' @param use.XX whether or not to compute X'X and return it
+#' @param use.Xy whether or not to compute X'y and return it
+#'
+#' @return A "sparsestep" object is returned, for which predict, coef, methods 
+#' exist.
+#'
+#' @export
+#'
+#' @examples
+#' data(diabetes)
+#' attach(diabetes)
+#' object <- sparsestep(x, y)
+#' plot(object)
+#' detach(diabetes)
+#'
+sparsestep <- function(x, y, lambda=1.0, gamma0=1e6, 
+		       gammastop=1e-8, IMsteps=2, gammastep=2.0, normalize=TRUE,
+		       intercept=TRUE, force.zero=TRUE, threshold=1e-7,
+		       XX=NULL, Xy=NULL, use.XX = TRUE, use.Xy = TRUE)
+{
+	call <- match.call()
+
+	nm <- dim(x)
+	n <- nm[1]
+	m <- nm[2]
+	one <- rep(1, n)
+
+	if (intercept) {
+		meanx <- drop(one %*% x)/n
+		x <- scale(x, meanx, FALSE)
+		mu <- mean(y)
+		y <- drop(y - mu)
+	} else {
+		meanx <- rep(0, m)
+		mu <- 0
+		y <- drop(y)
+	}
+
+	if (normalize) {
+		normx <- sqrt(drop(one %*% (x^2)))
+		names(normx) <- NULL
+		x <- scale(x, FALSE, normx)
+		cat("Normalizing in sparsestep\n")
+	} else {
+		normx <- rep(1, m)
+	}
+
+	if (use.XX & is.null(XX)) {
+		XX <- t(x) %*% x
+	}
+	if (use.Xy & is.null(Xy)) {
+		Xy <- t(x) %*% y
+	}
+
+	# Start solving SparseStep	
+	gamma <- gamma0
+	beta <- matrix(0.0, m, 1)
+
+	it <- 0
+	while (gamma > gammastop)
+       	{
+		for (i in 1:IMsteps) {
+			alpha <- beta
+			omega <- gamma/(alpha^2 + gamma)^2
+			Omega <- diag(as.vector(omega), m, m)
+			beta <- solve(XX + lambda * Omega, Xy)
+		}
+		it <- it + 1
+		gamma <- gamma / gammastep
+	}
+
+	# perform IM until convergence
+	epsilon <- 1e-14
+	loss <- get.loss(x, y, gamma, beta, lambda)
+	lbar <- (1.0 + 2.0 * epsilon) * loss
+	while ((lbar - loss)/loss > epsilon)
+	{
+		alpha <- beta
+		omega <- gamma/(alpha^2 + gamma)^2
+		Omega <- diag(as.vector(omega), m, m)
+		beta <- solve(XX + lambda * Omega, Xy)
+		lbar <- loss
+		loss <- get.loss(x, y, gamma, beta, lambda)
+	}
+
+	# postprocessing
+	if (force.zero) {
+		beta[which(abs(beta) < threshold)] <- 0
+	}
+	residuals <- y - x %*% beta
+	beta <- scale(t(beta), FALSE, normx)
+	RSS <- apply(residuals^2, 2, sum)
+	R2 <- 1 - RSS/RSS[1]
+
+	object <- list(call = call, lambda = lambda, R2 = R2, RSS = RSS,
+		       gamma0 = gamma0, gammastop = gammastop,
+		       IMsteps = IMsteps, gammastep = gammastep,
+		       intercept = intercept, force.zero = force.zero,
+		       threshold = threshold, beta = beta, mu = mu,
+		       normx = normx, meanx = meanx,
+		       XX = if(use.XX) XX else NULL,
+		       Xy = if(use.Xy) Xy else NULL)
+	class(object) <- "sparsestep"
+	return(object)
+}
+
+get.loss <- function(x, y, gamma, beta, lambda)
+{
+	Xb <- x %*% beta
+	diff <- y - Xb
+	b2 <- beta^2
+	binv <- 1/(b2 + gamma)
+	loss <- t(diff) %*% diff + lambda * t(b2) %*% binv
+	return(loss)
+}
diff --git a/R/path.sparsestep.R b/R/path.sparsestep.R
new file mode 100644
index 0000000..40467b8
--- /dev/null
+++ b/R/path.sparsestep.R
@@ -0,0 +1,163 @@
+#' @title Path algorithm for the SparseStep model
+#'
+#' @description Fits the entire regularization path for SparseStep using a 
+#' Golden Section search.
+#'
+#' @param x matrix of predictors
+#' @param y response
+#' @param max.depth maximum recursion depth
+#' @param ... further arguments to sparsestep()
+#'
+#' @export
+#'
+sparsestep.path <- function(x, y, max.depth=10, intercept=TRUE,
+			    normalize=TRUE, ...)
+{
+	call <- match.call()
+
+	nm <- dim(x)
+	n <- nm[1]
+	m <- nm[2]
+	one <- rep(1, n)
+
+	if (intercept) {
+		meanx <- drop(one %*% x)/n
+		x <- scale(x, meanx, FALSE)
+		mu <- mean(y)
+		y <- drop(y - mu)
+	} else {
+		meanx <- rep(0, m)
+		mu <- 0
+		y <- drop(y)
+	}
+
+	if (normalize) {
+		normx <- sqrt(drop(one %*% (x^2)))
+		names(normx) <- NULL
+		x <- scale(x, FALSE, normx)
+		cat("Normalizing in path.sparsestep\n")
+	} else {
+		normx <- rep(1, m)
+	}
+
+	XX <- t(x) %*% x
+	Xy <- t(x) %*% y
+
+	iter <- 0
+
+	# First find the smallest lambda for which all coefficients are zero
+	lambda.max <- 2^25
+	fit <- NULL
+	while (1) {
+		last.fit <- fit
+		fit <- sparsestep(x, y, lambda=lambda.max, normalize=FALSE, 
+				  XX=XX, Xy=Xy, ...)
+		iter <- iter + 1
+		if (all(fit$beta == 0)) {
+			lambda.max <- lambda.max / 2
+		} else {
+			lambda.max <- lambda.max * 2
+			break
+		}
+	}
+	cat("Found maximum value of lambda: 2^(", log(lambda.max)/log(2), ")\n")
+	iter <- iter + 1
+	if (is.null(last.fit)) {
+		betas.max <- fit$beta
+	} else {
+		betas.max <- last.fit$beta
+	}
+
+	# Now find the largest lambda for which no coefficients are zero
+	lambda.min <- 2^-12
+	fit <- NULL
+	while (1) {
+		last.fit <- fit
+		fit <- sparsestep(x, y, lambda=lambda.min, normalize=FALSE, 
+				  XX=XX, Xy=Xy, ...)
+		iter <- iter + 1
+		if (all(fit$beta != 0)) {
+			lambda.min <- lambda.min * 2
+		} else {
+			lambda.min <- lambda.min / 2
+			break
+		}
+	}
+	cat("Found minimum value of lambda: 2^(", log(lambda.min)/log(2), ")\n")
+	iter <- iter + 1
+	if (is.null(last.fit)) {
+		betas.min <- fit$beta
+	} else {
+		betas.min <- last.fit$beta
+	}
+
+	# Run binary section search
+	have.zeros <- as.vector(matrix(FALSE, 1, m+1))
+	have.zeros[1] <- TRUE
+	have.zeros[m+1] <- TRUE
+
+	left <- log(lambda.min)/log(2)
+	right <- log(lambda.max)/log(2)
+
+	l <- lambda.search(x, y, 0, max.depth, have.zeros, left, right, 1, 
+			   m+1, XX, Xy, ...)
+	have.zeros <- have.zeros | l$have.zeros
+	lambdas <- c(lambda.min, l$lambdas, lambda.max)
+	betas <- rbind(betas.min, l$betas, betas.max)
+	iter <- iter + l$iter
+
+	ord <- order(lambdas)
+	lambdas <- lambdas[ord]
+	betas <- betas[ord, ]
+	betas <- scale(betas, FALSE, normx)
+
+	object <- list(call=call, lambdas=lambdas, betas=betas,
+		       iterations=iter)
+}
+
+lambda.search <- function(x, y, depth, max.depth, have.zeros, left, right,
+			  lidx, ridx, XX, Xy, ...)
+{
+	cat("Running search in interval [", left, ",", right, "] ... \n")
+	nm <- dim(x)
+	m <- nm[2]
+
+	betas <- NULL
+	lambdas <- NULL
+
+	middle <- left + (right - left)/2
+	lambda <- 2^middle
+	fit <- sparsestep(x, y, lambda=lambda, normalize=FALSE, XX=XX, Xy=Xy, 
+			  ...)
+	iter <- 1
+
+	num.zero <- length(which(fit$beta == 0))
+	cidx <- num.zero + 1
+	if (have.zeros[cidx] == FALSE) {
+		have.zeros[cidx] = TRUE
+		betas <- rbind(betas, as.vector(fit$beta))
+		lambdas <- c(lambdas, lambda)
+	}
+
+	idx <- rbind(c(lidx, cidx), c(cidx, ridx))
+	bnd <- rbind(c(left, middle), c(middle, right))
+	for (r in 1:2) {
+		i1 <- idx[r, 1]
+		i2 <- idx[r, 2]
+		b1 <- bnd[r, 1]
+		b2 <- bnd[r, 2]
+		if (depth < max.depth && any(have.zeros[i1:i2] == F)) {
+			ds <- lambda.search(x, y, depth+1, max.depth,
+					    have.zeros, b1, b2, i1, i2, XX, Xy,
+					    ...)
+			have.zeros <- have.zeros | ds$have.zeros
+			lambdas <- c(lambdas, ds$lambdas)
+			betas <- rbind(betas, ds$betas)
+			iter <- iter + ds$iter
+		}
+	}
+
+	out <- list(have.zeros=have.zeros, lambdas=lambdas, betas=betas,
+		    iter=iter)
+	return(out)
+}
diff --git a/R/predict.sparsestep.R b/R/predict.sparsestep.R
new file mode 100644
index 0000000..5f8898e
--- /dev/null
+++ b/R/predict.sparsestep.R
@@ -0,0 +1,4 @@
+predict.sparsestep <- function(object, newx)
+{
+	NULL
+}
diff --git a/R/print.sparsestep.R b/R/print.sparsestep.R
new file mode 100644
index 0000000..1de1370
--- /dev/null
+++ b/R/print.sparsestep.R
@@ -0,0 +1,24 @@
+#' @title Print the fitted SparseStep model
+#'
+#' @description Prints a short text of a fitted SparseStep model
+#'
+#' @param obj a "sparsestep" object to print
+#'
+#' @export
+#'
+#' @examples
+#' data(diabetes)
+#' attach(diabetes)
+#' object <- sparsestep(x, y)
+#' print(object)
+#' detach(diabetes)
+#'
+print.sparsestep <- function(obj, ...)
+{
+	cat("\nCall:\n")
+	dput(obj$call)
+	cat("R-squared:", format(round(rev(obj$R2)[1], 3)), "\n")
+	zeros <- length(which(obj$beta == 0))/length(obj$beta)*100.0
+	cat("Percentage coeff zero:", format(round(zeros, 2)), "\n")
+	invisible(obj)
+}
diff --git a/R/sparsestep.R b/R/sparsestep.R
new file mode 100644
index 0000000..05cf9ff
--- /dev/null
+++ b/R/sparsestep.R
@@ -0,0 +1,5 @@
+#' sparsestep.
+#'
+#' @name sparsestep
+#' @docType package
+NULL