initial commit

14 files changed, 551 insertions, 0 deletions

diff --git a/DESCRIPTION b/DESCRIPTION
new file mode 100644
index 0000000..ae36cf2
--- /dev/null
+++ b/DESCRIPTION
@@ -0,0 +1,11 @@
+Package: sparsestep
+Type: Package
+Title: SparseStep Regression
+Version: 0.1
+Authors@R: person("Gertjan", "van den Burg", email="sparsestep@gmail.com",
+		  role=c("aut", "cre"))
+Description: Implements the SparseStep regression model.
+Depends: R (>= 2.10)
+License: GPL-2
+URL: http://sparsestep.github.com
+LazyData: true
diff --git a/INDEX b/INDEX
new file mode 100644
index 0000000..e69de29
--- /dev/null
+++ b/INDEX
diff --git a/MD5 b/MD5
new file mode 100644
index 0000000..e69de29
--- /dev/null
+++ b/MD5
diff --git a/NAMESPACE b/NAMESPACE
new file mode 100644
index 0000000..0e14978
--- /dev/null
+++ b/NAMESPACE
@@ -0,0 +1,6 @@
+# Generated by roxygen2 (4.1.1): do not edit by hand
+
+S3method(print,sparsestep)
+export(sparsestep)
+export(sparsestep.cd)
+export(sparsestep.path)
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
diff --git a/man/print.sparsestep.Rd b/man/print.sparsestep.Rd
new file mode 100644
index 0000000..6680788
--- /dev/null
+++ b/man/print.sparsestep.Rd
@@ -0,0 +1,22 @@
+% Generated by roxygen2 (4.1.1): do not edit by hand
+% Please edit documentation in R/print.sparsestep.R
+\name{print.sparsestep}
+\alias{print.sparsestep}
+\title{Print the fitted SparseStep model}
+\usage{
+\method{print}{sparsestep}(obj, ...)
+}
+\arguments{
+\item{obj}{a "sparsestep" object to print}
+}
+\description{
+Prints a short text of a fitted SparseStep model
+}
+\examples{
+data(diabetes)
+attach(diabetes)
+object <- sparsestep(x, y)
+print(object)
+detach(diabetes)
+}
+
diff --git a/man/sparsestep.Rd b/man/sparsestep.Rd
new file mode 100644
index 0000000..3d4a311
--- /dev/null
+++ b/man/sparsestep.Rd
@@ -0,0 +1,68 @@
+% Generated by roxygen2 (4.1.1): do not edit by hand
+% Please edit documentation in R/fit.sparsestep.R, R/sparsestep.R
+\docType{package}
+\name{sparsestep}
+\alias{sparsestep}
+\alias{sparsestep-package}
+\title{Fits the SparseStep model}
+\usage{
+sparsestep(x, y, lambda = 1, gamma0 = 1e+06, gammastop = 1e-08,
+  IMsteps = 2, gammastep = 2, normalize = TRUE, intercept = TRUE,
+  force.zero = TRUE, threshold = 1e-07, XX = NULL, Xy = NULL,
+  use.XX = TRUE, use.Xy = TRUE)
+}
+\arguments{
+\item{x}{matrix of predictors}
+
+\item{y}{response}
+
+\item{lambda}{regularization parameter}
+
+\item{gamma0}{starting value of the gamma parameter}
+
+\item{gammastop}{stopping value of the gamma parameter}
+
+\item{IMsteps}{number of steps of the majorization algorithm to perform for
+each value of gamma}
+
+\item{gammastep}{factor to decrease gamma with at each step}
+
+\item{normalize}{if TRUE, each variable is standardized to have unit L2
+norm, otherwise it is left alone.}
+
+\item{intercept}{if TRUE, an intercept is included in the model (and not
+penalized), otherwise no intercept is included}
+
+\item{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}
+
+\item{threshold}{threshold value to use for setting coefficients to
+absolute zero}
+
+\item{XX}{The X'X matrix; useful for repeated runs where X'X stays the same}
+
+\item{Xy}{The X'y matrix; useful for repeated runs where X'y stays the same}
+
+\item{use.XX}{whether or not to compute X'X and return it}
+
+\item{use.Xy}{whether or not to compute X'y and return it}
+}
+\value{
+A "sparsestep" object is returned, for which predict, coef, methods
+exist.
+}
+\description{
+Fits the SparseStep model for a single value of the
+regularization parameter.
+
+sparsestep.
+}
+\examples{
+data(diabetes)
+attach(diabetes)
+object <- sparsestep(x, y)
+plot(object)
+detach(diabetes)
+}
+
diff --git a/man/sparsestep.cd.Rd b/man/sparsestep.cd.Rd
new file mode 100644
index 0000000..94cf48a
--- /dev/null
+++ b/man/sparsestep.cd.Rd
@@ -0,0 +1,13 @@
+% Generated by roxygen2 (4.1.1): do not edit by hand
+% Please edit documentation in R/cd.sparsestep.R
+\name{sparsestep.cd}
+\alias{sparsestep.cd}
+\title{Coordinate descent algorithm for SparseStep}
+\usage{
+sparsestep.cd(x, y, lambdas = NULL, epsilon = 1e-05, intercept = TRUE,
+  ...)
+}
+\description{
+Coordinate descent algorithm for SparseStep
+}
+
diff --git a/man/sparsestep.path.Rd b/man/sparsestep.path.Rd
new file mode 100644
index 0000000..43625d5
--- /dev/null
+++ b/man/sparsestep.path.Rd
@@ -0,0 +1,23 @@
+% Generated by roxygen2 (4.1.1): do not edit by hand
+% Please edit documentation in R/path.sparsestep.R
+\name{sparsestep.path}
+\alias{sparsestep.path}
+\title{Path algorithm for the SparseStep model}
+\usage{
+sparsestep.path(x, y, max.depth = 10, intercept = TRUE, normalize = TRUE,
+  ...)
+}
+\arguments{
+\item{x}{matrix of predictors}
+
+\item{y}{response}
+
+\item{max.depth}{maximum recursion depth}
+
+\item{...}{further arguments to sparsestep()}
+}
+\description{
+Fits the entire regularization path for SparseStep using a
+Golden Section search.
+}
+